Impact of a tropical forest blowdown on aboveground carbon balance

Field measurements demonstrate a carbon sink in the Amazon and Congo basins, but the cause of this sink is uncertain. One possibility is that forest landscapes are experiencing transient recovery from previous disturbance. Attributing the carbon sink to transient recovery or other processes is challenging because we do not understand the sensitivity of conventional remote sensing methods to changes in aboveground carbon density (ACD) caused by disturbance events. Here we use ultra-high-density drone lidar to quantify the impact of a blowdown disturbance on ACD in a lowland rain forest in Costa Rica. We show that the blowdown decreased ACD by at least 17.6%, increased the number of canopy gaps, and altered the gap size-frequency distribution. Analyses of a canopy-height transition matrix indicate departure from steady-state conditions. This event will initiate a transient sink requiring an estimated 24–49 years to recover pre-disturbance ACD. Our results suggest that blowdowns of this magnitude and extent can remain undetected by conventional satellite optical imagery but are likely to alter ACD decades after they occur.

; Supplementary Fig. 2). By comparing these measurements to previously collected airborne lidar and inventory plot data, we identified multiple lines of evidence that indicate the disturbance was unprecedented in the previous 20 years of annual field measurements 20 . These measurements are used to address three main questions: 1) what is the impact of the blowdown disturbance on ACD and forest structure? 2) what is the estimated recovery time of the forest to its pre-disturbance ACD?
3) to what extent can similar events be detected by conventional remote sensing observations for large scale monitoring applications?

Results and discussion
Impact of disturbance on ACD and forest structure. We found that this blowdown disturbance caused substantial losses to ACD and changes to forest structure. The blowdown decreased mean ACD by 17.6% in comparison to pre-blowdown conditions, and average ACD loss was similar for old-growth (18.3%) and secondary forest areas (17.2%) in our study area (Supplementary Table 1). Carbon losses were heterogeneous throughout the forested landscape-some locations decreased in ACD by up to 63.1% while locations that were not impacted by the blowdown increased by up to 19.8% ( Supplementary Fig. 1). We explored two potential causes of this spatial heterogeneity. First, we found no evidence that ACD loss was greater along trails ( Supplementary  Fig. 3), which is inconsistent with expectations that tree damage and mortality are greatest along forest edges 21 . This discrepancy is possible because local trails are likely not wide enough to increase wind exposure. Second, we found a weak correlation between ACD prior to the blowdown and ACD loss ( Supplementary Fig. 4), which is consistent with previous research that found larger trees are more susceptible to damage in blowdowns 22 . The blowdown also significantly increased the size and frequency of canopy gaps, more than doubling the total forested area in gaps extending to ≤ 8 m height (Fig. 2a). The proportional increase in gap area was greater in secondary forests (arithmetic mean 82% increase) than in old-growth forest (mean 74% increase). The canopy gap size-frequency distribution was well described by a power-law probability distribution both before and after the blowdown, where most canopy gaps were small in size ( Supplementary Fig. 5). The power-law scaling exponent ( ) of the gap size-frequency distribution was consistently smaller after the blowdown for all gap height thresholds between 2 and 20 m, indicating an increased relative prevalence of larger gaps after the blowdown event. However, the decrease in was only significantly different for mid-canopy gaps extending to between 4 and 12 m height (Fig. 2b). Surprisingly, the size-frequency distribution of classically defined canopy gaps ≤ 2 m in height did not change significantly in response to the blowdown event 23 . Although classically defined gaps are useful because they can be easily quantified in ground surveys, our findings underscore the importance of high-resolution three-dimensional information to characterize the impacts of disturbance events in forest canopies (Fig. 1). Gap size-frequency results were qualitatively similar when old-growth and secondary forests were analyzed separately ( Supplementary Fig. 6).
Previous analyses demonstrated that canopy height and ACD changes in the old-growth forest landscape were consistent with steady state expectations 24,25 . We projected the equilibrium canopy height distribution using observed canopy height changes before and after the blowdown. The projected equilibrium is the expected distribution of canopy height if the canopy height transition probability matrix associated with the blowdown event persists indefinitely. By comparing the projected equilibrium distribution to the observed distribution prior to the blowdown event, we can test the hypothesis that the distribution of canopy height change associated with the blowdown event is a departure from steady-state conditions previously reported. Our projection of the   24 . Projected equilibrium heights were qualitatively similar when old-growth and secondary forests in our study area were analyzed separately ( Supplementary Fig. 7). Our interpretation assumes that changes in the distribution of canopy height and ACD between 2009 and 2019 were driven by the 2018 blowdown event. Interval length can bias estimates of carbon stock change 26 . Here we likely underestimate ACD losses and structural changes attributable to the blowdown for two reasons. First, we expect that the landscape accumulated aboveground carbon between 2009 and the 2018 blowdown event. We expect this is true because our study area contains secondary forests, and because field measurements indicate that old-growth forests in this area steadily increased in ACD between 1997 and 2017 20 . We are also likely to underestimate ACD losses from the blowdown because post-blowdown lidar data were collected a year after the disturbance, and initial recovery can be rapid 11,21,27 . Consequently, we consider our conclusion regarding carbon loss and increased gap area to be conservative.  www.nature.com/scientificreports/ Estimated duration of ACD recovery following disturbance. We used existing information about forest dynamics at this site to quantify the duration of the carbon sink (from forest regrowth) expected from this blowdown disturbance. We estimated recovery time by combining data from long-term field measurements at this site describing two processes: transient, rapid recovery following localized carbon losses (treefalls, Supplementary Fig. 8) and the long-term growth rate of ACD 20,28 . We estimate that carbon losses associated with this event will require 24-35 years to recover pre-disturbance mean ACD. Field measurements indicate an absolute increase in ACD of 0.7 Mg C ha −1 yr −1 between 2009 and 2016. Applying this increase to our estimate of preblowdown ACD increases the magnitude of carbon loss to 24% and increases the estimated recovery time to 38-49 years. Notably, these timeframes are longer than most field records of carbon accumulation in tropical forests. Two recent summaries of carbon balance in Amazonian and African tropical forests reported a mean monitoring period of 11-12 years 1,2 .
Our projected ACD recovery times assume that old-growth and secondary forests will regrow at the same rate. Previous studies of secondary forests in this landscape found that ACD is similar in old-growth forests and secondary forests older than 21 years, with rapid ACD recovery likely facilitated by the relatively fertile soils and proximity to intact forest 29,30 . Tree species richness and composition also stabilized within this time 29,31 . Our study area included only secondary forests older than 21 years at the beginning of our study ( Supplementary Fig. 2); secondary forest ACD was within 93% of old-growth forest ACD at the beginning of our study (Supplementary Table 1). It is possible that secondary forest recovery rates will be slightly faster due to remaining compositional differences, such as lower liana abundance 32 .
Potential for detection with conventional remote sensing methods. Although this disturbance reduced ACD by at least 17.6%, the blowdown was considerably smaller in size than larger reported blowdowns in the Amazon-the largest gaps we report are ~ 1 ha in size ( Supplementary Fig. 5) while other studies report single gaps over 10 times as large 17 . We asked whether conventional optical satellite methods used to characterize larger blowdowns are also useful for detecting more moderate blowdowns. To determine whether the event was detectible using conventional optical satellite data, we identified cloud-free Landsat observations before and after the blowdown and computed the change in non-photosynthetic vegetation (ΔNPV) in each pixel. ΔNPV is used to detect blowdowns because tree mortality increases the fraction of woody vegetation exposed to the sensor 9,11,21,22 . As expected, ΔNPV values were significantly positive-i.e. the fraction of NPV increased after the blowdown (t = 16.2, DF = 206, P < 0.001). However, the mean and maximum ΔNPV in our data (0.02 and 0.09, respectively) were smaller than ΔNPV values reported in other studies-in fact, these values are approximately one order of magnitude smaller than those reported for large blowdowns 9 , and the mean value is within the 95% confidence interval for no mortality in one previous study 21 . Further, although we observed a statistically significant relationship between ΔNPV and ACD loss, the explanatory power of this relationship was very low (r 2 = 0.03) ( Supplementary Fig. 9). We attribute the difficulty of detecting this blowdown using Landsat to a combination of factors. Relatively low temporal frequency associated with high cloud cover resulted in cloud-free observations 6 and 7 months before and after the blowdown event, respectively. Tropical forests regrow quickly after blowdown events-even severe blowdowns are not detectable after < 2 years 11 .
Prospective importance for tropical forest ACD dynamics. Results from this study are consistent with the idea that increasing ACD in tropical forest plots could, in part, reflect recovery from natural disturbance events that predate the onset of forest monitoring and remain cryptic to historical satellite records. This explanation does not preclude other potential causes of increasing ACD, such as carbon fertilization from rising atmospheric CO 2 and widespread prior anthropogenic disturbance. The magnitude of ACD loss and estimated ACD recovery reported here are specific to this event and not directly generalizable to other areas because blowdown intensity and recovery from disturbance are both highly variable. Instead, this case study is valuable for highlighting that moderate episodic events can cause ACD losses of at least 17.6% while remaining undetectable using conventional measurement methods.
We anticipate that advances in satellite remote sensing-like Planet's constellation of cube satellites and NASA's Global Ecosystems Dynamics Investigation (GEDI) spaceborne lidar-will provide new opportunities to thoroughly characterize how other moderate blowdowns contribute to carbon cycling in global tropical forests [33][34][35] . For instance, mortality from the blowdown described here was clearly apparent in high-resolution (3 m) Planet data collected 8 days before and seven days after the blowdown event ( Supplementary Fig. 1). The methodological framework presented here, in which a targeted airborne remote sensing campaign produced detailed measurements of ACD loss and structural changes, will be valuable for developing and benchmarking spaceborne measurements. Airborne lidar are increasingly available-such as > 100,000 ha of available tropical forest lidar data from the Sustainable Landscapes Brazil Project 36,37 , including multi-temporal lidar for some areas 38 -and will provide essential baseline data for future efforts. Our results suggest that incorporating widespread high-resolution monitoring of forest disturbance could increase the estimated aboveground carbon loss attributed to infrequent disturbance events in comparison to previous analyses limited to field records or historical satellite data 15 .

Methods
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 39  www.nature.com/scientificreports/ 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 40 . We focused on an area 1.4 km 2 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 40 , we limited our analysis to lidar returns with scan angle ± 15 degrees and randomly subsampled data to a homogenous resolution of 10 pts m −2 . Previous research demonstrates that lidar data collected above densities of 1 pts m −2 have similar predictive power for determining many forest properties (including tree height, tree density, and basal area) 41 ; 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 2 = 0.994, RMSE = 1.85 m; Supplementary Table 2) 42 . 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 43 . 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 44 . 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 45 -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 28 .
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 46 . 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 47 . 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 48 . To convert AGBD predictions into ACD, we multiplied AGBD by 0.47. 49 . 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 10 . 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 50 . We defined gaps according to Brokaw's classic definition: any contiguous area ≤ 2 m in height 23 , 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: www.nature.com/scientificreports/ where ζ ( ) is the Riemann zeta function. The parameter is a power-law exponent describing the sizefrequency distribution of gaps-small values of 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 10 .
We estimated using the Metropolis-Hastings algorithm within a Markov chain Monte Carlo (MCMC) procedure. Our analysis produces a Bayesian estimate of . We used an uninformative prior (uniform distribution between 1.01 and 5); this conservative potential range for 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 and standard deviation equal to 0.1; proposals of < 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, . 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 24 . 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 x is the right-hand eigenvector, and is the eigenvalue (not to be confused with the power-law exponent used to characterize the gap size frequency distribution). Here, 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 51 . 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 52 . 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 43 . 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 20 using the allometry of 55 . Here we report the long-term annual gain of 0.49 Mg C ha −1 yr −1 using the allometry of 44  www.nature.com/scientificreports/ 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 −1 . 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 20 . We calculated that the average annual ACD gain between 2009 and 2016 was 0.74 Mg C ha −1 yr −1 . 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.

Data availability
Canopy height data and code underlying this analysis are available at https:// github. com/ kccus hman/ LaSel vaBlo wdown; the repository contents at the time of submission are also permanently archived through the Brown University Dataverse (hosted by the Harvard Dataverse) at https:// doi. org/ 10. 7910/ DVN/ VHW3XR.