Quantifying the effect of forest age in annual net forest carbon balance

We used a global dataset of 126 EC forest sites ranging from 0 to ≥300 year-old stands (table S1 and figure S1 is available online at stacks.iop.org/ERL/13/124018/mmedia). The sites were part of both version 2 of the LaThuile FLUXNET and the FLUXNET 2015 datasets (https://fluxnet.fluxdata.org) of the FLUXNET network (Baldocchi et al 2001, Baldocchi 2008). Five vegetation types were considered in the study, including 76 evergreen forests, 27 deciduous forests, 11 mixed forests, seven shrublands, three savannas and two wetlands.

We aggregated daily NEP, GPP, ER and the associated uncertainties into annual sums (i.e. site-years) and computed an annual-average from all available years per site (i.e. site-average) (see supplementary information for details on EC data processing). One relevant aspect to consider is that the observation-derived GPP is determined via the measured night-time NEP (Reichstein et al 2005). This challenges the statistical independence of both variables, therefore risking a spurious correlation between GPP and NEP at annual scales (see supplementary materials). To avoid any spurious relationship between NEP and GPP, we used a proxy for GPP, i.e. GPP', which was determined as the ratio between latent heat flux (LE) and the square root of vapor pressure deficit $\left({{\rm{GPP}}}^{{\prime} }=\tfrac{{\rm{LE}}}{\sqrt{{\rm{VPD}}}}\right),$ (see supplementary materials for more details and figure S3) as based on physiological principles of the coupling between the C and water cycles at the leaf level (Chen et al 2002, Katul et al 2009, 2009). While LAI and FPAR could have been used as proxies for GPP, GPP' is the integrated response of phenology and physiology, therefore a direct metric of primary productivity, while the former are mostly phenology driven.

We considered forest age as the time since forest establishment or as the time since the occurrence of the last stand-replacing disturbance (see supplementary information for more details on the definition of forest age). Sites were selected based on the availability of information about forest age, disturbance history and management that would allow for an appropriate definition of a meaningful site stand age. These included a range of young and old growth forest sites (figure S2) that were established after complete, nearly complete or substantial removal of forest vegetation (e.g. harvest, fire, wind-throw, insect outbreaks), followed by reforestation/succession/afforestation activities within the flux tower footprint. We did not consider sites with ambiguous historical information or those that had experienced only low to partial intensity disturbances, which would not allow the determination of whole stand forest age. For uneven-aged stands, we followed Spies and Franklin (1991) and estimated the age of a stand as the age of the largest 10% of trees. Undisturbed old-growth forests where age information was available were also included. The information for each site was obtained from the literature, provided by the site principal investigators or from the Biological, Ancillary, Disturbance and Metadata database (table S1). In general, the wide span in stand age among sites and the multi-year record of observations per site permit evaluating the effect of age on both the mean and the interannual variability in NEP.

In addition to the C and water fluxes, we also obtained the following variables as statistical covariates for model development for each EC site:

• (i)  
  local microclimatic variables from in situ observations (i.e. mean annual air temperature (MAT), total annual precipitation (P), and global radiation (Rg));
• (ii)  
  information on nutrient availability (NA) divided into three classes: low NA (n = 67 sites); medium NA (n = 41 sites); and high NA (n = 18 sites) (based on (Fernández-Martínez et al 2014) study and/or expert knowledge, tables S1);
• (iii)  
  additional information on soil texture up to 2 m depth from the SoilGrids1km dataset (Hengl et al 2014);
• (iv)  
  information on forest management based on Campioli et al (2015), Luyssaert et al (2008) datasets, and indications from the PIs (managed forests n = 44 sites; and unmanaged forests n = 81 sites). Managed sites were dominated by human activity while unmanaged sites were undisturbed or experienced low human impacts;
• (v)  
  gridded monthly temperature and precipitation observations from the Climate Research Unit (CRU, http://cru.uea.ac.uk) (Harris et al 2014) to determine long-term linear climate trends and anomalies; and
• (vi)  
  local total atmospheric nitrogen deposition (N_deposition) from in situ observations. We collected N_deposition estimates from the gridded emissions dataset (Wang et al 2017) at 0.5° × 0.5° resolution when they were not available at site level.

The development of an NEP statistical model principally aimed to provide a data-driven estimate of the several factors that control the temporal and spatial variability in NEP, and further to estimate the relative contribution of the different predictive variables—especially age and GPP'—to variation in NEP. To do so, we used the aforementioned statistical covariates (i.e. GPP', microclimatic data, nutrient availability, soil texture, N_deposition, forest age, and forest management) to train and evaluate the ability of a random forest (RF; (Breiman 2001)) algorithm to explain NEP variability.

GPP' and ER are co-determinants of NEP (figure S3), therefore both were initially compared to represent site level NEP. However, given the higher correlation of NEP with GPP' compared to the relationship with ER, and the aforementioned statistical dependence between NEP and the gross C fluxes (i.e. ER and GPP), we chose to discard ER in the statistical analysis.

The role of forest age as an explanatory variable of NEP was additionally evaluated with a published nonlinear model (from now on identified as f(age)) to represent the NEP-age relationship (equation (1) and figure 1)

$$\begin{eqnarray}&&{\rm{NEP}}=e+a\cdot (1-{{\rm{\exp }}}^{b\cdot {Age}})\end{eqnarray} \tag{ 1 }$$

adapted from Amiro et al 2010, where the model parameters a, b and e (offset) were estimated.

Zoom In Zoom Out Reset image size

Although the NEP-age model (equation (1)) was originally developed to represent the temporal patterns of annual GPP-to-ER ratio in forest chronosequences, here we used it to describe the dependency of spatio-temporal NEP variability on forest age. The selection of the Amiro et al (2010) model to describe NEP spatial temporal dynamics assumes that the age-related patterns in GPP/ER are qualitatively similar to NEP, and is supported by a comparison to two other different empirical models presented in the literature (Coursolle et al 2012, Tang et al 2014). These two models were also tested but showed poorer model performance than the Amiro et al (2010) model in the multivariate analysis (table S2).

Estimation of forest age model parameters: The model parameters of the forest-age model (equation (1)) were first estimated in a leave-one-site-out cross-validation (CV) mode for the entire dataset based on a generalized nonlinear least squares (gnls) model using R software (Team 2015). To account for temporal auto-correlation of the observations, we combined the gnls model with an auto-regression moving average model. We minimized the sum of squared residuals weighted for the uncertainty of the observations (Richardson and Hollinger 2005). The standard errors of model parameters were estimated using a bootstrapping algorithm (N = 500 random resamplings). The model output (i.e. f(Age)) was further included as a covariate in the training of the RF algorithms as a nonlinear formulation of age effects on spatio-temporal variability of NEP.

RF algorithm and variables selection: We tested and assessed a RF algorithm using the caret R package (Kuhn et al 2008). This is a non-parametric technique, i.e. it makes no assumption about the residuals of the model. A priori, we used a regression algorithm (i.e. the Boruta algorithm (Kursa and Rudnicki et al 2010)), to determine the best set of predictive variables for NEP among the aforementioned variables. The Boruta method relies on an RF method and determines relevance of each variable by comparing the relevance of the original predictors to that of the randomized variables. It iteratively removes the features that are shown by a statistical test to be less relevant than random probes. We decided to keep the five best variables to improve the accuracy of predictions.

Model performance and model sensitivities: The performance of the statistical model was evaluated by directly comparing model estimates with observed values of NEP for each site-average or site-year in a leave-one-site-out CV mode. In other words, we excluded one site at a time in every training set to predict the mean NEP (site-average) or the annual variations in NEP (site-year) at the excluded site. The statistics used to analyze the results included the coefficient of determination (R²), Nash–Sutcliffe model efficiency coefficient (NSE), root mean squared error (RMSE), and mean absolute error (MAE) (Omlin and Reichert 1999). To quantify the importance of each predictive variable, we performed a model sensitivity analysis by removing a predictor from each regression analysis, then refitted and re-assessed the model without the variable left out. For the site-average analysis, the statistical model was trained using the site-years dataset. The same predictions were further averaged per site and compared to the site-average observations.

The statistical dependence of NEP on forest age supports the nonlinear NEP-age relationship (figure 2 and equation (1)) in that NEP increased rapidly with age followed by stabilization with forests aging (figure 2). This finding reflects expected age-related change in the size and the dynamics of the C pools (Noormets et al 2015, Woodall et al 2015, Zhang et al 2015, Gray et al 2016). However, NEP-age dynamics appeared to depend of climatic conditions, because biomass accumulation rates of regrowing forests vary with climate (Chazdon et al 2016, Anderson-Teixeira et al 2013). This covariation partly explains the differences of the timing when a maximum NEP is reached and then gradually decreases as forest ages among different environmental conditions and the substantial scatter of observations around the model response due to inter-site variability.

Zoom In Zoom Out Reset image size

The low correlation coefficient between NEP and age (table 1) (R² = 0.09 and 0.2 for site-years and site-average, respectively) could be attributed to the substantial contribution of other environmental factors to the spatial and temporal variability of NEP (figure 2). A model based on forest age alone is unable to capture such dependencies and warrants the need to include additional factors in a regression analysis. Although the regression does not show a substantial correlation, the fitted function showed a strong statistical significance, mostly because of the initial curve inflection attributed to the large effect of disturbances on the NEP fluxes in the first years of a recovery process (Amiro et al 2010, Coursolle et al 2012, Tang et al 2014).

Based on the aforementioned feature selection criterion (table S3), the RF algorithm accounted for the effects of forest age (i.e. age and f(Age)), GPP', MAT, and N_deposition. Both site-years and site-average variability were well captured by the different RF models (NSE = 0.62 for site-years and NSE = 0.71 for site-average) (figure 3), suggesting that the structure of the models is suitable for reproducing the spatio-temporal patterns of annual NEP. In addition, for both scenarios (i.e. site-years and site-average), we found that a model including only forest age and GPP' as predictive variables (i.e. NEP = f(age, f(Age), GPP')) had a good predictive capacity for both site-years and site-average (NSE = 0.60 for site-years and NSE = 0.67 for site-average). Although we depicted some high values in the residuals across-site (maximum = 454.4 gC m⁻² y⁻¹; minimum = −537.4 gC m⁻² y⁻¹), we found no significant patterns of residuals against covariates (figure S4).

Zoom In Zoom Out Reset image size

Model sensitivity tests whereby predictors were sequentially removed (table 2) supported the importance of forest age for explaining NEP variability. Whenever we removed forest age from the RF models, model performance decreased, while there were only small changes in model performance when removing either GPP', MAT or N_deposition.

Based on theoretical principles of the C cycle at ecosystem scale, forest age is expected to play a significant role in NEP. Consistent with the early forest dynamics theory on net primary productivity (NPP) trajectories with forest age (Odum 1969), we empirically found strong support for a nonlinear relationship between NEP and forest age, although an age effect was not evident when looking at a univariate relationship, due to spatial variability of other local covariates (figure 2). Hence, we followed a multi-variate approach (figure 3) that accounts for the co-variarying effects of other factors that change in space and time in order to assess the role of age in explaining NEP spatio-temporal variability. Furthermore, RF models have no prior assumption on the functional response between dependent and independent variables, therefore the relevance of forest age (i.e. age + f(age)) was also addressed by contrasting the model performance when removing variables (table 2) and by looking at the model residuals across age class (figure 4). Forest age emerges as the variable that explains most of the spatial and temporal variability in NEP, despite including information on climate and environmental conditions. Photosynthesis and respiration processes drive the link between NEP and forest age, therefore having long term time series of all component fluxes—enabling to establish individual curves per site—and observationally independent estimates of GPP/ER/NEP, could help disentangling whether the NEP-age dynamics are driven by the links between GPP and age or between ER and age.

Zoom In Zoom Out Reset image size

Previously, the effect of stand age on the temporal variability of NEP has been demonstrated via the control of age on NPP using a global dataset (Pregitzer and Euskirchen 2004, Tang et al 2014). Similarly, forest age plays a dominant role in explaining spatial variability in NEP in the East Asian monsoon region (Yu et al 2014, Gao et al 2016). Unlike previous studies, we tested the effect of several drivers on both site-average and site-years, allowing us to evaluate both spatial and interannual NEP variability. Both analyses (i.e. site-average, site-years) showed that forest age was one of the main drivers of NEP variability (table 2). However, some factors are temporally invariant at yearly scales (e.g. soil texture) or do not change over time due to data limitations (e.g. nutrient availability), while others could change (e.g. forest age, GPP', MAT). Therefore, the lack of temporal variability in these factors could reduce their contribution to the NEP in the site-years analysis.

GPP has been suggested as one of the main drivers of NEP spatial variability (Fernández-Martínez et al 2014). We found that excluding GPP' from the RF algorithm decreased the model efficiency by 0.02 and 0.04 for site-years and site-average, respectively (table 2). GPP' emerged as superior predictor compared to climate and soil properties (table S3), likely because it is more closely coupled with NEP, whereas climate and soil properties had variable effects depending on site characteristics. The statistical relation of GPP' to NEP appeared to be significantly stronger for stands younger than 20 years, than for intermediate-aged/old-growth forests (≥20 years) (figure S5). In the initial successional pathway, most of the year-to-year variability in NEP is explained by the changes in GPP and climate. However, as forest ecosystems mature and the autotrophic and respiratory processes start to balance each other out, the variations in NEP become more a function of forest age, or time since disturbance, rather than of individual variations in GPP or ER. Having the full representation of stand development stages is important for representing forest spatio-temporal C dynamics after stand-replacing disturbances more realistically (figure 4). This means that the controls of GPP on NEP are strongly dependent on the distribution of forest age, which emphasizes the relevance of age class distributions for understanding the dynamics of biosphere-atmosphere fluxes. The interactions between forest age and local conditions (e.g. GPP) suggest to move beyond stand age in reflecting changes in plant and soil pools, but also in appropriately parameterizing forest age-related changes in ecophysiological mechanisms both at plant and soil levels. Still, we have limited knowledge on the disturbance effects on detrital pools (and thus heterotrophic respiration, Rh), the type of transition between previous land cover/use (Carvalhais et al 2010) followed by different regeneration types (e.g. regrowth, plantation on pasture, former agricultural lands, and afforestation), and site history. These ecosystem conversions may strongly influence ecosystem C balance (Kutsch and Kolari 2015) and could explain the current bias present for the young forests (<20 years) (figure 4(a)).

While several environmental factors exert controls on NEP, we found that their statistical effect was minor in comparison to forest age. In many cases, sensitivities and even the sign of the relationship between environmental factors and NEP differ among case studies in the literature. For instance, NEP can be positively correlated with MAT in space (Fernández-Martínez et al 2014), whereas other studies find only very weak relationships (Law and Falge 2002, Reichstein et al 2007, Piao et al 2008). In boreal regions, air and soil temperature are the main factors affecting interannual NEP variability in old stands, while climatic conditions could not explain temporal patterns of NEP of young stands (Coursolle et al 2012). Here, MAT had a modest contribution to explaining NEP variability in the final model (table 2).

Rather than mean temperature, temperature changes in the recent past significantly influence current spatial variability of the forest C sink (Piao et al 2009). We tested both annual climate anomalies and climate trends (i.e. from 1960 to 2012 based on the CRU dataset) in the final models, but found that they had limited effect on explaining NEP variability (table S3). Nevertheless, future increases in temperature, changes in precipitation patterns and more extreme events will likely have significant effects on the C budgets of forest ecosystems (Thuiller et al 2011, Trumbore et al 2015).

Soil characteristics and fertility may play an important role in the spatial and interannual variability of NEP (Oren et al 2001, Janssens et al 2010, Fernández-Martínez et al 2014). Ecosystem C exchanges are generally limited by nutrient availability (often N) that may increase following disturbance (e.g. stand replacement, harvest). N mineralization increases available N while N uptake decreases N availability (Thornton et al 2002). However, we found that nutrient availability and clay content were not considered statistically strong drivers in explaining NEP variability. This was in contrast to previous studies concluding that nutrient availability is even more important than forest age in explaining across-site forest NEP variability (Fernández-Martínez et al 2014). This apparent contradiction emerged from the fact that earlier studies used a linear relationship, while we also included a nonlinear relationship. We showed that removing a nonlinear relationship between NEP and stand age (i.e. f(age)) results in a substantial loss in overall model performance and a significant reduction in the apparent importance of forest age for NEP (table S4).

Nevertheless, the apparent low contribution of climate, nutrient status, and soil properties to explaining NEP variability can be explained by the fact that their information is already embedded either in forest age or GPP'. The latter is clearly climatically driven, whereas the GPP-NEP dynamic is strongly controlled by forest nutrient availability (Fernández-Martínez et al 2014). Forest age is likely a superior predictor of spatio-temporal variation in NEP because it integrates relevant ecological information not captured by other single variables. In fact, forest age is rather a composite measure of numerous drivers that are more directly mechanistically coupled with C cycling processes.