Key drivers of ecosystem recovery after disturbance in a neotropical forest

Background

In tropical forests, natural disturbances caused by the death of one or more trees are the dominant forms of forest regeneration as the creation of canopy openings continuously reshapes forest structure (Goulamoussène et al. 2017). The immediate increase in light intensity allows the sunlight to penetrate the understorey (Goulamoussène et al. 2016) and light-demanding trees (Denslow et al. 1998) to establish and grow, thus contributing to the maintenance of biodiversity that shapes forest functionning (Liang et al. 2016). Another effect of canopy gaps is the local modification of the forest nutrient balance due to the large amounts of dead leaves and wood that decompose and mineralize (Brokaw and Busing 2000) and that shapes in turn the small-scale spatial variations in forest carbon balance (Feeley et al. 2007; Guitet et al. 2015; Rutishauser et al. 2010). In this way, the natural disturbance regime is a fundamental component of the functioning of tropical forests (Sheil and Burslem 2003).

With ongoing global (i.e land-use and climate) changes, tropical forests are currently facing deep and rapid changes in disturbance regimes that may hamper their recovering capacity (Hérault and Gourlet-Fleury 2016; Brienen et al. 2015). Human-induced disturbances may encompass a wide range of perturbations from long-lasting ones such as land-cover changes for industrial agriculture, slash-and-burn agriculture or mining (Dezécache et al. 2017a, b) to more insiduous modifications such as selective logging that may not affect the forest cover but modify forest functioning (Rutishauser et al. 2015). An even more insiduous perturbation is climate change (Hérault and Gourlet-Fleury 2016). Global circulation models have shown high probabilities of significant precipitation decrease for tropical areas with a risk of transition from short-dry-season rainforest to long-dry- season savannah ecosystems (Davidson et al. 2012). For instance, after the intense 2005 drought in Amazonia, the forest suffered an additional mortality, leading to a huge loss of live biomass (Phillips et al. 2009) with similar mortality events observed in Panama (Condit 1995), in China (Tan et al. 2013) or in South-East Asia (Slik 2004).

To our opinion, the drivers of the post-disturbance system trajectory may first be defined based on their origin: endogeneous and exogeneous. (1) Endogenous drivers refer to the internal properties of the system that may influence its post-disturbance behavior. A significant example for ecological systems is the species composition that partially informs on the immediate potential of the system to recover after disturbance. In that respect, the species identity is far less important than the functional signature of the species assemblage (Kunstler et al. 2016): for example, an assemblage of light-demanding species will respond differently to disturbance from an assemblage of shade-tolerant understorey species (Herault et al. 2010). Structural characteristics of the pre-disturbance species community (stem density, average size, live biomass and so on) may also be of primary importance because they are core indicators of the silvigenetic stage of the forest (Pillet et al. 2017). (2) Exogenous drivers refer to external constraints or forces that limit the possible system trajectories. They can be grouped into two broad categories: drivers that vary in space and those that vary in time. The local environment, i.e. the physical characteristics of the abiotic environment, is here defined in space but not in time. On the contrary, external conditions such as climatic stress are here considered to vary in time but not in space.

This study draws upon the long-term disturbance experiment of Paracou, French Guiana, to develop a modeling approach in order to mechanistically link the endo- and exogenous ecosystem drivers to the ecosystem recovery trajectory after disturbance. More specifically, we ask the following questions: (ⅰ) Do regenerating forests recover faster than mature forests given the same level of disturbance? (ⅱ) Is the local topography an important predictor of the forest recovery rates? (ⅲ) Is the community functional composition, assessed with community weighted-mean functional traits, a good predictor of carbon stock recovery? (ⅳ) How important is the climate stress (drought and/or soil water saturation) to shape the rate of carbon recovery? To do so, we partition the contributions to post-disturbance ACS (Aboveground Carbon Stock) gain (from growth and recruitment of trees sup 10 cm DBH) and ACS loss (from mortality) of survivors and recruited trees to detect the main drivers and patterns of ACS recovery after disturbance. We model the trajectory of those post-disturbance ACS changes (Piponiot et al. 2016b) in a comprehensive Bayesian framework. We then quantify the effect of (ⅰ) endogeneous (forest structure and composition) and (ⅱ) exogeneous (local environment and climate stress) drivers on the rates at which post-disturbance ACS changes converge to a theoretical steady state. Summing these ACS changes over time gives the net post-disturbance rate of ACS accumulation, an indicator of the ecosystem recovery rate. Disentangling ACS recovery with a demographic process-based approach, i.e. by segregating ACS changes into cohorts (survivors and recruits) and demographic processes (growth, recruitment, mortality), as opposed to an all-in-one model in which only the ecosystem net ACS change is modeled without examination of demographic processes, has been shown to be essential to reveal mechanisms underlying ACS responses to disturbance and to make more robust predictions of ACS recovery (Piponiot et al. 2016b).

Methods

Study site

The study was conducted at the Paracou experimental site (5°18'N, 52°55'W), a lowland tropical rain forest near Sinnamary, French Guiana. The site receives nearly two-thirds of the annual 3041 mm of precipitation between mid-March and mid-June, and < 50 mm per month in September and October (Wagner et al. 2011). More than 700 woody species attaining 2 cm DBH (diameter at breast height) have been described at the site, with 150 - 210 species of trees > 10 cm DBH per hectare. The floristic composition is typical of Guianan rainforests with dominant families including Leguminoseae, Chrysobalanaceae, Lecythidaceae, Sapotaceae and Burseraceae (Guitet et al. 2014). In 1984, nine 6.25 ha plots, each one divided into 4 subplots of 1.56 ha each, were established for a complete inventory of all trees > 10 cm DBH. From October 1986 to May 1987, the plots underwent three disturbance treatments (details in Table 1 and in (Blanc et al. 2009)).

Table  1.  Disturbance treatments (T1, T2, T3) implemented on the Paracou plots in 1986-1987

	Timber logging	Fuelwood logging	Thinning	% ACS loss
T1	DBH ≥ 50 cm, mean of 10 trees ·ha−1	-	-	[12−33%]
T2	DBH ≥ 50 cm, mean of 10 trees ·ha−1	-	DBH ≥ 40cm, all non-valuable species, mean of 30 trees ·ha−1	[33−56%]
T3	DBH ≥ 50 cm, mean of 10 trees ·ha−1	40 cm ≤ DBH ≤ 50 cm, all non-valuable species, mean of 20 trees ·ha−1	DBH ≥ 40cm, all non-valuable species, mean of 15 trees ·ha−1	[35−66%]
The percentage of Aboveground Carbon Stock loss (% ACS loss) is defined as the difference between the pre-disturbance ACS and its minimum value reached during the 4 years after the disturbance treatments

 | Show Table

DownLoad: CSV

Input data

Aboveground Carbon Stock (ACS) computation In all plots, diameter at breast height (DBH) of trees > 10 cm DBH were measured every two years from 1982 to 2016 resulting in 18 forest censuses. Trees were identified to the lowest taxonomic level. To get wood density, we applied the following standardized protocol: (ⅰ) tree identified to the species level were assigned the corresponding wood specific gravity value from the Global Wood Density Database (GWDD) (Chave et al. 2009); (ⅱ) trees identified to the genus level were assigned a genus-average wood density and (ⅲ) trees with no botanical identification or that were not in the GWDD were assigned the subplot-average wood density. The aboveground biomass (AGB) was estimated taking all uncertainties into account using the BIOMASS package (Réjou-Méchain et al. 2017). Biomass was assumed to be 47% carbon.

Disturbance intensity After disturbance, the subplot's ACS decreases rapidly until it reaches its minimum value acs_min a few years later. This transition point determines the beginning of the recovery period. The difference between the averaged pre-disturbance ACS acs_pre and this post-logging minimum value acs_min reached at time t=t_min defines the disturbance intensity DIST. In other words, the disturbance intensity is defined as the amount of aboveground carbon lost in the forest ecosystem during the first years during and after the disturbance.

Structure drivers The pre-disturbance forest structure was assessed with three variables: the stem density S_N (from 483 to 727 ind ·ha^-1) and the basal area S_BA (from 27 to 36 m ²·ha^-1) of subplot j at t_pre, the year preceding the disturbance experiment.

Environment drivers Three environmental drivers were selected from a preliminary exploratory analysis to represent independent source of variation in the local forest physical conditions: the proportion of bottom-lands E_BOTTOM, the average topographical slopes of the plot E_SLOPE and the standard deviation, i.e the heterogeneity, of the altitudinal distribution E_HETE .

Composition drivers The pre-disturbance forest composition was assessed in a functional trait space to avoid local taxonomic variations in tree assemblages that are of little importance for forest functioning. The four chosen orthogonal traits FT represent key dimensions of the tree functional strategy (Baraloto et al. 2010): wood density T_WD, seed mass T_SEED, specific leaf area T_SLA and maximum diameter T_DBH95 estimated as the 95th percentile of the species DBH distribution in the Guyafor database. The community weighted means of these functional traits were calculated the year preceding the disturbance experiment.

Climate drivers We considered two main sources of climate stress: soil drought C_DROUGHT and soil water saturation C_WATER . These variables were quantified using a water balance model, developed and calibrated in Paracou (Wagner et al. 2011), that was run using precipitation and evapotranspiration as inputs over the 1982-2016 time period. C_DROUGHT was estimated as the number of days with REW, Relative Extractable Water, below 0.4 while the number of days with REW equal to 1, the soil is full of water, defined C_WATER . These two covariates were computed between 2 consecutive censuses and then standardized at a yearly time-step.

Modeling strategy

We define two cohorts of trees. First, recruits are all the trees (> 10 cm DBH) that have been recruited since the perturbation. Trees that, for a given census, first went through the 10 cm DBH are called new recruits. Thereafter, they are called, for the following censuses, recruits and may grow or may eventually die between 2 censuses. Second, survivors are trees that were present in the forest before the disturbance and that survived the disturbance event.

For each subplot j and census k, with t_k the time since the beginning of the recovery period, we thus define 5 ACS changes : new recruits' ACS (Rr_{j, k}) is the ACS of all trees < 10 cm DBH at t_k-1 and ≥ 10 cm DBH at t_k ; recruits' ACS growth (Rg_{j, k}) is the ACS increment of living recruits between t_k-1 and t_k ; recruits' ACS loss (Rl_{j, k}) is the ACS in recruits that die between t_k-1 and t_k ; survivors' ACS growth (Sg_{j, k}) is the ACS increment of living survivors between t_k-1 and t_k ; survivors' ACS loss (Sl_{j, k}) is the ACS of survivors that die between t_k-1 and t_k . ACS changes are subject to large stochastic variation over time: because we are less interested in year-to-year variations than in long-term ACS trajectories, we modeled the cumulative ACS changes over time. Cumulative ACS changes (Mg C ·ha^-1) were defined as follows:

where j is the subplot, k is the census number, t_k the time since t₀ (yr) and Change is the annual ACS change (Mg C ·ha^-1·yr^-1), either recruits' ACS (Rr), recruits' ACS growth (Rg), recruits' ACS loss (Rl), survivors' ACS growth (Sg), or survivors' ACS loss (Sl).

Survivors Survivors' cumulative ACS changes are null at t=0 and have a finite limit, attained once survivors have all died. We modeled survivors' cumulative ACS growth cSg as:

where j is the subplot, p the plot it belongs to, t_k is the time since t₀. $\alpha ^{Sg}_{p}$ is the finite limit of the cumulative ACS change, $\beta ^{Sg}_{j, k}$ the rate at which the cumulative ACS change converges to this limit and (σ^Sg)² the variance of the model. By choosing an exponential kernel, we assume that survivors' ACS growth at t_k is proportional to survivors' ACS growth at t_k-1.

Because of our nested design with subplots j within plots p, we modeled the $\alpha ^{Sg}_{p}$ values with a random plot effect of mean $\alpha _{0}^{Sg}$ and variance $\left (\sigma ^{Sg}_{\alpha }\right)^{2}$:

Parameter $\beta ^{Sg}_{j, k}$ is the rate at which survivors' ACS growth on plot j at time t_k converges to a finite limit after the disturbance: it reflects the response rapidity of survivors' ACS growth to disturbance. Because we are interested in predicting variations in $\beta ^{Sg}_{j, k}$, we expressed the latter as a function of covariates:

with $\beta ^{Sg}_{0}$ the model intercept, $\lambda ^{Sg}_{l}$ the vector of l parameters associated to the covariates V_{j, k, l} for which we looked at their effects on the post-logging rate $\beta ^{Sg}_{j, k}$ in subplot j at time t_k . The covariates are defined above and related to the disturbance intensity DIST, the structure of the forest before disturbance (S_N, S_DG, S_BA ), the functional trait composition of the forest before disturbance (T_SEED, T_SLA, T_WD, T_DBH95), the local environment (E_BOTTOM, E_SLOPE, E_HETE ) and the climate stress (C_DROUGHT, C_WATER ). Note that values of the two later covariates changed with times. All covariates are centered and standardized before the inference. When all survivors in plot p are dead, all the C gained by their growth ${\left ({cSg}_{j, p, \infty }=\alpha ^{Sg}_{p}\right)}$ plus their initial ACS (acsmin_j ) will have been lost ${\left ({cSl}_{j, \infty }=\alpha ^{Sl}_{j}\right)}$. We thus defined

with $\alpha ^{Sg}_{p}, \alpha ^{Sl}_{j}$ the finite limits of survivors' cumulative ACS growth and ACS loss respectively, and acsmin_j the ACS of the subplot j at t_min =t₀. Then the cumulative carbon loss is

where j is the subplot, p the plot it belongs to, t_k is the time since t₀. $\alpha ^{Sl}_{j}$ is the finite limit of the cumulative ACS change, $\beta ^{Sl}_{j, k}$ the rate at which the cumulative ACS change converges to this limit and (σ^Sl)² the variance of the model. And with

with $\beta ^{Sl}_{0}$ the model intercept, $\lambda ^{Sl}_{l}$ the vector of l parameters associated to the covariates V_{j, k, l} for which we looked at their effects on the post-logging rate $\beta ^{Sl}_{j, k}$ in subplot j at time t_k .

Recruits When survivors are all dead, newcomers or recruits will constitute the new forest. We made the assumption that the recruits' annual ACS changes will converge to constant values, with ACS gains compensating ACS losses. Because there are no recruits yet at t₀, recruits' annual ACS growth (Rg) and ACS loss (Rl) are zero, and progressively increase to reach their asymptotic values. Recruits' annual ACS growth and ACS loss can be modeled with the function:

where t the time since the beginning of the recovery period. In the same logic as survivors' cumulative ACS changes, α is the asymptotic value of recruits' annual ACS change (Mg C ·ha^-1·yr^-1), and β is the rate at which this asymptotic value is reached. Contrary to recruits' annual ACS growth and ACS loss, the ACS of new recruits (Rr, the ACS of tree reaching the 10 cm DBH threshold) is high at t₀ because of the competition drop induced by logging, but then progressively decreases to reach its asymptotic value. We modeled it with the following function:

where t is the time since disturbance. As stated before, we chose to model cumulative ACS changes instead of annual ACS changes. The general model for recruits' cumulative ACS changes (ACS growth Rg, ACS loss Rl and ACS of new recruits Rr) is obtained by mathematical integrating from t₀ to t_k annual ACS changes:

where j is the subplot, p the plot, t_k is the time since t₀, R is the annual ACS change, either Rr, Rg or Rl and (σ^R)² the variance of the model. When R is Rg or Rl, η=-1; when R is Rr, η=1. Because of our nested design with subplots j within plots p, we modeled the $\alpha ^{R}_{p}$ values with a random plot effect of mean $\alpha _{0}^{R}$ and variance $\left (\sigma ^{R}_{\alpha }\right)^{2}$:

When the dynamic equilibrium is reached, annual ACS gain (growth and recruitment) compensates annual ACS loss (mortality). We thus added the following constraint for every plot p:

Using the same logic as for survivors, we are interested in predicting variation in β^R as follows:

with R being Rg, Rl or Rr depending on the process we were interested in, with $\beta ^{R}_{0}$ the model intercept, $\lambda ^{R}_{l}$ the vector of l parameters associated to the covariates V_{j, k, l} for which we looked at their effects on the post-logging rate $\beta ^{R}_{j, k}$ in subplot j at time t_k .

Model inference

Bayesian hierarchical models were inferred through MCMC methods using an adaptive form of the Hamiltonian Monte Carlo sampling (Carpenter et al. 2017). Codes were developed using the R language and the Rstan package (Carpenter et al. 2017). A detailed list of priors is provided in Table 2.

Table  2.  List of priors used to infer ACS changes in a Bayesian framework

Model	Parameter	Prior	Justification
Sg	$\alpha ^{Sg}_{p}$	$\mathcal {U}(10, 200)$	Around 100 survivors/ha storing 0.1 to 2.0 MgC each
Sg	$\beta ^{Sg}_{j, t}$	$\mathcal {U}(0, 0.25)$	$12 < {t^{Sg}_{0.95}}^{*} < +\infty$
Sl	$\beta ^{Sl}_{j, t}$	$\mathcal {U}(0, \beta ^{Sg}_{j, t})$	$t^{Sg}_{0.95} < {t^{Sl}_{0.95}}^{*} < +\infty$
Rr	$\alpha ^{Rr}_{p}$	$\mathcal {U}(0.1, 1)$	TmFO observed values (Piponiot et al. 2016b)
Rr	$\beta ^{Rr}_{j, t}$	$\mathcal {U}(0, 0.75)$	$4 < {t^{Rr}_{0.95}}^{*} < +\infty$
Rr	$\alpha ^{Rg}_{p}$	$\mathcal {U}(0.1, 3)$	Amazonian values (Johnson et al. 2016)
Rr	$\beta ^{Rg}_{j, t}$	$\mathcal {U}(0, 0.5)$	$6 < {t^{Rg}_{0.95}}^{*} < +\infty$
Rr	$\beta ^{Rl}_{j, t}$	$\mathcal {U}(0, 0.5)$	$6 < {t^{Rl}_{0.95}}^{*} < +\infty$
All models M∗∗	$\lambda ^{M}_{l}$	$\mathcal {U}(-\beta ^{M}_{j, t}, \beta ^{M}_{j, t})$	avoid multicollinearity problems
Models are : (Sg) survivors’ ACS growth, (Sl) survivors’ ACS loss, (Rr) new recruits’ ACS, (Rg) recruits’ ACS growth, (Rl) recruits’ ACS loss * t0.95 is the time when the ACS change has reached 95% of its asymptotic value **M is one of the five models, either Sg, Sl, Rr, Rg or Rl

 | Show Table

DownLoad: CSV

Identifying the key drivers of the post-disturbance system recovery

To assess the importance of the pre- (forest structure, environment and composition) and post- (climate stress) disturbance forest conditions, we simulated different scenarios modifying the covariate values but keeping an averaged (set to 0) disturbance intensity DIST. Note that all model covariates V_{j, t, l} were standardized before modeling so that, for a given covariate, a - 2, 0 or 2 value respectively refers to a very low, average or very high observed value.

Forest structure The effects of a regenerating (high stem density S_N =1, low basal area S_BA =-1), intermediate (medium stem density S_N =0, medium basal area S_BA =0) and mature (low stem density S_N =-1, high basal area S_BA =1) pre-disturbance forest structure on ecosystem recovery were compared.

Forest environment The effect of three contrasted forest environment were compared: predominance of bottom-lands (high proportion of bottom-lands E_BOTTOM =2, medium slope values E_SLOPE =0, medium altitudinal heterogeneity E_HETE =0), predominance of slopes (medium proportion of bottom-lands E_BOTTOM =0, high slope values E_SLOPE =2, medium altitudinal heterogeneity E_HETE =0) and hilly landscapes (medium proportion of bottom-lands E_BOTTOM =0, medium slope values E_SLOPE =0, high altitudinal heterogeneity E_HETE =2).

Forest composition The effect of pre-logging forest community dominated by conservative tree species (high wood density T_WD =2, high seed mass T_SEED =2, low specific leaf area T_SLA =-2, high maximal stature T_DBH95=2), by a disturbed community (low wood density T_WD =-2, low seed mass T_SEED =-2, high specific leaf area T_SLA =2, low maximal stature T_DBH95=-2) and by a true pioneer community (very low wood density T_WD =-4, very low seed mass T_SEED =-4, very high specific leaf area T_SLA =4, very low maximal stature T_DBH95=-4). The values of the last scenario may appear extreme but note that the model was calibrated with mature forest stands only so that covariate values have to be set out of the calibration range to get a true pioneer community.

Climate stress The effects of a wetter (nor or a few seasonal drought C_DROUGHT =-2, high soil water saturation during the wet season C_WET =2), a drier (seasonal droughts C_DROUGHT =1, medium soil water saturation during the wet season C_WET =0) and a even drier (heavy seasonal droughts C_DROUGHT =2, medium soil water saturation during the wet season C_WET =0) climate on ecosystem recovery were compared.

Sensitivity analysis

To assess the sensitivity of the ecosystem recovery process to the pre- (forest structure, environment and composition) and post- (climate stress) disturbance forest conditions, we simulated the model for an average disturbance intensity DIST=0 and, for each group of covariates V_{j, t, l≠DIST}, varying the values within a group of covariate while setting the other covariates to 0. In a nutshell, for each group of covariates Climate, Composition, Environment and Structure (ⅰ) we independently sampled covariate values from $\mathcal {U}(-2;2)$ while the covariates from the 3 other groups are set to 0, (ⅱ) we ran the model using the sampled covariate values for a set of 100 parameter values drawn from the posterior chains, (ⅲ) we estimated, after 30 years of simulation, the net carbon balance and (ⅳ) we did the procedure 1000 times per group of covariates. Doing so, the variability of the net carbon balance after 30 years reflects the sensitivity of ecosystem recovery to the varying group of covariates.

Results

Given that all the covariates were standardized before modeling, the absolute values of their associated parameters give the weight of each variable in shaping the rates β at which the ACS changes reach their asymptotic state. Negative covariate values indicate slowing and positive values indicate accelerating rates. The values of the disturbance intensity DIST parameters always ranked among the highest absolute values, with negative values for Survivors^′ACSgrowth and Newrecruits^′ACS (Fig. 1a, c) and positive ones for the other three cumulative fluxes (Fig. 1b, d, e).

Figure  1.  Effect of covariates on the rate at which post-disturbance ACS changes converge to a theoretical steady state. Covariates are the disturbance intensity DIST (in green), the structure of the forest before disturbance (the number of individual trees S_N, the basal area S_BA, in pink), the functional trait composition of the forest before disturbance (seed mass T_SEED, specific leaf area T_SLA, wood density T_WD, maximum size T_DBH95, in maroon), the local environment (proportion of bottomlands E_BOTTOM, average slope E_SLOPE, altitudinal heterogeneity E_HETE, in blue) and the climate stress (drought intensity C_DROUGHT, days with water saturation C_WATER, in red). Covariates are centred and standardized. Posterior distribution (median and 0.95 credibility intervals) are reported. Negative covariate values indicate slowing and positive values indicate accelerating rates. a Survivors' ACS growth. b Survivors' ACS loss. c New recruits' ACS. d Recruits' ACS growth. e Recruits' ACS loss

DownLoad: Full-Size Img PowerPoint

Pre-disturbance forest structure

The two variables chosen to describe the pre-disturbance forest structure, i.e. basal area S_BA and stem density S_N, have contrasted behaviors. Basically, associated parameters get their highest absolute values for the two growth models, both Survivors' and Recruits' (Fig. 1a and d), while being close to zero for the other three models (Fig. 1b, c, e). Contribution of Survivors' growth to ACS recovery is higher but declines quicker with high pre-disturbance S_BA values and low S_N values (Fig. 1a). Contribution of Recruits' growth to ACS recovery declines slowly with high S_N values (Fig. 1d). Mature forests (high S_BA, low S_N ) thus recovers faster than regenerating ones (Fig. 2). The relative importance of the pre-disturbance forest structure on the variability of ACS recovery rates is low, with 30 to 50 MgC ·ha^-1 recovered after 30 years (Fig. 6).

Figure  2.  Predicted contribution of annual ACS changes in three contrasted scenarios of pre-logging forest structure : Regenerating, Mature and Intermediate are defined with standardized covariates S_N and S_BA respectively set to [1, -1, 0] and to [-1, 1, 0]. The white line is the net annual ACS recovery, i.e. the sum of all annual ACS changes: survivors' ACS growth Sg and loss Sl, new recruits' ACS Rr and recruits' ACS growth Rg and loss Rl. Dotted lines are out of the calibration period (0-30 year). Maximum-likelihood predictions for ACS stocks (bottom-right) are projected within their credibility intervals (areas with higher levels of transparency

DownLoad: Full-Size Img PowerPoint

Figure  6.  Estimating the relative importance of the Climate stress, the pre-disturbance forest Composition, Environment and Structure in driving ecosystem recovery 30 years after disturbance. The violin plots represent the variability of the distribution of the net carbon balance when covariates within a given group are independently and randomly drawn from $\mathcal {U}(-2;2)$ while the covariates from other groups are set to 0 (the longest the boxplot the highest the sensitivity of ecosystem recovery to this given group of covariates)

DownLoad: Full-Size Img PowerPoint

Local environment

Three variables inform on the forest local environment, i.e. altitudinal heterogeneity E_HETE, proportion of bottomlands E_BOTTOM and the average slope E_SLOPE . The latter never significantly contributes to the β variability (Fig. 1). E_BOTTOM is important in defining recruits ACS fluxes, with positive parameter values for New recruits' ACS and negative ones for growth and loss (Fig. 1c, d, e). All fluxes together, recovery rates do not differ markedly with environmental conditions (Fig. 3) so that the importance of the local forest environment on the variability of ACS recovery rates is quite low, with 30 to 50 MgC ·ha^-1 recovered after 30 years (Fig. 6).

Figure  3.  Predicted contribution of annual ACS changes in three contrasted scenarios of forest environment : Bottomlands, Slopes and Hilly environment are defined with standardized covariates E_BOTTOM, E_SLOPE and E_HETE respectively set to [2, 0, 0], [0, 2, 0] and [0, 0, 2]. The white line is the net annual ACS recovery, i.e. the sum of all annual ACS changes: survivors' ACS growth Sg and loss Sl, new recruits' ACS Rr and recruits' ACS growth Rg and loss Rl. Dotted lines are out of the calibration period (0-30 year). Maximum-likelihood predictions for ACS stocks (bottom-right) are projected within their credibility intervals (areas with higher levels of transparency

DownLoad: Full-Size Img PowerPoint

Pre-disturbance forest composition

Four orthogonal functional traits (i.e. specific leaf area T_SLA, maximum stature T_DBH95, seed mass T_SLEED and wood density T_WD ) have been retained to summarize differences in pre-disturbance forest composition. All these traits have been found to influence post-disturbance ACS recovery rates (Fig. 1). For Survivors, the contribution of growth to ACS recovery is higher but declines quicker with high T_SLA and low T_WD (Fig. 1a) while losses declines slowly with high T_SLA and high T_SEED (Fig. 1b). For Recruits, the contribution of growth to ACS recovery is higher but declines quicker with high T_DBH95 and low T_SLA and T_SEED (Fig. 1d) while losses declines slowly with low T_WD (Fig. 1e). Forests, for which the pre-disturbance composition is dominated by conservative ecological strategies (high wood density, seed mass, maximal stature and low specific leaf area) recovers faster than disturbed forests dominated pioneer species (very low wood density, very low seed mass, very high specific leaf area and very low maximal stature) (Fig. 4). The relative importance of the pre-disturbance forest composition on the variability of ACS recovery rates is high, with 0 to 100 MgC ·ha^-1 recovered after 30 years (Fig. 6).

Figure  4.  Predicted contribution of annual ACS changes in three contrasted scenarios of pre-logging forest composition: Conservatists, Pioneers and Pioneers++ communities are defined with standardized covariates T_WD, T_SEED, T_SLA and T_DBH95 respectively set to [2, -2, -4], [2, -2, -4], [-2, 2, 4] and [2, -2, -4]. Pioneers++ refer to a true pioneer community (very low wood density, very low seed mass, very high specific leaf area and very low maximal stature). The white line is the net annual ACS recovery, i.e. the sum of all annual ACS changes: survivors' ACS growth Sg and loss Sl, new recruits' ACS Rr and recruits' ACS growth Rg and loss Rl. Dotted lines are out of the calibration period (0-30 year). Maximum-likelihood predictions for ACS stocks (bottom-right) are projected within their credibility intervals (areas with higher levels of transparency

DownLoad: Full-Size Img PowerPoint

Post-disturbance climate

Two variables related to climate stress occurring during the recovery process were tested. The intensity of the dry season C_DROUGHT decelerates the decline of Survivors' ACS changes (Fig. 1a, b) while it accelerates the decline of all the recruits' ACS changes (Fig. 1c, d, e). High soil water saturation during the wet season has a similar effect but with very high β values for recruits' ACS loss (Fig. 1e). The forest recovers faster in the driest climate scenarios and this is mainly due to the important Survivors' ACS losses in the wettest scenarios (Fig. 5). The effect of climate stress on the variability of ACS recovery rates ranks 2nd among the 4 groups of covariates, with 10 to 80 MgC ·ha^-1, depending on the post-disturbance climate conditions, recovered after 30 years (Fig. 6).

Figure  5.  Predicted contribution of annual ACS changes in three contrasted climate scenarios: Wet, Dry and Dry++ climates are defined with standardized covariates C_DROUGHT and C_WATER respectively set to [-2, 1, 2] and [2, 0, 0]. Dry++ refer to an extremely-dry climate (very high seasonal drought in the dry season, medium soil water saturation in the rain season). The white line is the net annual ACS recovery, i.e. the sum of all annual ACS changes: survivors' ACS growth Sg and loss Sl, new recruits' ACS Rr and recruits' ACS growth Rg and loss Rl. Dotted lines are out of the calibration period (0-30 year). Maximum-likelihood predictions for ACS stocks (bottom-right) are projected within their credibility intervals (areas with higher levels of transparency

DownLoad: Full-Size Img PowerPoint

Discussion

In this study, we modeled the post-disturbance ACS fluxes in a neotropical forest and found that by testing a few variables that are related to the main endogenous (forest structure and composition) and exogeneous (local environment and climate stress) drivers of ecological community dynamics, we could successfully predict ecosystem trajectories in a wide range of pre- and post-disturbance conditions. Modeling separately the surviving and recruited cohorts was confirmed to be an important methodological choice (Piponiot et al. 2016b), given that the highlighted drivers did not overlap, whether for the growth (Fig. 1a and d) or the loss (Fig. 1b and e) processes. This suggests that our methodological approach, deciphering ecosystem fluxes by demographic processes, could be very useful to predict the long-term trajectories in highly diverse tropical forests for which precise demographic data may be lacking, but aggregative forest dynamic censuses are available from forest inventories. In this study, the disturbance intensity gradient was induced by combining logging to thinning operations (Table 1). Because of its economic value and implications for forest management, selective logging experiments were set up very early on, and the data gathered by these experiments are unique in terms of experiment duration and spatial extent. Despite the particular nature of logging operations (focus on large and commercially-valuable trees even though logging damage concerns all DBH classes), we believe that our study gives clues on the key drivers of ecosystem recovery after large ACS losses induced by other disturbances (e.g. droughts, fire) that are expected to increase in intensity with ongoing global changes (Bonal et al. 2016).

On disturbance intensity

Disturbance intensity DIST remains, by far, the first predictor of the post-disturbance system trajectory (Fig. 1).

Survivors

High disturbance intensities obviously reduce the residual survivors' ACS so that ACS changes from survivors' growth is lower at the beginning but, because of the lower competition between survivors, they tend to live longer and reach the asymptotic state slowly (lower β, Fig. 1a). The positive β for Survivors' ACS loss, meaning that survivors tend to die faster after high levels of disturbance, may look surprising because this goes against the growth result. We believe that a high tree mortality, due to the low survival of damaged trees in highly disturbed systems, in the early post-disturbance years may have resulted in increased β values. However, those losses should rapidly decrease after a decade (Thorpe et al. 2008).

Recruits

High disturbance intensities alleviate competition, and this is probably why recruits' ACS growth is high just after disturbance in the enhanced growth conditions (Herault et al. 2010) and then quickly decrease (high βs, Fig. 1d and e). In these disturbed forests, intense self-thinning (Feldpausch et al. 2007) may explain the fast but limited-in-times ACS losses from survivors' mortality (Fig. 1e).

Endogeneous drivers

Pre-disturbance Forest structure All else being equal, mature forests (high S_BA, low S_N ) recover faster than regenerating ones. And this is mainly due to the higher ACS incoming fluxes from Survivors' growth (Fig. 2). Regenerating forests are composed of short-living fast-growing small species. These species are poorly efficient at carbon accumulation because of their limited growth response to canopy openings and competition alleviation. Indeed, it has been shown in the Paracou forests that species with the highest inherent growth rate (in the absence of disturbance) have the lowest growth response when a disturbance occurs (Herault et al. 2010). On the contrary, large mature trees are, despite their low numbers in many forests, key elements of carbon storage (Lindenmayer et al. 2012) and dynamics (Sist et al. 2014). Previously disturbed, logged or secondary forests, for which forest structure is characterized by low S_BA, and high S_N, may thus be far less resilient to new disturbance than natural undisturbed forests. This also means that post-logging ACS recovery that is currently estimated from the first logging rotation (Rutishauser et al. 2015) may be overestimated for the following logging rotations (Rutishauser et al. 2016). Finally, despite these clear outcomes, the importance of the pre-disturbance forest structure on the variability of ACS recovery rates, as compared to the other covariates (forest composition, environment and climate), remains low (Fig. 6) in the sensitivity analysis. However, we should keep in mind that all the disturbed plots were established in a natural undisturbed forest area so that the model has been parameterized with low pre-disturbance forest structure variability. We thus suggest that, in landscapes with contrasted and tumultuous history, the role of the forest structure in shaping the post-disturbance system trajectory would be much higher.

Pre-disturbance Forest composition The importance of the pre-disturbance forest composition on the variability of the post-disturbance ACS recovery rates is unexpected, with 0 to 100 MgC ·ha^-1, depending on the initial species assemblage, recovered after 30 years (Fig. 6). All the studied functional traits are implied in shaping one or more of the investigated ACS changes. When comparing two typical forest composition, i.e. an assemblage of conservative trees (high wood density, seed mass, maximal stature and low specific leaf area) and an assemblage of pioneer trees (exact opposite trait composition), the pioneer assemblage recovers much slower. This difference is mainly due to highly contrasted survivors' ACS changes. Both ACS growth and loss increase rapidly in surviving pioneers. This result is consistent with the acquisitive strategy of species with a high carbon budget (Sterck et al. 2011), that are well-known to have very fast turn-over rates (Aubry-Kientz et al. 2013; Hérault et al. 2011; Flores et al. 2014). Even if both survivors ACS growth and loss are boosted in pioneer assemblages, it is quite remarkable that the ACS balance in time is mainly under the control of survivors' ACS loss Sl. Survivors' ACS loss can cancel survivors' ACS growth in pioneer-dominated communities, resulting in a null ACS balance (zero is included in the simulations results, see Fig. 6). Why do those pre-disturbance pioneer communities have such high ACS losses in the post-disturbance times? Fast-growing pioneers are generally both poor competitors and poor stress-tolerant trees (He et al. 2013). A first possible explanation is thus that the stress induced by disturbance may be too high for these species that undergo, after disturbance, heavy losses. An alternative explanation is that the higher ACS growth mechanically induces, after a while, higher ACS losses. If pioneers grow faster as a result of growth-stimulating disturbance, they will pass through their natural life span faster, resulting in a transitory gain in carbon storage followed by a massive carbon release when these pioneers get older (Körner 2017). Introducing a time lag for the ACS loss models would be the only way to test the last hypothesis.

Exogeneous drivers

Local environment The forest local environment defined by the altitudinal heterogeneity, the proportion of bottomlands and the average slope, have been found to be of low importance in shaping variability of ACS recovery rates (Fig. 3). This result is quite surprising, given that the local environment is very often referred to as a driver of ecological processes in tropical forests (Grau et al. 2017), from fine pairwise interactions between individual trees (Kraft et al. 2008) to regional variation in community assemblages (Fayad et al. 2016). For instance, in the Paracou forest, the proportion of bottomlands have been found to be of primary importance for forest dynamics: treefall rates are twice as high as on hilltops and tree recruitment and growth rates are higher, leading to a lower basal area and ACS (Ferry et al. 2010). Nearly three fourths of the Paracou taxa are locally distributed as a function of relative elevation, with seasonally inundated bottomlands and well-drained plateaus revealing contrasted species associations (Allié et al. 2015). Despite the relative importance of E_BOTTOM in defining recruits' ACS changes, with positive β parameter values for new recruits' ACS and negative ones for recruits' ACS growth and loss (Fig. 1), all in all ACS recovery rates differ very little from hilly or sloppy plots (Fig. 3). On the one hand, the low-stress conditions of bottomlands (no seasonal drought, less wind) should induce a faster ecosystem recovery. On the other hand, the lower final ACS (Ferry et al. 2010) may mechanically lead to lower absolute carbon storage during recovery. All together, the two processes could be canceling each other, explaining why absolute carbon recovery is similar between bottomlands and hilltops. We also should keep in mind that the disturbance experiment was made by logging. During logging operations, bottomlands are avoided and logging is preferentially conducted in easier-to-access hilltop areas, whatever their proportion in the plot. This may have artificially reduced the environmental difference between logged plots and, in turn, the ACS recovery trajectories.

Post-disturbance climate stress Two seasonal climate stresses were studied: soil water saturation in the wet season and drought intensity in the dry season. The importance of the post-disturbance climate stress on the variability of ACS recovery rates was very high with, depending on the climate scenarios, 10 to 80 MgC ·ha^-1 recovered after 30 years (Fig. 6). The 2 driest scenarios recover initial ACS very quickly, i.e. in less than 60 years while the wettest one would reach a new asymptotic values, far below the initial system ACS (Fig. 6). The main difference between the 3 scenarios lies in the absolute values of the Survivors' ACS loss, with very high values for the wet scenario. This result may look strange given that drought has often been identified as one of the main climate drivers of tropical forest dynamics (Bonal et al. 2016; Wagner et al. 2012, 2013, 2014, 2016), with large mortality events among tropical trees during El Nino years for instance (Phillips et al. 2009) that have not only immediate but also long-term and cumulative impacts on the carbon cycle (Doughty et al. 2015). Those large mortality events are associated with tree hydraulic traits, the most susceptible species being those having a low hydraulic safety margin (Anderegg et al. 2016). In this context, why do the most intense dry season generate the lowest carbon losses? In our training dataset, the natural variability of the total rainfall from 2700 to 3100 mm ·yr^-1 is quite low and, moreover, is far above the 1500 mm ·yr^-1, the evapotranspiration threshold. This means that our experimental forest, located in the Guiana Shield, is not water-limited at all (Stahl et al. 2013), just like the Northern part of the Amazonian basin (Wagner et al. 2017). One may also expect that, because of hydraulic failure, standing death is more frequent during the driest years but, when plotting the tree mode of death registered at Paracou for each dead tree against the drought estimator, no evidence was observed for a potential trend (Aubry-Kientz et al. 2015). Our results thus reinforce the idea that the dominant seasonal climate stress in the Paracou forest is not drought during dry season but water saturation during wet season. This confirms the hypothesis that waterlogged soils in space or in time are risky for trees (Ferry et al. 2010). Moreover, during the rainy season, strong rainfall events often come with strong winds that may reinforce ACS losses (Toledo et al. 2011) and we know, from the Paracou dataset, that the highest total precipitation leads to the highest proportion of tree-fall deaths (Aubry-Kientz et al. 2015). Global climate models converge to simulate, at least for the Amazonian region, a change in precipitation regime over the coming decades (Malhi et al. 2009). Seasonal droughts are expected to become longer and stronger in the future (Joetzjer et al. 2013). Our simulations would suggest that post-disturbance forest recovery would be faster with these new climate conditions. However, we should keep in mind that our simulations were based on a model calibrated with data from a natural, undisturbed forest (Fargeon et al. 2016). With increasing mortality rates due to increasing drought occurrence and severity, the new tree community may be richer in post-disturbance pioneer species. And we have already seen that these new assemblages will slowly recover (Fig. 4) so that recurrent climate-stress in time would not lead to faster recovery rates, but rather to pioneer-rich forest communities with slow recovery rates.

Conclusion

More than half of the tropical forest area are currently designated by National Forests Services as production forests (Blaser et al. 2011) and they consequently play a key role in the tropical forest carbon balance (Piponiot et al. 2016a; Sist et al. 2015). In the Amazon, forest logging and degradation combined to climate change would render up to 80% of the forest area susceptible to major disturbance events in the coming decades (Asner et al. 2010). We have shown that the pre-disturbance forest composition and the post-disturbance climate conditions are of primary importance to predict the recovery potential of tropical forest ecosystems. From the Paracou long-term experiment, it becomes increasingly clear that highly-disturbed forests, because they contain a lot of pioneer species (Baraloto et al. 2012), will be less able to cope with (ⅰ) new disturbance such as logging and (ⅱ) the drier conditions induced by climate change. In other words, already-disturbed forests are likely to be the most vulnerable systems in the current global change context. Forest managers should thus (ⅰ) encourage the development of Reduced-Impact Logging techniques in order to minimize disturbance intensity and (ⅱ) pay a deep attention when drawing management plans to avoid logging pioneer-rich forest units. In the context of increasing disturbances on tropical forests, the lower capacity of disturbed forests to recover is not good news in our fight against climate change.

Acknowledgements

We are in debt with all technicians and colleagues who helped setting up the plots and collecting data over years. Without their precious work, this study would have not been possible and they may be warmly thanked here.

Availability of data and materials

Datasets supporting the conclusions of this article are available upon request to the scientific director of the Paracou Station (https://paracou.cirad.fr).

Authors' contributions

BH and CP drew the concept of the paper, developed the modeling mathematical framework and worked on the parametrization of the model. BH worked on the first draft and both authors improved and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.