Landsat NIR band and ELM-FATES sensitivity to forest disturbances and regrowth in the Central Amazon

Forest disturbance and regrowth are key processes in forest dynamics but detailed information of these processes is 15 difficult to obtain in remote forests as the Amazon. We used chronosequences of Landsat satellite imagery to determine the sensitivity of surface reflectance from all spectral bands to windthrow, clearcutting, and burning and their successional pathways of forest regrowth in the Central Amazon. We also assess whether the forest demography model Functionally Assembled Terrestrial Ecosystem Simulator (FATES) implemented in the Energy Exascale Earth System Model (E3SM) Land Model (ELM), ELM-FATES, accurately represents the changes for windthrow and clearcut. The results show that all spectral 20 bands from Landsat satellite were sensitive to the disturbances but after 3 to 6 years only the Near Infrared (NIR) band had significant changes associated with the successional pathways of forest regrowth for all the disturbances considered. In general, the NIR decreased immediately after disturbance, increased to maximum values with the establishment of pioneers and earlysuccessional tree species, and then decreased slowly and almost linearly to pre-disturbance conditions with the dynamics of forest succession. Statistical methods predict that NIR will return to pre-disturbance values in about 39 years (consistent with 25 observational data of biomass regrowth following windthrows), and 36 and 56 years for clearcut and burning. The NIR captured the observed successional pathways of forest regrowth after clearcut and burning that diverge through time. ELMFATES predicted higher peaks of initial forest responses (e.g., biomass, stem density) after clearcuts than after windthrows, similar to the changes in NIR. However, ELM-FATES predicted a faster recovery of forest structure and canopy-coverage back to pre-disturbance conditions for windthrows compared to clearcuts. The similarity of ELM-FATES predictions of 30 regrowth patterns after windthrow and clearcut to those of the NIR results suggest that the dynamics of forest regrowth for these disturbances are represented with appropriate fidelity within ELM-FATES and useful as a benchmarking tool. https://doi.org/10.5194/bg-2019-451 Preprint. Discussion started: 10 December 2019 c © Author(s) 2019. CC BY 4.0 License.

Forests in the Central Amazon affected by windthrow, clearcut, and cut+burn were addressed in this study. Windthrows (Mitchell, 2013) in the Amazon are caused by strong descending winds that uproot or break trees (Garstang et al., 1998;Negrón-Juárez et al., 2018;. In clearcut areas, forests are cut and cleared and in cut+burn areas forest are cleared and burned (Mesquita et al., 2001;Mesquita et al., 2015;Lovejoy and Bierregaard, 1990). The windthrow, clearcut, and cut+burn sites used in this study were selected based on the following conditions: (a) prior to disturbance they were upland 105 (no flooding) old-growth forest and located in the same region, with similar climatic, edaphic, and floristic differences; (b) long-term records of satellite imagery and corresponding field data before and after disturbance are available; and (c) no subsequent disturbance has occurred.  Marra et al., 2018). In this region 93% of stems are between 10 and 40 cm in DBH (Higuchi et al., 2012) and the annual tree mortality is 8.7 trees ha −1 for trees ≥ 10 cm in DBH (Higuchi et al., 1997). 150

Landsat satellite data and procedures
The Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Schmidt et al., 2013;Masek et al., 2006;Masek et al., 2013;Masek et al., 2008) surface reflectance (SR) from Landsat 5 Thematic Mapper (TM) was used in this study to 155 characterize the type of disturbance and their subsequent pathways of forest regrowth over our study areas. LEDAPS was developed to ensure that spectral changes in Landsat are associated with regrowth dynamics (Masek et al., 2012;Schmidt et al., 2013) and to facilitate robust studies of land surface changes at different temporal and spatial scales in tropical forests (Kim et al., 2014;Valencia et al., 2016;Alonzo et al., 2016). LEDAPS SR Landsat 5 TM (L5 hereinafter) is generated by the United States Geological Survey (USGS) using the Second Simulation of a Satellite Signal in the Solar Spectrum (6S) that corrects 160 for the influences of, among others, water vapor, ozone, aerosol optical thickness, and digital elevation on spectral bands (USGS, 2017;Vermote et al., 1997). L5 bands are derived using per-pixel solar illumination angles and generated at 30-meter spatial resolution on a Universal Transverse Mercator (UTM) mapping grid (USGS, 2017). LEDAPS in Landsat 7 Enhanced Thematic Mapper Plus (ETM+) sensor (L7, hereinafter) were also used to corroborate our predictions (described below).
Though L5 and L7 use the same wavelength bands they are different sensors and differences in surface reflectance may exist, 165 especially in tropical forests due to high atmospheric effects (Claverie et al., 2015). Landsat 8 was not used since comparison between Landsat 8 and both L5 and L7 is not straightforward due to differences in the spectral bandwidth of these sensors. We used LEDAPS since a long time series of data is available with high spectral performance (Claverie et al., 2015) and it is suitable for ecological studies in the Amazon (van Doninck and Tuomisto, 2018;Valencia et al., 2016). L5 and L7 are available in Google Earth Engine (Gorelick et al., 2017), which we used to retrieve and analyze these data. 170 7/12/19877/12/ , 8/2/19897/12/ , 7/20/19907/12/ , 8/8/19917/12/ , 7/31/19947/12/ , 6/21/19977/12/ , 7/26/19987/12/ , 7/13/19997/12/ , 7/24/20037/12/ , 8/4/20077/12/ , 8/6/20087/12/ , 7/27/20107/12/ , and 8/31/20117/12/ . The dates of L7 images used were 8/7/20117/12/ , 6/22/20127/12/ , 6/12/20147/12/ , 8/2/20157/12/ and 8/7/2017 The spectral characteristics of old-growth forests and their changes after disturbances were investigated using 19 cells of 3×3 185 pixels (Figure 1 b,c,d). The average of each cell was used in our analysis. Spectral characteristics for old-growth forests for each site were determined from cells located in the same position of the disturbance but previous to disturbance, and/or from adjacent areas. Five cells of old-growth forests were located from 1 to 2 km away from the windthrow site. Though closer distances may also represent old-growth forests, we were conservative since Landsat is not sensitive to clusters of downed trees comprising fewer than 8 trees (Negrón-Juárez et al., 2011). For the clearcut and cut+burn sites the spectral characteristics 190 of their respective old-growth forests (control) were studied from 3 cells per site located 500 to 800 meters away from the edge of the disturbance to minimize edge effects that are relevant in the first 100 m (Lovejoy et al., 1986;Laurance et al., 2007;Mesquita et al., 1999). The spectral characteristics for the windthrow were acquired from two cells containing the highest level of SWIR1 in 1987 that is associated with the maximum disturbance (Negrón-Juárez et al., 2011;Magnabosco Marra et al., 2018;. For the clearcut site three cells were located 400-500 m from the edge and for the cut+burn three 195 other cells distant from 100-300 m with respect to the edge. For the clearcut site we also selected four areas: A1, A2, A3, and AT (AT = A1+A2+A3), shown in Figure 1c and for the cut+burn site, three areas: A1, A2, and AT (AT =A1+A2), shown in Figure   1d L5 data for the windthrow, clearcut, and cut+burn sites encompass a period of 28 years with 13 years of missing data due to 200 cloud-cover or lack of image. In order to assess the forest regrowth to spectral levels similar to old-growth forests (control), we applied a gap filling method (Gerber, 2018) of time series to obtain estimates for missing years using the R package "zoo" (Zeileis et al., 2018). The gap-filled datasets were analyzed using the smoothing spline technique (R, 2017). To determine whether L5 bands were sensitive to regrowth, we analyzed changes in the slope (β) of the bands across our chronosequence.
A t-test on the slope coefficient was used to test the null hypothesis that β is zero (H0:β=0) against the alternative hypothesis 205 (H1:β≠0) at a 5% significance level (α=0.05). If the computed test statistic (t-stat) was inside the critical values then the H0 was not rejected. The critical values (±t1-α/2, n-2, n is the number of data points) were obtained from statistical tables (Neter et al., 1988). Forests in the Manaus region affected by windthrows are dominated by tree species from genera Cecropia and Pourouma in about 3-5 years (Magnabosco Marra et al., 2018;Nelson and Amaral, 1994) and the clearcut and cut+burn sites were dominated by Cecropia and Vismia about 6 years after the disturbances (Mesquita et al., 1999;Mesquita et al., 2001). 210 The slopes of the time series were determined after these periods, i.e. 1991, 1987, and 1990 for windthrow, clearcut, and cut+burn sites, respectively.
A comparison of successional pathways of forest regrowth among studied disturbances was conducted that was feasible due to the similar conditions of climate, soils, and structure and composition of the old-growth forests. Time series of L5 bands regression ANOVA (analysis of variance) model. Calculations were conducted on the R 3.5.2 software platform (R, 2017) using the package gss (general smoothing splines) (Gu, 2018). We calculated the smooth spline (using the cubic fit algorithm) of observed data and the associated standard errors, from which we calculated Bayesian 95% confidence intervals. Predictions of the time after disturbance needed to reach old-growth forests values are based on these data using the function "ssanova" 220 (Fitting Smoothing Spline ANOVA Models) of the R package "gss" (General Smoothing Splines), Version: 2.1-9. The predictions were compared with published field observations (Section 2.1) where data were available and L7 images were used to assess the reliability of our predictions.

Forest regrowth simulation in ELM-FATES 225
Time series of L5 bands sensitive to disturbances and the pathways of forest regrowth were compared with modeling results from FATES (Fisher et al., 2015;Fisher et al., 2010;Holm et al. 2020). The underlying model structure and concepts in FATES are based on the Ecosystem Demography (ED) concept (Moorcroft et al., 2001), and is described in detail at https://github.com/NGEET/fates. A major development is the modularization of the model structure in FATES so that 230 boundary conditions and vegetation can be coupled with ESM land models. FATES is integrated into the E3SM Land Model (ELM) (Riley et al., 2018;Zhu et al., 2019) and within the Community Land Model (CLM) (Fisher et al., 2019;Lawrence et al. 2020) coupled to the Community Earth System Model (Hurrell et al., 2013). In this study we used ELM-FATES. ELM- based on parameter and demography sensitivity analysis at a site 40 km from the BDFFP (Holm et al., 2020), at the ZF2 research station (Magnabosco Marra et al., 2014). Holm et al. (2017) found that with the improved parameterization, ELM-FATES closely matched observed values of basal area, leaf area index (LAI), and mortality rates but underestimated stem density for a Central Amazon old-growth forest near the BDFFP.

245
Model simulations were driven by climate-forcing data derived from measurements collected between the years 2000 to 2008 at the K34 flux tower located at (2.6°S, 60.2°W) (de Araujo et al., 2002) about 40 km from the BDFPP, at the ZF2 research station. ELM-FATES (using the git commit "4a5d626" and the version corresponding to tag 'sci.1.0.0_api.1.0.0') was run and spun-up for 400 years until a stable biomass equilibrium was reached within the modeled forest. We then simulated a one-time logging treatment of a near-complete mortality of all trees (98% "clearcut") with a remaining of 2% consisting of only small 250 trees > 5cm DBH to aid in recruitment (Mesquita et al., 2015). Second, we simulated a one-time windthrow disturbance that killed 70% of trees ("windthrow") as was reported in a recent observational study on windthrows in the same region (Magnabosco Marra et al., 2018). All dead trees were 'removed' and therefore did not enter modeled soil pools. The fire module in ELM-FATES is currently under final development and testing and therefore burned simulations are not included in this study. The old-growth forest simulated by ELM-FATES, used as a pre-disturbance metric, was based on previously 255 validated tropical parameterization and sensitivity testing in the same region (Holm et al., 2020); and see supplemental material from Fisher et al. (2015) for description of plant functional type specific carbon allocation and allometry schemes and updates from the ED model framework. Simulations of disturbance and subsequent vegetation regrowth were initiated from this oldgrowth forest state. The model design used here only allows for simulating intact forests with natural disturbances (e.g., gap dynamics or windthrows) or harvested forests, but not both at the same time or in adjacent patches. Accounting for distance to 260 intact forests was excluded due to the current limited understanding of seed dispersal mechanisms (i.e., spatial variability, dispersal limitation, etc.) in tropical forests (Terborgh et al., 2019). We use a more general form of seed production, such that the individual cohorts in ELM-FATES use a targeted fraction of net primary production (NPP) during the carbon allocation process (after accounting for tissue turnover and storage demands), which adds to the site-level seed pool for recruitment of new cohorts. Field data was not used to simulate or calibrate the modeled forest regrowth post disturbance. 265 To account for uncertainty in the representation of plant physiology within tropical evergreen forests, we analyzed an ensemble of 20 simulations varying in targeted plant functional traits. We prescribed each ensemble with a single tropical evergreen plant functional type (PFT) that varied in wood density (0.44 to 1.06 g cm -3 ) and maximum rate of carboxylation (Vcmax; 42 to 55 umol m 2 s -1 ) ( Table 1), via random sampling. To evaluate changes in canopy coverage of the forest stand each PFT 270 additionally varied by an allometric coefficient (1.35 to 1.65) determining crown area to diameter ratio, and a leaf clumping index (0.59 to 1.0 out of 0-1 fraction) that determines how much leaf self-occlusion occurs and decreases light interception, and the direct and diffuse extinction coefficients in the canopy radiation calculations. The default values for these parameters are based on, or derived from references given in Table 1. Each ensemble member represents a single PFT across the spectrum of fast-growing 'pioneer' PFTs and slow-growing 'late successional' PFTs to provide a reasonable spread across the trait 275 uncertainty when assessing regrowth from disturbance. We characterized pioneer plants in our simulations as having low wood density (Baker et al., 2004) and high Vcmax based on the inverse relationship between these two plant traits, as well as a low crown area coefficient and low leaf clumping factor; i.e., monolayer planophile distribution (Lucas et al., 2002). These correlated relationships were applied in the ensemble-selected traits ( Figure 2). The opposite relationship was applied for slowgrowing, 'late successional' PFTs (e.g., high wood density, low Vcmax, high crown area coefficient, and high leaf clumping 280 factor).  In order to evaluate ELM-FATES performance during forest regrowth we compared NIR, the most sensitive band to regrowth (see results), with ELM-FATES outputs of aboveground biomass (AGB, Mg ha -1 ), total stem density of trees ≥10 cm DBH (stems ha -1 ), leaf area index (LAI, one-sided green leaf area per unit ground surface area, m 2 m -2 ), and total live crown area (m 2 m -2 ) since these variables directly influence the surface reflectance (Ganguly et al., 2012;Lu, 2005;Masek et al., 2006;Powell et al., 2010;Ruiz et al., 2005). We suggest that testing an array of modeled forest variables (e.g., biomass structure, density coverage of vegetation, and proportion of the tree crown that has live foliage) provides a robust comparison for comparison to NIR due to multiple forests characteristics contributing to and affecting NIR reflectance (Ollinger, 2011), and 305 reduces model unknowns and biases that can rise when using only one model variable. The usage of different stand structure and canopy processes can be helpful when evaluating ELM-FATES during different phases of forest regrowth. In addition, we averaged modeled outputs of crown area, stem density, and LAI since each of these variables influence the reflectance of forests, and defined this average as the modeled 'canopy-coverage'. Measurements of forest canopy cover have been used to analyze plant growth and survival, and it is an important ecological parameter related to many vegetation patterns (Ganey and 310 Block, 1994;Jennings et al., 1999;Paletto and Tosi, 2009). Modeled diameter growth rates (cm y -1 ) for trees with DBH ≥10 cm are also shown to provide information on the successional dynamics within ELM-FATES.

L5 bands and disturbances
All L5 bands showed an increase in surface reflectance immediately after windthrow, clearcut, and cut+burn sites except NIR which decreased (with higher decrease after burning) (Figure 3 a, b and c). This decrease in NIR was due to exposed woody material and dry leaves, typical after windthrow (Negrón-Juárez et al., 2010a;Negrón-Juárez et al., 2011) and clearcutting 320 (Chapter 7 in Adams and Gillespie, 2006) or the dark surface following burning (Pereira et al., 1997). representing the standard deviation of all pixels from respective cells. About one year after the disturbance the bands that experienced increases in surface reflectance showed a decrease in surface reflectance (the opposite for NIR) due to the 325 increases in vegetation cover. A similar response is expected for the clearcut that occurred in 1982 and therefore before the beginning of our available data (L5 imagery are available from 1984, Figure 3e). The similarity of spectral signatures for the control forests previous to the disturbances suggests comparable structure and floristic composition.

340
About six years after the disturbance, NIR reached a maximum value after and then decreased slowly with time showing a significant negative trend (Table 2). SWIR1 also showed a significant negative trend with time but only for the clearcut site (Table 2). In general, GREEN, BLUE, RED, SWIR1, and SWIR2 bands returned to pre-disturbance values (control) about six years after the disturbance (Figure 3d, e, f and Table 2). Therefore, we used NIR (which remained higher than predisturbance values throughout the time series, and is potentially sensitive to ecosystem properties of re-growing forest) to 345 investigate the regrowth dynamics in comparison to our control forests.
We used the relationships presented in Figures 4, 5, and 6 to determine the time that NIR from the disturbance sites became similar to control NIR. The average control NIR was 28±1%. For the windthrow site the NIR became similar to control levels after about 39 years (range 32 to 57 years). For the clearcut and cut+burn sites, this period was estimated to be 36 years (range 350 31 to 42 years) and 56 years (range 42 to 93 years), respectively. From Figure 4-6 it is evident that the type of disturbance has a clear effect on the pathways of NIR recovery. L7 data, in general, are within the 95% CI of predictions.   Figures 3d-f. The critical values (t0.975,8 and t0.975,12 )  During the first 12 years following the windthrow, the spline curve fitted to the NIR data decreased by ~0.13% y -1 after which the rate of decrease doubled (0.26% y -1 , Figure 4). For clearcutting, NIR decreased faster, i.e. ~0.4% y -1 . The decrease of NIR for the clearcut site appears to be independent of the distance from the edge of the disturbance since the changes in NIR of all 365 selected areas (A1, A2, A3 and AT) are similar ( Figure 5). For the cut+burn site, the rate of change of NIR to values similar to the control forests was the slowest among all disturbances considered (~0.15% y -1 ) ( Figure 6). The cut+burn site showed differences with respect to the border of the disturbance (areas A1, A2, and AT), which may be related to the spatial heterogeneity of burnings and forest responses.  Figure 1c) and prediction

of NIR to pre-disturbance values (in blue). The plots show the data (circles), the fit (solid line), and the 95% confidence interval (CI, dashed lines). Grey bar represents the control (old-growth forests) NIR of 28±1% and the black horizontal 380
dashed line is 28%.

FATES model and regrowth from forest disturbance
To address our goal of improving the connection between remote sensing, model benchmarking and the fidelity of future 390 predictions of forest regrowth processes, we examine the representation of such processes within ELM-FATES. The average of the ELM-FATES 20-member ensemble predicted a continuous, and almost linear, regrowth of biomass (Figure 7a) after clearcut and windthrows. The modeled recovering biomass returned to modeled old-growth forest values quicker for windthrows (37 years, range 21 to 83 years) compared to clearcuts (42 years, range 27 to 80 years). However, the annual rate of change of biomass regrowth over 50 years was faster in the clearcut simulation (2.5 Mg ha -1 yr -1 ) than the windthrow 395 simulations (2.0 Mg ha -1 yr -1 ), which was due to the clearcut site recovering from initial biomass of near-zero and proportionally greater contribution of fast-growing pioneer species.
The model simulation of stem density, LAI, and crown area are shown in Figures 7b-d, respectively. For stem density, ELM-FATES (black line) predicted an average of up to eight years before any new stems reached ≥10 cm DBH (a stand developing 400 period, Figure 7b) for the clearcut. The simulated stem density for old-growth forests (Figure 7b; green line) was ~200 stems ha -1 (≥10 cm at 1.3 m), ~ 408 trees lower than observations. The model ensembles with typical early successional traits predicted a forest with many fast-growing, small-diameter stems <10 cm, with maximum early successional stem densities reaching 1,560 and 1,414 stems ha -1 for clearcut and windthrows, respectively. Once the canopy closes and self-thinning dominates (average of 15 years after disturbance), there are declines in stem density as trees gain biomass, and canopy closure 405 forces some trees into the understory, where they die at faster rates due to shading. The modeled forests returned to old-growth stem density conditions 39 and 41 years after windthrows and clearcut, respectively (Table 3). Though the two disturbance types had very similar times to return to pre-disturbance conditions, they differ in the speed of recovery. ELM-FATES predicts faster diameter growth increments (1.3 cm yr -1 ) and canopy closure for a forest composed of all 'pioneer' type PFTs and slower (0.5 cm yr -1 ) diameter growth and more open canopies for 'late successional' forest type (Figure 8). Diameter growth is an 410 emergent model feature of dynamical plant competition for light and stand structure, and in coherence with observational studies of secondary forests growing through succession (Brown and Lugo, 1990;Winter and Lovelock, 1999;Chapin III et al., 2003).
The LAI of the modeled old-growth forest (4.0 m -2 m -2 ), prior to disturbances, was close to observed LAI (4.7 m -2 m -2 ) 415 measured near our study sites (Chambers et al., 2004). Due to disturbance, the initial modeled LAI (Figure 7c) and total crown area ( Figure 7d) decreased, as expected. During regrowth from disturbance both LAI and total crown area rapidly recovered, and LAI even surpassed pre-disturbance values. This pattern resembles the initial NIR spike due to fast growing PFTs. These two canopy coverage attributes reached maximum values after 3 to 6 years, depending on the disturbance and response of forest attributes (Table 3). To evaluate model results against remote sensing observations, we compared the initial period after 420 the disturbance of the spikes in NIR to the 'canopy-coverage' metric (combination of LAI, stem density, total crown area) over the same modeled period. ELM-FATES predicted that after a windthrow the forest took 5.7 years to reach maximum values of canopy-coverage, which was sooner than the clearcut simulation (7 years). While the modeled timespan for this initial period was similar to that inferred from NIR, there was disagreement between which disturbance recovery occurred fastest (windthrow in ELM-FATES vs. clearcut in NIR), similar to disagreement in recovery of AGB (Table 3).

that represented a fast-growing 'pioneer' forest stand and a slow-growing 'late successional' forest stand from a clearcut disturbance (black and grey) and a windthrow disturbance (orange). Variations in diameter increment are a result 440
of differences in the following traits: wood density, Vcmax, crown area, and a leaf clumping index.

Table 3. Summary of different times of regrowth (years) to old-growth forest status after two disturbance types;
windthrows and clearcuts from ELM-FATES model results and remote sensing. As well as the time (years) it takes forest attributes to reach maximum values during regrowth, and the corresponding value at this maximum peak. AGB 445 (Mg C ha -1 ), stem density (stems ha -1 ), LAI, and crown area (m 2 m -2 ) refer to simulation results, as compared against NIR remote sensing. The average of AGB and stem density is characterized as modeled 'forest structure'. The average of crown area, stem density, and LAI is characterized as modeled 'canopy-coverage' in this study, and additionally compared against NIR. ELM-FATES provided a prediction of the values in each forest variable when the stand reached its production limit and full canopy closure, at which point there was a shift to a declining trend and decreases in forest attributes that outpaced any gains 465 (Table 3). At maximum peak recovery and carrying capacity limit, the highest forest values occurred in the clearcut simulation, matching the higher NIR from clearcut a few years after the disturbance (Figure 4-6). Over the longer self-thinning period modeled LAI decreased and returned to modeled old-growth values 26 years after clearcut, and gradually over 53 years for windthrows (Table 3). LAI was the only variable that had a noticeable faster recovery in the clearcut simulations. After both disturbances the total crown area permanently remained high (0.99 m 2 m -2 ) and slightly higher than the crown area of the 470 simulated old-growth forests (0.98 m 2 m -2 ), suggesting that disturbances can generate a denser canopy, as discussed below.

Discussion
Our results show that Landsat reflectance observations were sensitive to the initial changes of vegetation following 475 windthrows, clearcut, and cut+burn, three common disturbances in the Amazon. Specifically, a decrease in NIR and an increase in SWIR1 were the predominant spectral changes immediately (within a few years) following disturbances. The increase in SWIR1 was different among the disturbances with the maximum increase observed in the cut+burn, followed by clearcut and then the windthrow site. The highest increase in SWIR1 in cut+burn sites may be related to the highest thermal emission of burned vegetation (Riebeek, 2014). Likewise, the relatively higher moisture content of woody material in the windthrow site 480 decreases the reflection of SWIR2. On the other hand, in our control (old-growth) forests, we observed typically high NIR reflectance due to the cellular structure of leaves (Chapter 7 in Adams and Gillespie, 2006), absorption of red radiation by chlorophyll (Tucker, 1979), and absorption of SWIR1 by the water content in leaves (Chapter 7 in Adams and Gillespie, 2006Gillespie, ). 2010bNegrón-Juárez et al., 2008), we found that NIR was more sensitive to the successional pathways of regrowth for all the disturbances considered. NIR has also been associated with succession (Lu and Batistella, 2005) and regrowth (Roberts et al., 1998) in natural and anthropogenic disturbed tropical forests (Laurance, 2002;Chazdon, 2014;Magnabosco Marra, 2016;Laurance et al., 2011). Previous studies have shown that changes in NIR are related to leaf structure and surface characteristic (Roberts et al., 1998;Xiao et al., 2014) with youngest leaves having higher NIR with respect to fully formed and older leaves 490 (Roberts et al., 1998). Maximum values of NIR were observed about 6 years after clear cut, which is the time pioneers form a closed canopy (Mesquita et al., 1999;Mesquita et al., 2015;Mesquita et al., 2001) and characterized by a relative uniform distribution of tree diameter and heights (Vieira et al., 2003). This maximum in NIR was higher in the clearcut site dominated by species from the genus Cecropia and Pourouma (Mesquita et al., 2015;Massoca et al., 2012) than the site affected by cut+burn dominated by Vismia species (Mesquita et al., 2015;Laurance et al., 2018). The higher NIR values in Cecropia and 495 Pourouma is due to their monolayer planophile distribution of large leaves that produced high reflectance compared to Vismia that have a rougher and denser canopies that traps more NIR (Lucas et al., 2002). The high values in NIR might be related to the low leaf wax (Chavana-Bryant et al., 2017) from new trees and/or scattering related to leaf and canopy water (Asner, 2008).
NIR decreases with the dynamics of succession due to increase in the canopy roughness (Hallik et al., 2019).

500
After the establishment of pioneers, the NIR decreases with time but with different rates depending on the type of disturbances.
In windthrown areas, tree mortality and subsequent recruitment may continue for several decades, promoting changes in functional composition and canopy architecture (Magnabosco Marra et al., 2018). Cecropia and Pourouma trees grow relatively quickly and after closing the canopy they limit light penetration due to their large leaves creating a dark, cooler, and wetter understory (Mesquita et al., 2001;Jakovac et al., 2014). As a result, light levels in the understory decline faster with 505 time and thus allow the recruitment and establishment of shade tolerant species. The cohort of Cecropia and Pourouma species has relatively short lifespan and, ~20 years after disturbance, secondary and old-growth forest species start to establish (Mesquita et al., 2015). With the self-thinning of Cecropia and Pourouma, the growing understory traps more light and consequently albedo decreases (Roberts et al., 2004). This pattern is consistent with the decline of NIR and observed changes in canopy architecture (Mesquita et al., 2015), photosynthesis, and LAI (Saldarriaga and Luxmoore, 1991). In contrast, the 510 architecture of Vismia species that dominate cut+burn areas allows higher light levels in the understory and subsequent recruitment of Vismia or other genera with similar light requirements. As a consequence, species turnover and structural changes are relatively slower than in clearcut areas (as found by Jakovac et al. (2014) in a study conducted a few kilometers from the BDFFP) and windthrows, consistent with changes in NIR. In the course of succession, Vismia tends to be replaced by Bellucia, which is a species with similar leaf and canopy structure as Vismia (Mesquita et al., 2015). This pattern favors the 515 penetration of light through the canopy  for several decades before a more shaded understory allows the germination and establishment of old-growth species .
For the windthrow we estimate that the NIR should become similar to pre-disturbance conditions in about 39 years. This value agrees with the 40 years of biomass regrowth found using ground-based data in the Central Amazon (Magnabosco Marra et 520 al., 2018). This result also corroborates previous studies that NIR operates in the best spectral region to distinguish vegetation biomass (Tucker, 1979(Tucker, , 1980 and photosynthesis (Badgley et al., 2017). For clearcutting and cut+burn, the regrowth time was about 36 and 56 years respectively, but no ground-based estimates were available for comparison. Still, NIR showed that the pathways of regrowth from clearcut and cut+burn are divergent with time ( Figure 5 and Figure 6), which is consistent with observational studies (Mesquita et al., 2015). 525 In general, we found that NIR may be used as a proxy in modeling studies aimed at addressing forest regrowth after disturbances. Though NIR is useful to distinguish successional stages up to decades after the disturbance, it may not represent the whole successional processes. As soon as the forest canopy becomes structurally similar to that of the mature forest, NIR will no longer be sensitive to changes in vegetation attributes (Lucas et al., 2002). Though L5 NIR may be complemented with 530 current Landsat measurements (L5 NIR has comparable performance to the Landsat 8 NIR Operational Land Imager algorithm (OLI) (Vermote et al., 2016)), it is important to emphasize that our estimates of recovered reflectance and biomass in disturbed areas do not capture full recovery of diversity in floristic attributes and species composition that can take centuries (Rozendaal et al., 2019). The predominance of Cecropia, after clearcut, and Vismia, after cut+burn, have also been found in the Western (Gorchov et al., 1993;Saldarriaga et al., 1986) and the Southern (Rocha et al., 2016) Amazon suggesting that our findings are 535 applicable to other regions. However, an Amazon-wide study is beyond the scope of our work, which is to explore the sensitivity of Landsat to different disturbance types.
Our analysis demonstrates that this version of ELM-FATES has the capacity to reproduce initial response to disturbance and regrowth patterns after the clearcut and windthrow that occurred over similar time ranges compared to NIR. The strongest 540 agreement occurred when ELM-FATES predicted higher peaks of post disturbance stem density and LAI in clearcuts than in windthrows, which can be used for future benchmarking, consistent with the higher peak of NIR from clearcuts ( Figure 5 vs. Figure 4). This effect may be due to ELM-FATES having more homogeneous canopies after clearcuts as well as more open disturbed area for fast growing plants, which is also an observed trend. In addition, the average regrowth times to predisturbance values were close between ELM-FATES and NIR (Table 3), showing that pathways of forest regrowth in ELM-545 FATES are comparable to observed patterns in tropical forests. ELM-FATES predicted a continuous, and almost linear, regrowth of biomass for the first 50 years of simulation after both clearcut and windthrows, (Figure 7a) consistent with NIR results, and observational studies (Mesquita et al., 2015;Saldarriaga et al., 1988;Jakovac et al., 2014;Magnabosco Marra et al., 2018). In addition the changes in biomass rates predicted by ELM-FATES were similar to biomass observations recorded after clearcut (2.3 vs. 2.6 Mg ha -1 yr -1 ) (Mazzei et al., 2010), as well as a faster rate of AGB accumulation after clearcut 550 compared to windthrow, similar to a study reporting higher regrowth rates in more highly-disturbed sites (Magnabosco Marra et al., 2018).
Landsat showed a faster recovery of NIR to pre-disturbance conditions in clearcuts compared to windthrows. Faster growth is characteristic of anthropogenically driven secondary forests that reflect rapid colonization and monodominance of adapted 555 species and genera to changed environmental conditions; e.g., high growth rates, low self-competition, high leaf area index, low herbivory rates, etc. (Poorter et al., 2016;Mesquita et al., 2015;Rozendaal and Chazdon, 2015). Alternatively, ELM-FATES predicted a faster recovery of structural AGB and canopy-coverage to pre-disturbance conditions for windthrows (70% tree mortality) compared to clearcuts (98% tree mortality). A major contributing factor to this pattern resulted from larger modeled diameter increments after windthrows (0.92 cm yr -1 ) compared to clearcuts (0.82 cm yr -1 ) in the first 20 years after 560 the disturbance, setting the trajectory for faster regrowth to pre-disturbance after windthrows. Only LAI had a faster recovery to pre-disturbance values after clearcuts, which is expected due to the newly developed forest having simplified forest structure and the canopy being more homogeneous after a clearcut (Rosenvald and Lohmus, 2008). ELM-FATES predicted the timing of peak canopy-coverage was marginally sooner after windthrows compared to clearcuts, opposite to the NIR pattern. This discrepancy may be related to more biomass loss and open canopy coverage, followed by a lack of rapid colonization in the 565 modeled clearcut. Due to the higher stand disturbance that naturally occurs from clearcuts, and the diverse complexities in tropical forest composition, we emphasize that the dynamics of different competing PFTs in ELM-FATES requires further investigation. Emerging modeling studies that include plant trait trade-offs, for example, in leaf and stem economic spectrum or fast vs. slow growth strategies may help to better capture the drivers of forest productivity and demography, enabling improved modeled responses to global change scenarios (Fauset et al., 2019;Sakschewski et al., 2015). Here we test the basic 570 representation of biomass demographics, prior to the more challenging aspects of representing interacting functional diversity in recovering systems (Fisher et al., 2015;Powell et al., 2018).
ELM-FATES predicts the total stem density for a closed canopy forest to be very low compared to observations (200 simulated vs 600 stems ha -1 ), but modeled total AGB was close to that reported for the same region (110 simulated vs. ~150 MgC ha -1 575 (Chambers et al., 2013)). This discrepancy is due to ELM-FATES predicting a disproportionately high number of large trees ( Figure 9; with 8% of stems >60 cm and 4.5% of stems >100 cm), resulting in a crowded canopy, which out-compete smaller understory trees. Higuchi et al. (2012) reports 93% of trees ≥ 10 cm DBH in the study site to be below 40 cm DBH, while ELM-FATES predicted noticeably less trees below 40 cm DBH (Figure 9). Low stem density could be attributed to multiple model assumptions, such as high density-dependent mortality and self-thinning due to the marginal carbon economics of 580 understory trees, low branch-fall turnover, a need for greater limitation of maximum crown area than currently modeled, and increases in mortality rates with tree size (Johnson et al., 2018). Our findings here will guide future ELM-FATES and ecosystem modeling development efforts towards improving the representation of forests comprised of dense canopies, and how they shift during regrowth.
Land surface models do not typically simulate spectral leaf reflectance, but there is potential to include such output within radiative transfer schemes as is currently done in CLM. That addition would greatly assist our ability to compare with Earth observations datasets. In lieu of this development, we show that with successional aging modeled forest structure returns to pre-disturbed values (through canopy closure) with similar recovery time as inferred from NIR, occurring with the process of canopy closure, all which can be compared against remote sensing vegetation indices (see Supplementary Figure 1) and 590 metrics. Which vegetation index (e.g., Normalized Difference Vegetation Index (Rouse et al., 1973), Enhanced Vegetation Index (Huete et al., 2002), etc.) or metric properly represent the successional pathways following disturbances remains an important area of study. 595 Figure 9: Total stem density (stems ha -1 ) separated into six diameter (cm) size classes from Central Amazon field data located at the near-by ZF2 site (green bars) averaged from 1996-2011, and predicted by ELM-FATES (gray bars). The percentages represent the proportion of stems in each size class relative to the total stem density.

Conclusions
We tested the sensitivity of Landsat surface reflectance to windthrow, clearcut, and cut+burn forest in Central Amazon. NIR was more responsive to the successional pathways of forest regrowth years after the disturbance. NIR showed that pathways 605 of forest regrowth were different among the disturbances, with cut+burn being the most different in terms of spatial