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

Forests dominate carbon (C) exchanges between the terrestrial biosphere and the atmosphere on land. In the long term, the net carbon flux between forests and the atmosphere has been significantly impacted by changes in forest cover area and structure due to ecological disturbances and management activities. Current empirical approaches for estimating net ecosystem productivity (NEP) rarely consider forest age as a predictor, which represents variation in physiological processes that can respond differently to environmental drivers, and regrowth following disturbance. Here, we conduct an observational synthesis to empirically determine to what extent climate, soil properties, nitrogen deposition, forest age and management influence the spatial and interannual variability of forest NEP across 126 forest eddy-covariance flux sites worldwide. The empirical models explained up to 62% and 71% of spatio-temporal and across-site variability of annual NEP, respectively. An investigation of model structures revealed that forest age was a dominant factor of NEP spatio-temporal variability in both space and time at the global scale as compared to abiotic factors, such as nutrient availability, soil characteristics and climate. These findings emphasize the importance of forest age in quantifying spatio-temporal variation in NEP using empirical approaches.


Introduction
Forests cover about 30% of the Earth's terrestrial surface and store around 90% of terrestrial vegetation carbon (C) (Canadell et al 2000, Gower 2003, Le Quéré et al 2018, indicating their fundamental role in terrestrial C dynamics (Bonan 2008, Beer et al 2010, Pan et al 2011, Hicke et al 2012, Carvalhais et al 2014. However, the functioning of forest ecosystems is likely to be altered by changing climate (Ciais et al 2005, Xiao et al 2009, Zhao and Running 2010, Reichstein et al 2013, ecological disturbances (Chambers et al 2007, Bowman et al 2009, Amiro et al 2010 and management (Noormets et al 2015, Naudts et al 2016. Therefore, it is important to characterize current and future forest net ecosystem production (NEP) for regional to country-level assessments, and to evaluate mitigation strategies that minimize carbon dioxide (CO 2 ) emissions to the atmosphere (Becknell et al 2015, Trumbore et al 2015, Law et al 2018.
The overall NEP dynamic at a given site emerges from combined responses to factors that control both gross primary productivity (GPP) and ecosystem respiration (ER) (NEP=GPP-ER). At the ecosystem level, the forest NEP patterns following stand-replacing disturbance are mostly controlled by the timevarying dominance between autotrophic and catabolic processes. After disturbance, heterotrophic respiration (Rh) generally tends to increase because of an aboveground biomass transfer to the litter and soil organic matter C pools (Law et al 2003, Kurz et al 2008, Harmon et al 2011, Noormets et al 2012, Lindauer et al 2014, Paul-Limoges et al 2015, while GPP collapses due to an instantaneous reduction in leaf area, resulting in a net release of CO 2 to the atmosphere. On the one hand, a shift from C source to C sink occurs as canopy development supports GPP and net C accumulation in plants increase. On the other hand, Rh and ER decline due to a reduction in litterfall and substrate availability through decomposition. The resulting imbalance between GPP and ER persists until vegetation and soil C pools increase up to the point when ER comes into equilibrium with GPP (Schwarz et al 2004, Lindroth et al 2008, Luyssaert et al 2008, Tang et al 2014.
Several approaches used for assessing forest NEP include micro-meteorological and biometric techniques, process-based models, and/or satellite data. However, annual regional C stock assessments that account for age-related physiology, regrowth, and soil processes following disturbance are challenging due to lack of information in disturbance history or management practices (Zscheischler et al 2017). Flux sites provide the annual net uptake of CO 2 from the atmosphere (i.e. NEP) that can be used to calibrate empirical models for mapping annual NEP at regional scales. However, current empirical uspcaling exercises (Jung et al 2011, Tramontana et al 2016 do not directly include proxies that allow the dynamics of C fluxes with age to emerge, therefore it is not clear how well the aforementioned data-driven models captured such dynamics. Thereby, empirical estimation of annual NEP that explicitly accounts for disturbance and forest age effects are of relevance for regional C stock studies .
Despite the recognized effects of forest age in controlling spatial and interannual variability of NEP, there is still debate about the quantitative role of forest age in the empirical annual forest C estimates. In fact, the most recent observation-based synthesis studies tackling NEP spatio-temporal variability and its drivers reached diverging conclusions on the importance of forest age. While some authors have shown that forest age is a key factor controlling forest C balance (Chen et al 2002, Coursolle et al 2012, Yu et al 2014, Gao et al 2016, others have indicated that spatial and interannual variability of NEP is mainly controlled by nutrients availability and soil properties (Bhatti et al 2002, Janssens et al 2010, Vicca et al 2012, Fernández-Martínez et al 2014 or climate conditions (Thornton et al 2002, Amiro et al 2006, Coursolle et al 2006, although several authors report that the C budget in forest ecosystems is less sensitive to climatic conditions than expected in certain regions (Law and Falge 2002, Reichstein et al 2007, Yi et al 2010. Given the fundamental understanding of the role of forest age in NEP and the contrasting results from previous meta-analyses, we revisited the importance of forest age to the spatial and temporal variability in NEP based on a more up-to-date, larger, and higher quality eddy-covariance (EC) dataset including 126 forest ecosystem sites. We further expanded previous observation-based syntheses by exploring nonlinear empirical model formulations to incorporate forest developmental stage and environmental factors for calculating realistic NEP spatio-temporal variability. Such a model can eventually be used to estimate NEP at a global scale and infer likely limits to NEP variation and the future forest C sink as forests age.

Datasets
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 observationderived GPP is determined via the measured nighttime 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  figure S3) as based on physiological principles of the coupling between the C and water cycles at the leaf level (Chen et al 2002, 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°r esolution when they were not available at site level.

NEP model development
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 variablesespecially 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)  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 spatiotemporal 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).

Experiment design
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 autocorrelation 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 2 ), 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 siteyears dataset. The same predictions were further averaged per site and compared to the site-average observations.

Age-dependent forest carbon dynamics
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 . 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.
The low correlation coefficient between NEP and age (table 1) (R 2 =0.09 and 0.2 for site-years and siteaverage, 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

Random forest algorithm performance and model sensitivity
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 siteaverage) (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 siteaverage). Although we depicted some high values in the residuals across-site (maximum= 454.4gCm −2 y −1 ; minimum=−537.4gC m −2 y −1 ), we found no significant patterns of residuals against covariates ( figure S4). 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 .
4. Discussion 4.1. Forest age as a key driver of spatial and interannual variability in NEP 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  3(b)) site-average. The scatterplots of the observed versus modelled annual sums of NEP were grouped by age classes. NEP=f(age, f(Age), GPP′, MAT, N deposition ). Table 2. Changes in model performance caused by removing predictors from the best model set-up and then refitting the model without the left out variable(s). These results were computed by leave-one-site-out cross-validation. The (-) symbol means that the predictive variable(s) were removed from the Random Forest models. R 2 =coefficient of determination; NSE=Nash-Sutcliffe model efficiency coefficient; RMSE=root mean squared error; MAE = mean absolute error; total n=716 for site-years and n=126 for site-average. NEP=f(age, f(Age), GPP′, MAT, N deposition ). 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 fluxesenabling 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. 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 andEuskirchen 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. siteaverage, 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 siteyears analysis.

Site-years
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)).

Climate and soil properties controls on spatial annual NEP variability
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.

Implications
While GPP, climate and soil properties are significant factors influencing the variability of NEP across space and time, we conclusively demonstrate that forest age performs as a strong indicator of spatio-temporal variability in NEP and is a useful integrated proxy for ecological changes that constrain NEP at the global scale.
Many global ecosystem models rely on simple representations of forest age dynamics and few consider the role of successional changes in C cycling processes (Anderson-Teixeira et al 2013), which introduces uncertainties into long term simulations of forest C dynamics (Friend et al 2014, Friedlingstein et al 2014. Additionally, given the statistical power of the proposed model in comparison to other state-of-the art approaches (Jung et al 2011, Tramontana et al 2016, this study points out new directions towards further developments in bottom-up upscaling exercises based on EC data. Regardless of the modeling strategy, reliable annual maps of forest age distribution and/or disturbance history will be required in order to make accurate predictions of NEP in space and time. These will further support the design of sustainable forest management and climate change mitigation strategies that depend on the effect of forest aging and age class distribution (Pan et al 2011, Thuiller et al 2011, Trumbore et al 2015. This study emphasizes the need for increased focus on forest demography, which may amplify or exceed the importance of climate sensitivity for predicting the future of the terrestrial C cycle. USCCC. We acknowledge the financial support to the eddy covariance data harmonization provided by CarboEuropeIP, FAOGTOSTCO, iLEAPS, Max Planck Institute for Biogeochemistry, National Science Foundation, University of Tuscia, Université Laval, Environment Canada, US Department of Energy and the database development and technical support from Berkeley Water Center, Lawrence Berkeley National Laboratory, Microsoft Research eScience, Oak Ridge National Laboratory, University of California, Berkeley and University of Virginia. This research was funded by both the BACI project and the Independent Monitoring Project. We would like to thank Sönke Zaehle for providing the N deposition data. We thank Mirco Migliavacca and Ulrich Weber for the assistance with the processing of the EC and the CRU data and Andrew Durso for proofreading the manuscript. Finally, we would like to thank Tomislav Hengl for providing the SoilGrids1km dataset and Sara Vicca for helping us in determining the nutrient availability classification of some sites used in this study. This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic within the National Sustainability Programme I (NPU I), grant number LO1415. This research was funded by both the BACI project and the Independent Monitoring Project. We would also like to acknowledge NOVA grant UID/AMB/04085/2013 and the GlobBiomass project.