Impact of a tropical forest blowdown on aboveground carbon balance

Study site

This study was conducted at La Selva Biological Station, located in the lowland Atlantic forest of Costa Rica (10°26′ N, 83°59′ W). The mean annual temperature is 26 °C; mean annual precipitation is 4 m and all months have mean precipitation > 100 mm³⁹. La Selva has undulating topography, with elevation varying between 10 and 140 m above sea level. La Selva Biological Station includes multiple land uses; our analysis includes 103.5 hectares of forest, comprising 33.0 ha of old-growth forest and 70.5 ha of forests with past human disturbance (secondary forests, abandoned agroforestry, abandoned plantation, selectively-logged forests); here, we refer to all areas with past human disturbance as “secondary forests”. This study area does not include the full extent of old-growth or secondary forests at La Selva—we focused our drone data collection on this area because it contained the most severe apparent disturbance from the blowdown. Forests with past human disturbance have been naturally regenerating for a range of time (since 1955–1988); we excluded secondary forests with regeneration starting after 1988.

Lidar data

We use two airborne lidar datasets to quantify dynamics in canopy structure and ACD. Data were collected in 2009 and 2019 (Supplementary Table 2). Data from 2009 were collected by a fixed-wing aircraft over the entire reserve; data from 2019 were collected using the Brown Platform for Autonomous Remote Sensing⁴⁰. We focused on an area 1.4 km² in size that includes the region of most severe damage from the blowdown (Supplementary Fig. 1). Both lidar sensors were discrete-return systems. To minimize variation in lidar height estimates from variable laser beam divergence and detector characteristics, we only used data from first returns for all analyses. For the 2019 drone-based lidar with higher native point density and a wider scan angle range⁴⁰, we limited our analysis to lidar returns with scan angle ± 15 degrees and randomly subsampled data to a homogenous resolution of 10 pts m⁻². Previous research demonstrates that lidar data collected above densities of 1 pts m⁻² have similar predictive power for determining many forest properties (including tree height, tree density, and basal area)⁴¹; both lidar datasets in this study are above this density threshold. All lidar data were projected using EPSG 32,616.

For all lidar data, we calculated height above ground using a digital terrain model (DTM) created from lidar data collected in 2006 and validated using 4184 independent measurements within the old-growth forest (intercept =  − 0.406, slope = 0.999, r² = 0.994, RMSE = 1.85 m; Supplementary Table 2)⁴². We verified that the horizontal geolocation accuracy with < 1 m between lidar datasets by comparing lidar returns from building roofs present in all datasets; we also used data from roof lines to adjust for a 0.7 m positive bias in the 2009 measurements relative to the 2019 lidar measurements (Supplementary Table 2).

AGBD and ACD

We estimated aboveground biomass density (AGBD) for each lidar dataset using a model parameterized with 18 0.5 ha field plots established for the CARBONO project, in which plots were placed using a random stratified sampling design to characterize edaphic variation within La Selva⁴³. All trees, lianas, and palms with diameter ≥ 10 cm at 1.3 m height are measured annually in CARBONO plots; field data from 2009 were used to parametrize the lidar-derived AGBD model. For each field plot, aboveground biomass was estimated for each stem using an allometric model including a regional diameter-height relationship and wood density⁴⁴. We used wood density values at the most specific level possible (species, genus, family, or site-level mean). AGBD estimates do not include non-woody plant material, herbaceous plants, hemi-epiphytes, epiphytes, or woody stems smaller than 10 cm diameter. Omitted pools likely comprise < 15% of total AGBD⁴⁵—the proportion of omitted pools may be larger for secondary forests, but we do not expect this to be an important source of bias because our study area includes only older secondary forests. Data and detailed methods for plot measurements are publicly available²⁸.

We used a model relating top-of-canopy height (TCH) to AGBD using a power relationship:

where a and b are parameters fit using non-linear maximum likelihood analysis⁴⁶. Maximum likelihood analysis was performed using the quasi-Newton method in the ‘stats4’ package in R (version 4.0.3), and assuming normally distributed residuals⁴⁷. TCH was calculated by first computing the mean height of first lidar returns in every pixel of a 5 m × 5 m grid, and then computing the mean height of all pixels that fell within the boundaries of a single plot. This model explained 74% of variation in AGBD among field plots, with 9.2% RMSE (Supplementary Fig. 10). We assume that the relationship between TCH and AGBD was the same before and after the blowdown event. The distributions of model residuals showed no heteroscedasticity (Supplementary Fig. 10). Previous research indicates that a single relationship fits old-growth and secondary forests at La Selva⁴⁸. To convert AGBD predictions into ACD, we multiplied AGBD by 0.47.⁴⁹.

We applied this model to the 2009 and 2019 lidar datasets, using a 0.5 ha raster resolution corresponding to the field plot size. We used a Monte Carlo method, sampling 1000 times from the observed parameter values describing residual variation in the TCH to AGBD model, to quantify uncertainty (95% confidence intervals) in estimates of ACD and ACD change.

Gap size-frequency distribution

We quantified the canopy gap size-frequency distribution for each lidar dataset by creating a canopy height model (CHM) with 1.25 m pixels¹⁰. To ensure that every pixel had a height value (in the occasional case where a pixel had no lidar returns), we created the canopy height model by using a Delaunay triangulation of first returns, gridded to 1.25 m resolution⁵⁰. We defined gaps according to Brokaw’s classic definition: any contiguous area ≤ 2 m in height²³, and we also repeated the analysis for a range of height thresholds up to 20 m height (in 2 m increments). We included diagonal pixels in our calculation of contiguous area.

We characterized the gap size-frequency distribution using the Zeta distribution, which is a discrete probability distribution, defined for integers (k) ≥ 1, giving the probability that a gap contains k pixels:

where (zeta left(lambda right)) is the Riemann zeta function. The parameter (lambda ) is a power-law exponent describing the size-frequency distribution of gaps—small values of (lambda ) indicate increased frequency of large gaps. Previous research indicates that the power-law Zeta distribution is appropriate for comparing gap sizes in tropical forests with diverse disturbance regimes¹⁰.

We estimated (lambda ) using the Metropolis–Hastings algorithm within a Markov chain Monte Carlo (MCMC) procedure. Our analysis produces a Bayesian estimate of (lambda ). We used an uninformative prior (uniform distribution between 1.01 and 5); this conservative potential range for (lambda ) was chosen based on results from Kellner and Asner (2009). We used a random normal proposal density, with a mean equal to the previous iteration of (lambda ) and standard deviation equal to 0.1; proposals of (lambda ) < 1 were replaced with a new random proposal. We used a chain of length 100,000, discarded the first 5000 observations the burn-in period, and thinned the chain by using every 25^th value. We used the remaining 3800 values to determine the posterior median and 95% Bayesian credible intervals for the power-law exponent, (uplambda ). Results from our entire study area (old-growth and secondary forests combined) are reported in the main text; we also repeated this analysis for old-growth and secondary forest areas separately (Supplementary Fig. 6).

Canopy height change

We quantified the distribution of canopy height change between 2009 and 2019 using the 5 m resolution CHM, where canopy height was estimated using the mean height of all first lidar returns in a pixel; at this resolution, there were no pixels with no lidar returns. A forest in steady-state is expected to have mean canopy height change of approximately zero. A forest recovering from past disturbance is expected to have a mean canopy height change greater than zero, and a forest that experiences large disturbance during the interval between measurements is expected to have a mean canopy height change less than zero²⁴. We calculated the distribution of canopy height change between 2009 and 2019 by subtracting the initial height of a pixel from the final height of a pixel (Supplementary Fig. 11).

The steady-state canopy height distribution of a forest is the expected distribution of canopy height if observed transition probabilities do not change. To calculate the projected steady-state canopy height distributions for each time interval, we created a canopy height transition matrix, (A), with 54 rows and 54 columns. Each row and column of (A) corresponds to a single 1-m height class, and the maximum value (54 m) was selected from the maximum height of any pixel. An entry from row i and column j in (A), a_ij, represents the number of pixels that were in height class j at the beginning of the time interval and in height class i at the end of the time interval. The projected steady-state height distribution is then obtained using an eigenvector decomposition:

where (mathbf{x}) is the right-hand eigenvector, and ({varvec{uplambda}}) is the eigenvalue (not to be confused with the power-law exponent used to characterize the gap size frequency distribution). Here, (mathbf{x}) gives the distribution of canopy heights for which applying the canopy height transition matrix results in no overall change in the distribution of heights. We used a Bayesian procedure to quantify uncertainty in our projections of steady-state canopy heights. Specifically, we assume that height transition probabilities (i.e. the columns of (A)) follow a multinomial distribution. The multinomial distribution has a conjugate prior distribution, the Dirichlet distribution, from which the posterior distribution can be solved numerically⁵¹. We assumed an uninformative Dirichlet prior, which is equivalent to assuming that height transitions are uniformly distributed across prospective height classes. We sampled from the posterior distribution of each height class 10,000 times to determine the 95% Bayesian credible intervals of the projected steady-state canopy height distribution. Results from our entire study area (old-growth and secondary forests combined) are reported in the main text; we also repeated this analysis for old-growth and secondary forest areas separately (Supplementary Fig. 7).

Blowdown detection from Landsat imagery

To assess whether this blowdown event could be detected using Landsat imagery, we compared Landsat 8 images before and after the blowdown event. We manually chose Landsat images that were closest in time to the event without cloud cover over our area of interest. These images came from November 10, 2017 (190 days before the blowdown) and December 31, 2018 (226 days after the blowdown). We used the Landsat 8 Surface Reflectance Tier 1 data product, which is atmospherically corrected using United States Geological Survey Land Surface Reflectance Code. We then performed a Spectral Mixture Analysis (SMA) with endmembers for photosynthetic vegetation, non-photosynthetic vegetation (NPV), and shade; previous studies have shown that the change in proportion of NPV per pixel correlates with blowdown mortality and tree damage^9,17,21,22. Because no pure pixels of NPV were apparent in our image, we used the tropical forest Landsat endmembers published by Schwartz et al. (2017). SMA was performed using ENVI’s linear spectral unmixing tool, using the constraint that endmembers must sum to one⁵². We compared the change in NPV to the change in ACD across our study landscape (Supplementary Fig. 9).

Estimated recovery time from disturbance

We estimated recovery time required for the landscape to recover to its pre-blowdown ACD using the CARBONO project plot data⁴³. Tropical forest AGBD recovery following disturbance is highly variable and depends on multiple factors including the severity of disturbance, climate, and soils^53,54. Any approach to estimate recovery time will be highly uncertain due to unknowable factors such as future disturbance and climate, so we used two approaches that are likely to over- and under-estimate recovery time to approximate the range of likely recovery times.

First, we estimated recovery time using the long-term annual ACD gain in old-growth forests at our study site from 1997 to 2016. This long-term gain was previously reported in²⁰ using the allometry of⁵⁵. Here we report the long-term annual gain of 0.49 Mg C ha⁻¹ yr⁻¹ using the allometry of⁴⁴. Because our analysis indicates a mean ACD loss of 17.4 Mg C ha⁻¹ (Supplementary Table 1), the long-term ACD trend predicts 35 years to recover from the blowdown. We expect that this overestimates recovery time because ACD gain directly after the disturbance event should be greater than the long-term average.

Second, we estimate recovery time by incorporating fast initial ACD gain following large treefalls. We calculated the average ACD gain following the four observations in the CARBONO record where the annual ACD loss was > 11.5 Mg C ha⁻¹. For 5 years following these disturbance events, annual ACD gain was larger in magnitude than the long-term trend (Supplementary Fig. 8). We estimated recovery time assuming these elevated post-disturbance recovery rates for the first 5 years, and then assuming the long-term gain for subsequent years, resulting in a predicted recovery time of 24 years. We expect that this underestimates recovery time because not all parts of the landscape experienced large decreases in ACD, but this method assumes that the entire landscape will have elevated ACD gain following the disturbance.

Finally, we repeated both approaches above, additionally estimating a higher value for pre-blowdown ACD in 2018 because plot based observations at La Selva indicate that ACD increased between 2009 and 2016²⁰. We calculated that the average annual ACD gain between 2009 and 2016 was 0.74 Mg C ha⁻¹ yr⁻¹. To estimate ACD in 2018 immediately prior to the blowdown, we added this average rate to our 2009 ACD estimates for the period of 2009–2018. We subsequently re-estimated recovery time to the higher 2018 ACD estimate, resulting in recovery time between 38 and 48 years using the second and first approaches above, respectively.

Source: Ecology - nature.com