Vulnerability of Amazon forests to storm-driven tree mortality

Tree mortality is a key driver of forest community composition and carbon dynamics. Strong winds associated with severe convective storms are dominant natural drivers of tree mortality in the Amazon. Why forests vary with respect to their vulnerability to wind events and how the predicted increase in storm events might affect forest ecosystems within the Amazon are not well understood. We found that windthrows are common in the Amazon region extending from northwest (Peru, Colombia, Venezuela, and west Brazil) to central Brazil, with the highest occurrence of windthrows in the northwest Amazon. More frequent winds, produced by more frequent severe convective systems, in combination with well-known processes that limit the anchoring of trees in the soil, help to explain the higher vulnerability of the northwest Amazon forests to winds. Projected increases in the frequency and intensity of convective storms in the Amazon have the potential to increase wind-related tree mortality. A forest demographic model calibrated for the northwestern and the central Amazon showed that northwestern forests are more resilient to increased wind-related tree mortality than forests in the central Amazon. Our study emphasizes the importance of including wind-related tree mortality in model simulations for reliable predictions of the future of tropical forests and their effects on the Earth’ system.


Introduction
The productivity (22 PgC year −1 , or 35% of global terrestrial productivity) (Pan et al 2013, Fernandez-Martinez et al 2014 and aboveground biomass (AGB, 98 PgC, or 25% of terrestrial biomass) (Pan et al 2013, Malhi et al 2011 of the Amazon (representing 53% of global tropical forest area; supplementary figure S1 available at stacks.iop.org/ERL/13/054021/mmedia) are well quantified, but the mechanisms affecting their spatial patterns remain uncertain (Stephenson et al 2011, Malhi et al 2015, Negrón-Juárez et al 2015. Tree mortality is an important component of forest dynamics (Kellner and Asner 2009), and affects ecosystem processes (Frelich 2002, dos Santos et al 2016 as well as spatial patterns of productivity and biomass (Malhi et al 2015, Stephenson et al 2014, Negrón-Juárez et al 2015, Johnson et al 2016. For instance, the less-deforested and rainy northwestern Amazon (NWA) (no dry season, monthly rainfall > 100 mm (Sombroek 2001 figure 1(a)) between 2013 and 2015 represented for dots of different color and sizes described in the legend. The thin black contour lines are included to better visualize the spatial patterns of windthrows. The thick blue contour line represents the area without dry season. Inset shows the association between windthrow size and rainfall from TRMM 3B42 V7 3 Hour 0.25 • × 0.25 • data (y = 33.32x-2.19, r 2 = 0.63, p < 0.001). The dashed lines represent the 95% confidence intervals. (c) Mean annual number of rainfall events ≥10 mm hr −1 in the Amazon (base period 1998-2016) using TRMM 3B42 data. biomass (Baker et al 2004, Malhi et al 2006 than the Central Amazon (CA), which has three months of dry season (consecutive months with rainfall ≤ 100 mm, (Sombroek 2001)). The higher tree mortality in NWA may help explain these patterns (Stephenson and van Mantgem 2005, Galbraith et al 2013, Malhi et al 2015. Yet, the variation in mortality rates across the Amazon is not fully explained by factors intrinsic to forests (i.e. factors influenced by or resulting from the ecosystem (Wenger 1984)), such as competition, growth, defense strategies, and soil texture and nutrients (Malhi et al 2015, Chao et al 2008, Stephenson et al 2011, Quesada et al 2012, Laurance et al 1999, Higgins et al 2011. Processes producing tree mortality are not fully represented in CMIP5 ESMs (Coupled Model Intercomparison Project Phase 5 Earth System Models) (Taylor et al 2012), and may constitute an important and unassessed source of uncertainty.
One major contributor of tree mortality in the Amazon are downbursts (strong descending winds) associated with severe convective systems as squall lines (Garstang et al 1998, Negrón-Juárez et al 2017, Cohen et al 1989. Downbursts create gaps of uprooted or broken trees, windthrows (Mitchell 2013). Previous studies focused over Brazil used Landsat imagery to identify windthrows ranging in size from a single pixel (30 m × 30 m) to hundreds of hectares (Nelson et al 1994, Chambers et al 2013, Espirito-Santo et al 2010, Negrón-Juárez et al 2011. Our objectives in this study are to use: (i) Landsat images between 2013 and 2015 to identify the spatial occurrence of windthrows (<1 yr old) across the whole Amazon; (ii) surveys of tree mortality over windthrown areas to determine tree mortality rates and floristic composition; and (iii) chronosequences of Landsat images to compare the occurrence of windthrows over NWA and CA. Our final objective is to use these data to inform a modeling experiment of forest response to increased wind-related tree mortality rates likely to accompany projected increases in extreme rainfall events in the Amazon (sections 2.6.2 and 14.8.5 in IPCC 2013).

Research sites and climate
The study areas are the regions of Iquitos (located in NWA) and Manaus (located in CA) (figure 1(a)). In Iquitos MAR (mean annual rainfall) is ∼3000 mm, MAT (mean annual temperature) is 25.9 • C, and June through August are the months with the lowest precipitation (∼183±10 mm month −1 , figure S2(a)). In Manaus, the MAR and MAT are 2300 mm and 26.2 • C, respectively, and the region experiences three consecutive months of dry season (June through August, ∼65 ± 10 mm month −1 , figure S2(b)). The research sites in Iquitos and Manaus (described in sections 2.2 and 2.5) are lowland old-growth forests with trees ∼30-40 m tall where windthrows have occurred.

Landsat imagery and windthrows
To determine the spatial distribution of recent windthrows ≥ 25 ha across the Amazon (figure 1(a)) 229 Landsat 8 images (L8, each covering ∼3.4 × 10 4 km 2 , 30 m × 30 m of pixel resolution) with low cloud cover (<20%) between 2013 and 2015 were used. Windthrows were identified by their distinctive fan-shape diverging from a central area with radiant corridors separated by standing forest (Nelson et al 1994) and spectral characteristics (Negrón-Juárez et al 2010). To assess the severity of windthrows, spectral mixture analysis (SMA) (Adams et al 1995, Shimabukuro andSmith 1991) was applied to bands 2-7 in L8 scenes. SMA quantifies the per pixel fraction of endmembers which sums to match the full pixel spectrum of the image (Adams et al 1995). Image-derived endmembers of photosynthetic (green) vegetation (GV), non-photosynthetic vegetation (NPV), and shade were used. The fractions of GV and NPV were then normalized without shade (Adams and Gillespie 2006) as GV/(GV+NPV) and NPV/(GV+NPV). The normalized NPV images were used to calculate the total area of the selected windthrows (supplementary text S1). With this data we produced a windthrown area weighted 2D Kernel Density Estimate (KDE) with the SAGA geospatial analysis library (www.saga-gis.org) using a quartic kernel with a 100 km radius, resampled to a 1 km resolution. The KDE is plotted as a contour plot showing the probability that a windthrow greater than one hectare will occur within each 1 km 2 pixel ( figure 1(b)).
Chronosequences of Landsat imagery were used to calculate the return frequency of windthrows in Iquitos and compared with frequencies from our previous study over Manaus (Chambers et al 2013). Landsat 5 Thematic Mapper images (L5) were obtained from the United States Geological Survey (http://glovis.usgs.gov), and atmospherically corrected and converted to reflectance using the Atmospheric CORrection Now (ACORN) software (ImSpec LLC, Boulder, CO). The Carlotto (Carlotto 1999) technique, which corrects for haze and smoke contamination, was applied over the scenes as needed. All anthropogenic areas, water bodies, clouds and theirs shadows were masked out. We used ten L5 images for five intervals (pairs) for Iquitos (Landsat scene P006R063, 1988(Landsat scene P006R063, -1985(Landsat scene P006R063, , 1991(Landsat scene P006R063, -1989(Landsat scene P006R063, , 1999(Landsat scene P006R063, -1995(Landsat scene P006R063, , 2002(Landsat scene P006R063, -2001(Landsat scene P006R063, , 2009(Landsat scene P006R063, -2005 each on approximately the same day of the year. Paired images were intercalibrated (old image, e.g. 1985, with respect to new image, e.g. 1988) band by band by regressing the encoded radiances using temporally invariant targets (Negrón-Juárez et al 2011), and the normalized (without shade) NPV was calculated. Changes in NPV for each period (ΔNPV, a quantitative measure of the changes in dead biomass associated with disturbance (Negrón-Juárez et al 2011)) were calculated (normalized NPV newimage -normalized NPV oldimage ).
The ΔNPV was also used to determine the location of plots (∼3 ha and divided in 100 subplots of 0.03 ha) for conducting forest inventories over windthrown forest (covering a gradient of disturbance severity) and adjacent undisturbed forests. For each plot, we measured tree diameters for all trees ≥ 10 cm in DBH (diameter at the breast height = 1.3 m above ground or above buttress) and measured trees were identified. We studied two windthrows in the Iquitos region (Allpahuayo, in the Allpahuayo-Mishana National Reserve and Nauta, in the Nauta Research Station, managed by National University of the Peruvian Amazon-UNAP) (Rifai et al 2016)

Tropical Rainfall Measuring Mission data
Data from the 3 Hour Multi-satellite Precipitation Analysis (3B42) of Tropical Rainfall Measuring Mission (TRMM) (Huffman et al 2007) version 7 for the period 1998 to 2016 was used to determine areas with heavy rainfall as well as the association between heavy rainfall and windthrows (text S2). We performed these estimates on two data cases: (1) the remotely sensed data only (windthrow from Landsat), and (2) the remotely sensed data plus inventory plots data (Landsat+plot). These data was multiplied by the total analyzed forested area across the five Iquitos Landsat scene-pairs (31 900 km 2 ) described in section 2.2. Tree mortality data from Landsat included windthrows with more than eight dead trees based on the average mortality rate as function of gap size (figure S3) as well as data from eight-tree events from single-pixel disturbances (Negrón-Juárez et al 2011). Plot data include mortality of single trees (Chao et al 2008). The OLS method used log-binned data with the starting bin center at event of eight trees in size, while the MLE method used the data directly. An extensive description of the power law fit method used is described in Di Vittorio et al (2014). To directly compare this average mortality event-size distribution with those in Manaus (Chambers et al 2013) (that uses the same method) we aggregated the distribution into the same set of discrete event-size bins which were defined in logarithm with base 10 for Manaus, plus one extra bin to cover the larger clusters detected in Iquitos. Probability distribution-function of event size and dead trees is presented. Further details are presented in text S3.

Comparison of Iquitos and Manaus windrelated tree mortality
The Nauta windthrow (November 2, 2009) in Iquitos was compared with the previously studied ZF2 windthrow in Manaus (Negrón-Juárez et al 2010). For the Nauta windthrow, L5 scenes (Path 006/Raw063) from June 27, 2008, and December 07, 2009 were used to calculated ΔNPV. The ΔNPV image was used to randomly select thirty disturbed pixels covering the whole gradient of disturbance. Each pixel was located using a handheld GPS (60CSx GPS, Garmin Ltd.) and a plot of 20 m × 20 m was installed. The number of live and dead trees with DBH ≥10 cm was recorded in the surveys. The survey was done 10 months after the occurrence of the windthrow. Dead trees were identified by the direction of fall similar to the predominant direction of the windthrows, presence of leaves, preserved trunk bark, and fresh fiber characterizing snapping. An ANCOVA analysis (using SAS 9.2, SAS Institute Inc., North Carolina, USA) was used to compare the association between tree mortality and ΔNPV in Iquitos and Manaus. The associated downburst velocities for the Nauta windthrow were estimated as 19-31 m s −1 (68-112 km h −1 ) (figure S4).

Simulation of forest and windthrows dynamics
We used the ZELIG-TROP model (Holm et al 2014) to simulate forest and tree mortality interactions due to increased windthrow events. ZELIG-TROP is an individual-based demographic gap model with trees of varying size and structure, as well as a dynamic vegetation model in that regeneration, growth, and mortality rates all dynamically vary as a function of climate, stand density, canopy openness, and competition for resources. The data used to parameterize the species traits and demographic attributes (such as maximum DBH, age, growth-rate scaling coefficient, stress tolerances, and recruitment rates) was obtained from our inventory plots, RAINFOR data ( (see table  S2).
We simulated 50 spatially explicit plots, each at the patch-dynamic scale of 400 m 2 (totaling two hectares). The simulations started from bare ground (as typical for successional-based gap models), thus simulating dynamic and stochastic individual treerecruitment, growth, and mortality within competing environmental conditions for Manaus and Iquitos. Model simulations were run for 400 yr until the forest reached a mature steady state (compared to field data in table 1), and all results are averaged over the final 50 yr of simulation. The tree mortality rate (% stems yr −1 ) was simulated to be 2.1% for Iquitos (i.e. control), which is very close to observed values (i.e. ∼2.3%) (Chao et al 2008). The control simulations agreed with observed values of forest attributes (table 1) and predictors of tree mortality (text S4 and figure S5) from a previous tree pulling study (Ribeiro et al 2016).
To study how increases in windthrows may affect community composition, the current background mortality (M) in Iquitos (∼2%) and Manaus (∼1%) was near doubled (2M) (table 1), assuming that all increase in mortality is associated with windthrows. A justifiable value of increased mortality associated with windthrows is not currently possible since model predictions of future climate do not include windthrows as a form of disturbance. To increase mortality to 2M the model randomly selected trees ≥10 cm DBH to die and are removed from the forest on an annual basis. Control simulations were denoted as Iquitos M and Manaus M, and doubled-mortality simulation as Iquitos 2M and Manaus 2M.
Sixty eight species belonging to 46 genera were used in the Iquitos simulations (table S2). For Manaus simulation, 90 species belonging of 51 genera were used (Holm et al 2014). Eighteen genera (figure S6) were found common in both Amazon regions. These 18 common genera accounted for 75% and 47% of the total basal area (BA) in Iquitos and Manaus, respectively. These overlapping genera were used to compare shifts in community composition between the M and 2M simulations in these regions (figure S6 (a) and (b)). Using the 18 common genera allowed for a stricter and standardized basis for comparing compositional shifts and alleviate biases that would manifest when including genera that are only present in one region and not the other.
A non-metric multidimensional scaling (NMDS) ordination analysis (Shepard 1962, Kruskal 1964a, 1964b, Faith et al 1987 on genera abundance classified by BA (m 2 ha −1 ) was used to analyze the plot data and model results. We used this technique to reduce the multidimensionality of our studied communities Table 1. Model results for Iquitos. Control outputs over Iquitos (± standard deviation) compared against field data and disturbance model outputs. The tree mortality rate for Iquitos is the exponential mortality (2.3%) (Chao et al 2008) * , which can be considered the true annual mortality since both differ only by about 0.02% (Sheil et al 1995). Since ZELIG-TROP is an individual-based model that prognoses mortality mechanistically, mortality is reported as individual level annual mortality. ANPP = aboveground net primary production. into a two dimensional-ordination space and assess possible differences in genera composition. The NMDS from 124 inventory plots distributed in Iquitos and Manaus (supplementary table S1) was compared against simulations. For M and 2M simulations, the NMDS analysis was based on 20 plots per region. Stress scores (that represents the difference between the reduced and the complete multidimensional space) < 0.3 indicates a good representation of data in the reduced dimensions.

Results
Between 2013 and 2015, the areas of the Amazon basin without dry season (figure 1(a)) were most affected by windthrows ( figure 1(b)). These areas, which include northwestern Peru, southern Colombia, southern Venezuela, and the central and western regions of the Brazilian Amazon, also experience more frequent and intense hourly rainfall rate events (figure 1(c)). A positive association was found between hourly rainfall rates and windthrows (inset in figure 1(b)). The Iquitos power law exponents range from −2.26 to −2.55, with a mean value of −2.44 (figure S7) which are outside the range of the Manaus exponent (−2.71 and −2.97) (Chambers et al 2013, Di Vittorio et al 2014. These power law exponents do not change substantially between methods (OLS and MLE) or data used (Landsat and Landsat+ plot data) (figure S7) and show that windthrows are more frequent in Iquitos. In Iquitos 37% of annual tree mortality associated with windthrows occurred in gaps formed from the death of eight or more trees (bins ≥ 8 trees) (table 2). This proportion is more than double the rate previously estimated for CA (17%) (Chambers et al 2013).
The downburst speeds that produced the Nauta windthrow (figure 2(a)) were less intense (estimated as 68-112 km h −1 , figure S4) than over the ZF2 forest (estimated as 93-147 km h −1 ) (Negrón-Juárez et al 2010) that was expected since slower moving system (12 m s −1 in Nauta (figure S4) vs 19 m s −1 in ZF2 (Negrón-Juárez et al 2010)) produce lower downburst velocities (Garstang et al 1998). A strong linear relationship (p < 0.001 and r 2 = 0.8; figure 2(b)) was found between observed tree mortality and ΔNPV in Nauta and ZF2. The ANCOVA analysis showed that the slopes (112 ± 6.0 (SE) and 99 ± 5.6 for Nauta and ZF2) were not statistically different (F = 2.38, p = 0.12), even for the unforced through origin case (F = 0.16, p = 0.69). Despite the fact that downburst speeds were lower during the Nauta windthrow than the ZF2 windthrow, the tree mortality per area associated with these storms was similar suggesting that trees in Iquitos had a higher sensitivity to winds.
The control simulations predicted the community composition in both regions accurately (figure S6(a) and figure S6(b) compared to figure 3(a) and (b)), with an exception of the lower dominance of Eschweilera (Lecythidaceae) in Manaus. The control simulations also reproduced the differences in genera composition observed between Manaus and Iquitos (figure 3(c)), as indicated by the distances between genera similarity derived from the NMDS ordination (i.e. separation of red and blue points in figure 3(d)). Iquitos 2M resulted in a significant, yet smaller decrease in biomass compared to the biomass decrease in Manaus 2M (29.6%; two sample t test, t (99 1.97) = 186.20, p < 0.001 vs. 41.9%; two sample t test, t (99 1.97) = 108.98, p < 0.001, respectively).
Iquitos 2M did not lead to a significant shift in genera composition with respect to Iquitos M (figure 3(d)) (Wilcoxon rank sum, Z = 0.33, p = 0.74, and also seen in a cluster analysis in figure S8(a)). In contrast, Manaus 2M had significant shifts in the genera composition with respect Manaus M (Wilcoxon rank sum, Z = 2.28, p = 0.02, figure S8(b)). Furthermore, Manaus 2M had community similarities to Iquitos M (figure 3(d) and figure S8(c)) but differences exist (Wilcoxon rank sum, Z = 1.03, p = 0.3). Text 5 summarize model changes in AGB, NPP, BA and diversity index associated with doubling the mortality rate. In both Iquitos 2M and Manaus 2M the model predicted that emergent trees had the largest decrease in basal area (after weighting by abundance) ( figure S6(c)). This decrease was stronger in Manaus 2M (19%) compared to Iquitos 2M (14%).

Discussion
We found that windthrows are spatially more common in NWA (north of the Ecuadorian Amazon, northeast of Peru, and northwest of Brazil) ( figure 1(b)) consistent with and expanding on previous studies focused on Brazil (Nelson et al 1994, Espirito-Santo et al 2014, Espirito-Santo et al 2010. Our inventory plots vary in important attributes (e.g. biomass, size-distribution of trees, tree height, stem density and community mean wood density (Quesada et al 2012, Feldpausch et al 2011, Mitchard et al 2014, Baker et al 2004) that characterize the typical gradient across the Amazon driven by climatic, edaphic and ecological aspects (Quesada et al 2012). Thus, our inventory plots in Iquitos and Manaus may be considered representatives of NWA and CA respectively.
The power law analysis (figure S7) suggests that windthrows occur more frequently in NWA than in CA (table 1 in Chambers et al 2013) across all size classes. NWA and CA have distinctive community composition, which partially can be related to treemortality patterns, as discussed in previous studies (Chao et al 2008, ter Steege et al 2006, ter Steege et al 2013, Magnabosco Marra 2016. We found that NWA is more vulnerable to wind-related tree mortality that may be related to a number of factors including tree stem density, tree size, species composition, wood density, soil type, root architecture, and topographic exposure (Boose et al 2004, Vandermeer et al 2000, Chao et al 2008, de Toledo et al 2011, Magnabosco Marra 2016, Ribeiro et al 2016. The similarities in tree density (supplementary table S1) and tree size-distribution (supplementary figure S9) in NWA and CA suggest that these structural attributes are unlikely to be major contributors of the observed differences in tree mortality, in agreement with previous studies (Chao et al 2008).
Winds produce tree mortality by uprooting (more likely associated with soil characteristics) or breaking (more likely associated with mechanical failure) (Rankin-de Merona et al 1990, Korning andBalslev 1994). The soils in the NWA are dominantly Acrisols with low infiltration rates that limit root growth due to high soil compaction, while soils in CA are primarily Ferrasols which have physical properties that allow development of deeper root systems (Hengl et al 2014). These soil characteristics correspond with gradients of root depth across the region (Ichii et al 2007, Nogueira Lima et al 2012. Furthermore, in NWA light limitation, abundant rainfall and soil conditions (Nemani et al 2003, Malhi and Davidson 2009, Malhi et al 2004, ter Steege et al 2006, Baker et al 2004, Quesada et al 2012 promote rapid vertical growth (Stephenson et al 2011) and shallow root systems (Quesada et al 2012). Frequent extreme rainfall events in NWA (figure 1(c)) might also result in saturated soils (Foster 1988) and strong wind loads affecting the mechanical stability of trees by swaying (Gardiner et al 2016). These characteristics (individually or combined) limit the anchoring of trees into soil and may make NWA trees more sensitive to winds ( figure 2(b)). Tree-pulling experiments in Iquitos, similar to those we have conducted in Manaus (Ribeiro et al 2016), are needed to test these assumptions.
Our modeling results suggest that community composition in Manaus is more vulnerable to an increase in windthrows with respect to Iquitos, consistent with previous studies showing that forests with low mortality rates are more sensitive to increases in mortality rates (Schietti et al 2016, de Toledo et al 2013, Johnson et al 2016. Manaus 2M predicted community similarities to Iquitos M, but differences exist. These differences may be related to either individual or combined effects, such as (i) differences in soil characteristics and nutrients (Quesada et al 2012); (ii) local effects of disturbances; or (iii) model limitations, as the model cannot incorporate all complex and diverse plant physiological processes (Wright 2002), all of which are challenging data to collect in the tropics.
Our modeling results suggest that the community composition in Iquitos did not significantly change with increases in tree mortality associated with windthrows. The higher frequency of windthrows in Iquitos (with respect to Manaus) may have resulted in a forest that is more resilient (as defined by Holling 1973) to increases in tree mortality rates from windthrows. This resiliency was represented by the overlap of Iquitos M vs. Iquitos 2M points in figure 3(d), as well as the smaller decrease in biomass in Iquitos 2M. Iquitos might only see a significant change in community composition with a larger increase in mortality or different drivers of mortality.
Our results shows that Landsat and Landsat+plot data do not produce changes in the power law exponent suggesting that an important fraction of plot level tree mortality is produced by winds. However, we emphasize the need for additional field observations of different mortality agents for proper region-specific attributions of cause and effect, as well as vulnerability, resistance, and responses of different species to those agents. Such observations would improve understanding of forest dynamics beyond the importance of general tree mortality in explaining patterns of biomass and productivity in the Amazon. Our results directly address the need to incorporate the effects of wind-related tree mortality on ecosystem processes in ESMs to reduce uncertainties of carbon and climate projections.

Conclusions
We found that windthrows are common in the region extending from the northwest Amazon (northeastern Peru, southern Colombia and Venezuela and northwest Brazil) to central Brazil, with the highest occurrence of windthrows found in northwest Amazon. The more frequent extreme winds associated with more frequent severe convective systems in NWA may explain the higher tree mortality observed in this region. The higher frequency of windthrows in NWA may have resulted in a forest that is more adapted to these disturbances. Our model results suggest that increases in the occurrence of windthrows may produce a shift in composition in CA but not in NWA.