Large-scale variations in the dynamics of Amazon forest canopy gaps from airborne lidar data and opportunities for tree mortality estimates
Study area
The study area was the Amazon basin (Fig. 6)42. The basin was sub-divided in four regions with markedly different forest dynamics, geography and substrate origin, adapted from the classification of Feldpausch et al.33: West (parts of Brazil, Colombia, Ecuador and Peru), Southeast (Bolivia and Brazil), Central-East (Brazil) and North (Brazil, Guyana, French Guiana and Venezuela). The natural vegetation mainly corresponds to broadleaf moist forests and tropical seasonal forests, with both terra firme and seasonally flooded forests. Across the Amazon, there is a wide range of average monthly rainfall (100–300 mm) and dry season length (DSL) (0–8 months)43.
Figure 6
The Amazon in South America with colored regions, defined in Feldpausch et al.33, indicating faster (West and Southeast) and slower forest dynamics (Central-East and North). Small black lines represent single-date airborne lidar data acquisitions from the EBA project (n = 610 flight lines). Red triangles illustrate multi-temporal lidar data acquisition over five sites (BON, DUC, FN1, TAL and TAP). Circles indicate the location of field inventory plots (n = 181). R v4.0.2 was used to plot this figure32.
Full size image
The five sites selected for the multi-temporal assessment of the static and dynamic gaps relationship (red triangles in Fig. 6) were: Adolpho Ducke forest (DUC), Tapajós National Forest (TAP), Feliz Natal (FN1), Bonal (BON) and Talismã (TAL). These areas were chosen to represent distinct forest types, vegetation structure and biomass stocks. The predominant vegetation types consisted of dense rain forests (DUC and TAP), seasonal forests (FN1), and open rain forests (TAL and BON). DUC and FN1 are mostly undisturbed forests, while TAP underwent fire and/or selective logging in the past. TAL and BON were affected by a known fire occurrence in 2010. The sites cover a gradient of aboveground biomass (AGB) that increase, in average, from TAL (185 Mg ha−1), FN1 (235 Mg ha−1), BON (251 Mg ha−1), and DUC (327 Mg ha−1) to TAP (364 Mg ha−1)44.
Data acquisition and pre-processing
Airborne lidar data
Multi-temporal lidar data were obtained by an airplane at each of the five sites (red triangles in Fig. 6), as part of the Sustainable Landscapes Brazil project. The time-interval window was close to 5 years and was sufficient to measure the long-term aggregated dynamics of tree mortality. The area covered by lidar in the 2012–2018 period was ~ 43 km2, ranging from 480 ha at TAL site to 1200 ha at DUC site (Supplementary Table S5).
In addition to the multi-temporal datasets, 610 single-date airborne discrete-return lidar data strips (approx. 300 m wide by 12.5 km long; ~ 3.75 km2 each) were acquired during 2016 (acquisition dates in Supplementary Figure S6A) using the Trimble HARRIER 68i system at an airplane. The average flight height was 600 m above ground and the scan angle was 45° (dataset from the EBA project31).
For both lidar datasets, multiple lidar returns were recorded with a minimum point density of 4 points m−2. Horizontal and vertical accuracy ranging from 0.035 to 0.185 m and from 0.07 to 0.33 m, respectively.
Following the procedures described by Dalagnol et al.19, the lidar point clouds were processed into canopy height models (CHM) of 1-m spatial resolution. The steps of CHM processing included the: (a) classification of the points between ground and vegetation using the lasground, lasheight, and lasclassify functions from the LAStools 3.1.145; (b) creation of a Digital Terrain Model (DTM) using the TINSurfaceCreate function from FUSION/LDV 3.646; (c) normalization of the point cloud height to height above ground using the DTM; and (d) CHM generation by extracting the highest height of vegetation using the CanopyModel function from FUSION.
Environmental and climate data
To analyze the environmental and climatic drivers of gap dynamics, we used a spatialized set of variables and products for the whole Amazon, including: (a) HAND product at 30 × 30 m47; (b) slope calculated from the Shuttle Radar Topography Mission (SRTM) at 30 × 30 m48; (c) soil fertility proxied by SCC at 11 × 11 km37; (d) floodplain cover map at 30 × 30 m49; (e) forest degradation proxied by a non-forest distance map derived from the 30-m global forest change dataset v1.4 (2000–2016)50; (f) monthly mean rainfall (mm), climate water deficit (mm) and wind speed (m s−1), obtained from the TerraClimate dataset at 5 × 5 km (1958–2015)43; and (g) DSL at 28 × 28 km51. All variables and products, except HAND and slope, were resampled to the predominant spatial resolution of most datasets (5 km × 5 km), especially the climate data. We used the SRTM instead of the lidar DTM because the very narrow lidar DTMs (300–500 m) would not permit to determining the lowest point in the terrain to accurately calculate the HAND for every pixel.
Long-term field inventory data
We used data from 181 long-term field inventory plots from the RAINFOR network (Fig. 6)5. The data were collected at closed canopy mixed forests with vegetation structure preserved from fire and logging. All trees with diameter at breast height (DBH) ≥ 10 cm were measured at least twice5. These plots had 852 censuses from 1975 to 2013 with median plot size of 1 ha. The mean re-census interval was 3 years. Tree stem mortality rates (m; % year−1) were calculated as the coefficient of exponential mortality for each census interval and each plot52 (Eq. 1). The m estimates were then averaged by plot and were weighted by the censuses interval length, in years1.
$$m = left[ {lnleft( {N0} right) – lnleft( {Nt} right)} right]/t$$
(1)
where N0 and Nt are the initial and final number of trees, and t is the censuses interval.
Data analysis
Gap definition and static–dynamic gaps relationship
Dynamic gaps were detected using multi-date lidar data at the five study sites: DUC, TAP, FN1, BON, and TAL. We define here dynamic gaps as those opened between two periods of observation associated with canopy turnover events, including tree mortality. For this purpose, we calculated a delta height difference of 10 m between the two acquisitions (~ 5 years apart) and filtered for detections with area greater than 10 m2. This height difference was strongly correlated with tree loss at the canopy level in previous studies19, 20. Because standing dead trees do not necessarily generate gaps, we assume that the dynamic gaps are mostly related to the felled canopy trees associated with broken and uprooted mode of death.
Static gaps were delineated using the CHM from the second lidar acquisition at the five sites (Supplementary Material S1). We applied and compared two types of gap delineation: a traditional method based on a fixed height cutoff (H = 2, 5 or 10 m), and an alternative method based on the relative height (RH = 33, 50, and 66% maximum tree height) around a neighborhood (W = 5–45 m). Since the relative height method did not depend on absolute height values, it should better account for local canopy variability and lower stature vegetation, as opposed to the fixed height method. For both methods, we tested a variety of parameters in the search of an optimal calibration amongst the sites. We filtered gaps for a minimum area of 10 m2, which corresponded to an approximation of the mean canopy area of trees greater than 5-cm DBH in tropical forests21. We also filtered them for a maximum area of 1 ha to automatically exclude open areas that very likely did not correspond to small-scale disturbance from treefall gaps21.
The spatial match between each static and dynamic gap event was assessed by intersecting the detections and calculating metrics of precision (p), recall (r) and F1-score (F) (Eqs. 2–4) (more information at Supplementary Material S1). p represents the percentage of total correct detections, r represents the percentage of reference data correctly mapped, and F represents the harmonic mean between p and r, that is, a balance between commission and omission errors. Methods and parameters were compared to determine the optimal method for static gap delineation, i.e. higher F means greater agreement between static and dynamic gaps.
$$Precisionleft( p right) = true , positives/number , of , gap , polygons$$
(2)
$$Recallleft( r right) = true , positives/number , of , mortality , polygons$$
(3)
$$F1 – scoreleft( F right) = left( {2 times p times r} right)/left( {p + r} right)$$
(4)
Finally, considering the optimal gap delineation method, we modeled the relationship between static-dynamic gaps at the landscape scale using a linear regression. For this purpose, annualized dynamic gap fraction and static gap fraction (i.e., the area occupied by gaps in relation to the total area of the flight line) were calculated at the 5-ha scale. Following the strategy by Wagner et al.53, we defined this value after several simulation tests between variable estimates, change rates and plot area (Supplementary Figure S7). Data and residuals were inspected for normality, and variables were transformed to the logarithmic scale prior to the linear model fitting. To assess the model, we calculated the coefficient of determination (R2), absolute Root Mean Square Error (RMSE) and relative RMSE (%) (ratio of RMSE and the mean of observations). To obtain more reliable and unbiased estimates of the model predictive performance, we calculated the RMSE considering out-of-sample values with a leave-one-site-out cross-validation (CV) strategy. Thus, we fitted the model with four sites and calculated the RMSE with predicted and observed values for the site not used in the modeling. We repeated this process for all five sites. A 95% prediction interval described the variability of tree mortality estimates from the gap fraction.
Spatial variability of static gaps across the Brazilian Amazon
We delineated static gaps on the single-date airborne lidar datasets (n = 610 flight lines) using the optimal gap delineation method and parameters assessed in the previous section. To characterize the gaps variability across the region, we calculated the gap fraction and mean gap size for each site.
Assessment of landscape- and regional-scale drivers of static canopy gaps
To quantify the relationship between static gaps and landscape- and regional-scale predictors, we employed correlation matrices and generalized linear models (GLM). Binomial GLM and Gaussian GLM were applied for landscape and regional models, respectively (detailed information at Supplementary Material S2). Models were assessed using a tenfold CV approach with 30 repetitions. The gap data used in this analysis were those obtained from the 610 single-date lidar data. We defined landscape-scale drivers as those showing great heterogeneity intra-site such as the topography (HAND and slope variables). We defined regional-scale drivers as those having great variability across sites such as the rainfall (Mean_pr and SD_pr), wind speed (Mean_vs and SD_vs), climate water deficit (Mean_def and SD_def), DSL, SCC, floodplains, and non-forest distance.
Through the modeling we evaluated if gap occurrence (presence or absence) and gap size increased at valleys and steep terrains of the Amazon, represented by low HAND and high slope, respectively. As previously demonstrated with tree mortality ground observations, we also tested if gap fraction would increase with: (1) higher water stress, represented by low Mean_pr, and high SD_pr, Mean_def, SD_def, and DSL; (2) higher soil fertility, expressed by high SCC; (3) higher wind speed, proxied by high Mean_vs and SD_vs; (4) higher forest degradation/fragmentation, represented by low non-forest distance; and (5) areas of seasonally flooded forests, expressed by high floodplains cover. Model residuals were inspected in comparison to fitted values using also variogram and Moran’s I analyses to assess for potential biases and spatial correlation (detailed information in Supplementary Material S2). Static gap fraction and Nonforest_dist were transformed to log-scale due to non-normality data.
Amazon-wide dynamic gaps mapping and relationship with tree mortality
To obtain a map of dynamic gap estimates over the Amazon, we first applied the GLM model based on environmental and climate drivers to estimate static gap fractions for the whole region. We then applied the static–dynamic gaps relationship to estimate annualized dynamic gap fraction (% year−1). To explore the opportunities for tree mortality estimates based on gap dynamics, we compared the spatialized dynamic gap estimates with time-averaged mortality rates from long-term field inventory data using a linear model. The model was assessed using a tenfold CV approach with 30 repetitions and the RMSE calculated out-of-sample. We acknowledge that the comparison between field tree mortality and lidar gap estimates is not trivial. However, it is the best source available of independent mortality data to compare the results. Field plot-estimates located within the same 5-km cell of the lidar gap estimates were averaged, resulting in 88 pairs of lidar- and field-estimates samples for validation. The mean annualized dynamic gap fraction per Amazonian region (Fig. 6) was extracted and compared using one-way ANOVA and post-hoc Tukey–Kramer tests. More
