Ozone and haze pollution weakens net primary productivity in China

Atmospheric pollutants have both beneficial and detrimental effects on carbon uptake by land ecosystems. Surface ozone (O3) damages leaf photosynthesis by oxidizing plant cells, while aerosols promote carbon uptake by increasing diffuse radiation and exert additional influences through concomitant perturbations to meteorology and hydrology. China is currently the world’s largest emitter of both carbon dioxide and short-lived air pollutants. The land ecosystems of China are estimated to provide a carbon sink, but it remains unclear whether air pollution acts to inhibit or promote carbon uptake. Here, we employ Earth system modeling and multiple measurement datasets to assess the separate and combined effects of anthropogenic O3 and aerosol pollution on net primary productivity (NPP) in China. In the present day, O3 reduces annual NPP by 0.6 Pg C (14 %) with a range from 0.4 Pg C (low O3 sensitivity) to 0.8 Pg C (high O3 sensitivity). In contrast, aerosol direct effects increase NPP by 0.2 Pg C (5 %) through the combination of diffuse radiation fertilization, reduced canopy temperatures, and reduced evaporation leading to higher soil moisture. Consequently, the net effects of O3 and aerosols decrease NPP by 0.4 Pg C (9 %) with a range from 0.2 Pg C (low O3 sensitivity) to 0.6 Pg C (high O3 sensitivity). However, precipitation inhibition from combined aerosol direct and indirect effects reduces annual NPP by 0.2 Pg C (4 %), leading to a net air pollution suppression of 0.8 Pg C (16 %) with a range from 0.6 Pg C (low O3 sensitivity) to 1.0 Pg C (high O3 sensitivity). Our results reveal strong dampening effects of air pollution on the land carbon uptake in China today. Following the current legislation emission scenario, this suppression will be further increased by the year 2030, mainly due to a continuing increase in surface O3. However, the maximum technically feasible reduction scenario could drastically relieve the current level of NPP damage by 70 % in 2030, offering protection of this critical ecosystem service and the mitigation of long-term global warming. Published by Copernicus Publications on behalf of the European Geosciences Union. 6074 X. Yue et al.: Ozone and haze pollution weakens net primary productivity


Introduction
The large carbon flux between forest ecosystems and the atmosphere plays an essential role in the global carbon cycle [Beer et al., 2010;Xiao et al., 2011;Xue et al., 2011Xue et al., , 2016. This net carbon flux can be both determined by the rate of carbon assimilation (gross primary production (GPP)) and by the time of carbon being locked up in living plant tissue, i.e., woody residence time (τ w , years) before being released into the atmosphere [Lloyd and Farquhar, 1996;Malhi, 2012]. The τ w is an important physiological parameter that expresses the balance between forest recruitment/growth and mortality when it is in a steady condition . Many previous studies have found a close relationship between τ w and forest carbon stocks and a large amount of uncertainty in the estimation of τ w per se [Vieira et al., 2005;Zhou and Luo, 2008;Zhang et al., 2010;Zhou et al., 2012]. Therefore, quantifying τ w is important to gain a clearer understanding on the global forest carbon budget balance.
Dynamic global vegetation models (DGVMs) are useful tools for mapping global carbon stock and predicting its variation future. In most DGVMs, the global terrestrial ecosystems are represented by different plant functional types (PFTs), and much attention has been paid to improve the modeling of photosynthesis and GPP for each PFT. Typically, carbon stock or biomass has been simulated as the vegetation carbon pool including leaf, stem (including branch), and root produced from the allocation of net primary production (NPP). For forest ecosystems, the total aboveground biomass (AGB, Mg ha À1 ) is mainly determined by the stem biomass. To predict forest growth and biomass accumulation, most DGVMs use the τ w parameter to represent the time carbon remained in an ecosystem. The default PFT-specific value of τ w introduces little spatial variance and is assumed to be temporally invariant for model simulation.
Although DGVMs could simulate biomass at regional to global scales, an increasing number of recent studies showed that the assumed invariant τ w can induce large uncertainties in modeled biomass [Delbart et al., 2010;Castanho et   vegetation to the projected climate by 2100 using seven global vegetation models (including DGVMs) and found large variations in global carbon stocks among these models. This large variation cannot be solely explained by the differences in simulated NPP; when uncertainties in τ w were taken into account, an additional 30% of the variation was explained. In fact, trend of projected NPP were more consistent, while the temporal trend signs of τ w may even be different (negative or positive) [Friend et al., 2014]. These results demonstrated the necessity of accurately estimating τ w during model simulations. A recent study based on forest field observations also showed that the spatial patterns of biomass of lowland Amazonian forests can be better explained by the patterns of τ w rather than by the variation in forest productivity [Malhi et al., 2015]. These authors stated that more attention should be paid to the allocation of NPP and τ w aside from that for GPP [Malhi et al., 2015]. However, the default values of τ w for most current DGVMs are often not representative and may even be highly biased when compared with observed values .
When forest is near equilibrium, biomass production can be balanced by biomass loss, and therefore, τ w could be determined by forest mortality. At a regional or global scale, forest mortality may originate from either intrinsic self-thinning or natural and human disturbances, which makes it difficult to estimate τ w [Stephenson et al., 2011]. Previous studies have shown that τ w could be closely related to abiotic (e.g., temperature and precipitation) or biotic (e.g., stand density) factors [Lines et al., 2010;Galbraith et al., 2013;Malhi et al., 2015]. Lines et al. [2010] compiled forest mortality rates from forest inventory data and explored the influence of environmental and physiological variables for forests of the eastern United States. That study emphasized the influence of climate, soil, species, and size (stem diameter) on forest mortality rates at individual tree level. Thurner et al. [2016] estimated temperate and boreal forest residence time (including both aboveground and belowground forest carbon stocks) at the regional scale using remote sensing data. They found that the forest residence time was closely related to different climate drivers, depending on specific ecosystems. At the global scale, another term-carbon residence time (including both aboveground and belowground carbon stocks)-is used for all terrestrial ecosystems; researchers found that carbon residence time is highly correlated with temperature and precipitation [Carvalhais et al., 2014]. More confidence in projected spatiotemporal variation in forest biomass would be achieved if more accurate τ w or carbon residence time can be obtained by field observations or processes/empirical-based models.
The objective of this study is to investigate the determinants of forest τ w based on field collected values and its role in the improvement of DGVM simulations. Specifically, we initially collected field-derived τ w and analyzed the possible drivers of its change including environmental and biotic variables. Since global natural and anthropogenic disturbances are missing, we thus mainly focus on forest plots without disturbances for a long time. The collected field-derived τ w covers all forest biomes including boreal, temperate, and tropical forests. Second, we generated a global, gridded τ w map using the field-derived τ w values and other meteorological and physiological variables as predictors with the random forest (RF) method. Third, we used the resulting global τ w map to parameterize a DGVM-Integrated Biosphere Simulator (IBIS) [Foley et al., 1996;Kucharik et al., 2000], and based on an overview of different model descriptions of τ w , we assessed how field-derived τ w data can improve biomass simulations. Our results can improve our understanding of global patterns of τ w and thus the global carbon cycle under climate change.

Collection of Woody Residence Time From Field Data
Assuming that forest is near an equilibrium state, τ w can be calculated as follows : where M w is the mean AGB (Mg ha À1 ) and W p is the mean aboveground woody productivity (including stem and branch, Mg ha À1 yr À1 ). Since we mainly focus on aboveground τ w in the present study, AGB and aboveground W p are used, which are different from Thurner et al. [2016]. For forest ecosystems, τ w is closely related to, but not identical with, tree lifetime and is mainly determined by the lifetime of medium and large trees Malhi et al., 2015].

10.1002/2016GB005557
We compiled τ w values for a total of 1319 forest sites from published literature. These sites cover all climate zones from boreal to temperate and tropical forests ( Figure 1). To ensure the accuracy of our collected τ w , we only collected data from sites with no major human or natural disturbances for at least 100 years. This was done to assure to a large extent that forest plots retained were near the equilibrium state. Similarly, only "mature" or "old growth" forests were selected in our analysis; that is, secondary forest sites that have not reached a state of equilibrium were excluded. For most of these sites, τ w is not directly available, and we calculated τ w from the reported AGB and woody productivity (including stem and branch) using equation (1).
Many of our τ w values were obtained from four articles: Galbraith et al. [2013], Luyssaert et al. [2007], Keeling and Phillips [2007], and Guo and Ren [2014]. Most tropical forest sites are from Galbraith et al. [2013], who compiled the τ w values for tropical forests from published literature. The τ w values were calculated either by equation (1) or as the reciprocal of woody turnover time (see the supporting information of Galbraith et al. [2013] for details). Luyssaert et al. [2007] compiled aboveground biomass and NPP data for leaf, stem, and twig data for boreal, temperate, and tropical forests. We selected those forest sites that have historically lacked extensive human management for our analysis. Keeling and Phillips [2007] analyzed the relationship between forest productivity and aboveground biomass by compiling data from published literature. We calculated the τ w values using their data set provided in the supporting information. Guo and Ren [2014] collected plot-level forest productivity and biomass data of natural and planted forests in China from published literature with detailed information on the forest carbon pools and other physiological characteristics of forests such as forest age, and we calculated the τ w values for natural forests from the biomass and productivity data using equation (1). For each site, we also obtained meteorological data, i.e., annual precipitation (P, mm) and air temperature (T,°C) if these values were available in corresponding paper. For those sites where these variables were missing, we retrieved the values from global values of WorldClim (http://www.worldclim.org/).

Analysis of the Drivers of Woody Residence Time
We analyzed the global patterns of τ w according to P and T. Further exploration of the relationships between τ w and biotic variables of carbon use efficiency (CUE, ratio of NPP to GPP)/forest age was also conducted since these variables were considered to be closely related with W p and AGB. Since not all of the plots contain information on forest age (938 plots) and GPP (32 plots) from in situ observations, fewer plots were used during this analysis. To examine how forest τ w is determined by climate and forest age, we used linear mixed effect models (LMMs). The absolute values of soil nutrient contents were not available, and we thus assumed different random effects for local nutrient classification from Harmonized World Soil Database (http://webarchive. iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/). We tested the effects of P, T, and forest age on τ w by running all possible LMMs, which were evaluated by Akaike information criterion (AIC). Models with a ΔAIC ≤ 7 were selected for model averaging and coefficient estimation and then used to predict τ w in relation Global Biogeochemical Cycles 10.1002/2016GB005557 to forest age, P, and T. Forest age and τ w was log transformed as the nonlinear relationship between AGB and age. All explanatory variables were standardized by subtracting the mean from each value and dividing by the standard deviation to allow for easier interpretation of coefficients and improvement of model convergence [e.g., Martin et al., 2016]. All LMM analyses are conducted by R with the nlme [Pinheiro et al., 2015] and MuMIn [Barton, 2015] packages.

Estimation of Global Forest Woody Residence Time
In this study, to put τ w as model parameter for a DGVM, i.e., Integrated Biosphere Simulator (IBIS), we extrapolated the field collected τ w into spatial continuous layers by RF method [Breiman, 2001]. The RF method is a formalized nonparametric machine-learning algorithm and has been successfully used for parameter extrapolation [e.g., Simard et al., 2011;Su et al., 2016]. The main advantage of this method is that it does not require assumptions to be made regarding to the normality of covariables and can minimize the within-group variance. By the inherent unique tree "bagging" algorithm, the RF method can select a random subset of covariables at each candidate split and thus overcome the overfitting habit of decision tree algorithms [Breiman, 2001]. The RF extrapolation method was implemented based on the "randomForest" R package [Liaw and Wiener, 2002], which includes both classification and regression functions. We generated a regression "random forest" by the field collected τ w and possible predictors. Two sets of global forest τ w were generated: (1) present-day spatially variant τ w and (2) projected spatially and temporally variant τ w during 2006-2100 at the annual scale. Since we do not have the global wall-to-wall forest age and independent in situ woody NPP (and thus CUE), the predictors used in the RF method were different from those in the LMM analysis. The present-day τ w was calculated by five ancillary predictors: T, P, GPP, evapotranspiration (ET), and digital elevation model (DEM). In contrast, three predictors, T, P, and DEM, were included to build the RF model and project annual τ w by 2100 . The global GPP/NPP and ET data sets were both available from the Numerical Terradynamic Simulation Group website (http://www.ntsg.umt.edu/biblio). The MOD17 GPP/NPP data sets were at a 1 km resolution on monthly (MOD17A2) and annual (MOD17A3) scales and are available from 2000 to the present year. The Moderate Resolution Imaging Spectroradiometer (MODIS) ET data set (MOD16) was produced by Mu et al. [2011], using a calculation algorithm based on the Penman-Monteith equation [Monteith, 1965] at a resolution of 1 km on a global scale. In our case, the improved version of annual MODIS ET products (MOD16A3) was used in the calculation, which was from a consistent version of meteorological variables from the Global Modeling and Assimilation Office, and this product has generally been validated and has acceptable accuracy [Mu et al., 2011].The DEM data used were obtained from the NASA Shuttle Radar Topography Mission with a resolution of 1 km (http://srtm.csi.cgiar.org/). The projected spatially and temporally variant τ w was calculated annually from the meteorological data from Community Climate System Model version 4.0 (CCSM4.0) simulated climatic variables under the Representative Concentration Pathway 4.5 (RCP4.5) scenario (http://www.cesm.ucar.edu/models/ccsm4.0/).

Model Description of Woody Residence Time
Most DGVMs simulate three carbon pools: leaves, stems (for trees), and roots [e.g., Foley et al., 1996;Sitch et al., 2003]. For a given PFT, an invariant value is commonly used for carbon residence time during model simulation. In IBIS, NPP is allocated among the three carbon pools at an annual scale. In steady state conditions, the instantaneous changes in the biomass pool j of PFT i are represented as where a i, j is the fraction of annual NPP allocated to the biomass pool and τ i,j is the carbon residence time of that biomass pool, which is identical with τ w for the carbon pool of stems and branches. Note that a i, j is a fixed value in IBIS; however, some other DGVMs (such as Lund-Potsdam-Jena [Sitch et al., 2003]) estimate NPP allocation using allometric equations.

Model Simulation of AGB Using Estimated Woody Residence Time
We investigated how τ w would improve the simulation of global forest AGB using IBIS. We first examined the influence of τ w on simulated AGB at the site level using τ w values ranging between 25th and 75th percentiles of each field-based PFT data. We then resampled the two global forest τ w maps to a 0.5°× 0.5°scale to specify Global Biogeochemical Cycles 10.1002/2016GB005557 (1) parameterization with the temporally (and also spatially) dynamic τ w (projected annual τ w , only for 2006-2100), (2) parameterization with the spatially dynamic but temporally invariant τ w (present-day τ w ), and (3) parameterization with the IBIS model default τ w .
The climatic data to drive IBIS, including monthly mean air temperature, precipitation, relative humidity, cloudiness, diurnal temperature range, wind speed, and the number of wet days during 1948-2005 were obtained from the Climate Research Unit [Harris et al., 2014]. These data were used as DGVM inputs to reproduce the amounts of global historical forest biomass during the years 1948-2005. Independent field collected forest AGB data (a total of 716 plots ) with specific recorded measurement time were also compiled to validate the model simulations with improved and default τ w ( Figure S1 in the supporting information). To project the AGB to the end of 2100, the CCSM4.0 simulated climatic variables at the land surface under the RCP4.5 scenario (1.25°× 1.0°) were used to drive the DGVM. The CCSM4.0 data covering the period of 2006-2100 were downloaded and interpolated to a 0.5°× 0.5°resolution for model inputs.

Drivers of Woody Residence Time
The model-averaged coefficients indicate positive relationships between τ w and the logarithm of forest age (slope = 0.34 ± 0.01, importance value = 1, and p < 0.001) and mean annual temperature (slope = 0.07 ± 0.02, importance value = 1, and p < 0.001), which show negative relationship between mean annual precipitation (slope = À0.09 ± 0.02, importance value = 1, and p < 0.001). Furthermore, the slope related to forest age also increased (interaction term = À0.07 ± 0.01, importance value = 1, and p < 0.001) and decreased (interaction term = À0.1 ± 0.01, importance value = 1, and p < 0.001) with mean annual temperature and precipitation. The interaction term between temperature and precipitation is considered to be less important in characterizing forest τ w (interaction term = 0.01 ± 0.01, importance value = 0.32, and p = 0.48). The first two models included in the model-averaging process have great prediction power (conditional R 2 = 0.55; Table 1 and Figure 2), and this results in an improved prediction performance compared with models containing only age, precipitation, or temperature (Table 1).

Estimation of Global Forest Woody Residence Time
The percentage increase in mean squared error (%IncMSE) and node purity (IncNodePurity) was calculated to evaluate the importance of each predictor (Figure 3). The absence of the annual precipitation and GPP will significantly increase mean squared error and node purity and thus lower the predictive ability of the RF Figure 3. The mean importance of variables for 100 runs woody residency time random forest model, denoted by (a) the percentage increase of mean squared error (%IncMSE) and (b) the increase in node purity (IncNodePurity). Temp, Prep, GPP, ET, and DEM represent annual temperature, annual precipitation, gross primary production, evapotranspiration, and digital elevation model, respectively. %IncMSE means the percentage increase in mean squared error when the variable is randomized. IncNodePurity is a specific parameter for regression trees, which is measured by the residual sum of squares. The larger the %IncMSE and IncNodePurity of a variable are, the worse the model performs when the variable is randomized.
Global Biogeochemical Cycles 10.1002/2016GB005557 model. Therefore, these predictors are the two most important ones for the estimation of global τ w . Similarly, ET, temperature, and DEM were also found to be effective predictors for the model.
The τ w exhibits large spatial heterogeneity, with an average value of 66.7 years on a global scale (Figure 4). Large values of τ w were found in tropical areas, especially in the Amazonian forests. The τ w values in this area were generally above 50 years, and larger values above 70 years were observed in the eastern part of Amazonian forest, showing a west-east increasing gradient as highlighted by other authors . In contrast, τ w for central African forests has a moderate value of about 50 years. The estimated τ w for boreal forests in central Siberia also had a high value of greater than 100 years; meanwhile, the τ w for Figure 4. Spatial pattern of (a) the estimated woody residence time (τ w , years) based on collected field data and (b) relative difference of estimated τ w . The relative difference is estimated as the ratio of standard deviation of τ w to the resulting τ w . The standard deviation of τ w is calculated when 75% of the collected field plot data were randomly selected in each of 100 random forest simulation model runs. The estimated τ w has a 1 km resolution and was resampled to 0.5°× 0.5°when used as parameters for dynamic global vegetation model simulations in this study.
Global  Figure 5 shows the model performance of the simulated present-day AGB with model "default" and our estimated present-day τ w . Each point in Figure 5 indicates one or more validation plots that are located in the same IBIS grid (0.5°× 0.5°). The simulated AGB using default τ w showed large scattering, with underestimation for large amounts of AGB. Moreover, the default τ w resulted in many small amounts of AGB (close to 0 Mg ha À1 ), which indicates a nonforest PFT for the model simulation. The resulting AGB from the improved τ w has a relatively close relationship with plot values, even though overestimation and underestimation were observed for small and large AGB values (Figure 5b). Note that the simulated AGB using our estimated present-day τ w was also subjected to scattering, which may be caused by the scale difference between measured and simulated AGB (0.01°and 0.5°).

Simulation of AGB Based on Estimated Woody Residence Time
Ten Fluxnet sites, representing five different woody PFTs, were randomly selected to test the AGB uncertainties due to τ w forward in time under RCP4.5 climate scenarios (Table S1 in the supporting information and Figure 6). The projected AGB was found to increase in most sites when using the averaged τ w , but with large uncertainties caused by τ w . One exception was observed for a boreal coniferous forest in Canada, where the AGB was consistently decreasing during all the test runs ( Figure 6). The simulated AGB was shown to be sensitive to τ w for all sites, resulting in a large variation in τ w by the year of 2100. The projected AGB for Manaus and BANNA was found to increase fast when the averaged τ w was used and a twofold AGB was obtained by 2100 for BANNA. The four temperate PFT sites witnessed similar simulation results and an increasing trend in AGB by 2100 for the averaged τ w case, indicating that the temperate forests should benefit from the projected climate change in the future. Results from two of four boreal sites (US-WCr and US-Syr) revealed a generally constant AGB over the simulation period and a decreasing AGB for the remaining two of all τ w cases. However, our simulation results also show that the boreal forest may benefit from future climate change with an increasing AGB trend in certain test cases. This further necessitates the accurate estimation of τ w for model simulation. The green lines show the 1000 test runs using the random τ w data resulting between the one fourth and three fourth percentiles; the red line shows the result of the averaged τ w based on the collected field data. The abbreviations are defined as follows: TrB, tropical broadleaf forest; TeC, temperate coniferous forest; TeB, temperate broadleaf forest; BoC, boreal coniferous forest; BoB, boreal broadleaf forest. All the test sites were randomly selected from Fluxnet; details of the sites are provided in Table S1.

10.1002/2016GB005557
Difference in the reproduced historical variations of AGB increased for the three model runs from 2006 to 2100 (Figure 7). The two improved runs with our estimated τ w showed a consistent increasing trend for AGB, while the model run with default τ w results in a decreasing trend and AGB does not change much after 2060. Though similar results were found for the two improved model runs, simulated AGB with present-day τ w was slightly larger than that with projected annual τ w . The simulated global forest AGB was 171.1 and 170.0 Mg ha À1 in 2100 for the two cases, respectively. Because the projected τ w showed a decreasing trend by 2100 ( Figure S2), assuming a constant τ w may overestimate AGB for long-term simulations. Furthermore, the simulated spatial patterns of AGB in 2100 under the RCP4.5 scenario also showed large differences for the three runs ( Figure S3). AGB from the two model runs with estimated τ w is much larger in tropical forests and temperate forests areas compared with results with default τ w . Another major difference was found for the transition from boreal coniferous evergreen forest to broadleaf cold deciduous forests in the Siberian areas, indicating the essential role of τ w in mapping PFT over long-term simulation ( Figure S4). Actually, the simulated forest PFT mapping showed large difference among the three simulation runs ( Figure S4). Furthermore, the climate changes from cold and dry prior to 2040 to warm and wet later in the projected scenario. The interactions of changes in climate and the different woody residence values may result in the overall different variations in AGB for the simulation periods.

Abiotic and Biotic Influences on Woody Residence Time
Few studies have focused on the mapping of τ w at a global scale. A recent study investigated the global carbon residence time using field data and model simulations [Carvalhais et al., 2014]. Because carbon residence time used in that study is different from τ w in the present study in that the former calculates the residence time for both soil and vegetation residence time (including leaf, stem, and litter), comparing their values and our results directly would be difficult. However, their results show carbon residence time increases as temperature decreases in forests located in tropical to boreal areas, which is different from our results. This may be because decomposition occurs more slowly in boreal areas than in tropical areas because of the low temperature, resulting in the largest carbon residence time in boreal forests [Carvalhais et al., 2014]. In our case, we did not consider soil carbon residence time and this could result in relatively small τ w for boreal forests compared with other PFTs. This difference indicates that when calculating total ecosystem carbon stocks for both vegetation and soil, modelers should pay more attention to carbon residence time in soil of ecosystems in cold areas [Yang et al., 2014;Chen et al., 2015].
Globally, we observe a negative influence of P on τ w according to LMM analysis; while T has positive influence on τ w . However, the influence of P and T on τ w may be different for a given PFT ( Figure S5). T is found to be positively related to τ w for boreal (p = 0.54) and temperate PFTs (p < 0.05), while it is negatively related to τ w for warm temperate (p = 0.63) and tropical PFTs (p < 0.05). In cold areas, the forest may suffer from frost damage, thus resulting in large forest mortality [Thurner et al., 2016]. Therefore, forest mortality in these areas is more controlled by temperature because of the inhibition of carbon assimilation for the low temperature and short grown season [van Dijk et al., 2005;Beer et al., 2010]. Thurner et al. [2016] found that the forest carbon residence time for temperate forests in Northern Hemisphere was more controlled by P, which is partly Global Biogeochemical Cycles 10.1002/2016GB005557 supported by our results ( Figure S5). We found a quadratic relationship between P and τ w for temperate forest. In detail, τ w is found to decrease with P when P is below 1500 mm, while it increases with P when P is above 1500 mm. Most our collected temperate forest sites with P above 1500 mm are located in west North America with Mediterranean climate (Figure 1). The forests in these areas are found with large mortality in recent years most probably due to droughts [van Mantgem and Stephenson, 2007;van Mantgem et al., 2009]. Even though the mechanism of forest mortality due to droughts is still controversial, drought-induced forest mortality has been widely observed globally [Allen et al., 2015]. The temperate forest plots with P below 1500 mm are mostly located in coastal areas of Europe with high latitude and mountainous areas of China (Figure 1). Forest growth in these areas may be limited by the surplus water due to the limited evaporation . Thus, our results show that determinants of τ w may be different even for the same PFT depending on the local climate [Carvalhais et al., 2014;Malhi et al., 2015;Thurner et al., 2016]. Furthermore, several studies show that accumulation of woody biomass was more controlled by mortality due to extreme disturbance events instead of normal meteorological conditions Thurner et al., 2016]. Therefore, the relationships using only T and P as predictors may not be enough for the global estimation of τ w (Table 1 and Figure S6).
Our LMM analysis shows the important influence of interaction term between biotic term (i.e., age) and T (or P) on τ w . Over large scales, τ w may be determined more by forest structure and ecosystem succession [Kobe et al., 1995;Purves et al., 2008], climate [Allen et al., 2015], and responses to climate change and disturbances [Kurz et al., 2008]. Therefore, the estimation of forest mortality or τ w , especially at a large scale, is difficult and also has large uncertainty [e.g., Carvalhais et al., 2014;Thurner et al., 2016]. Malhi et al. [2015] suggested that a direct link may exist between forest mortality (i.e., τ w ) and carbon use efficiency and/or woody allocation. Based on our collected sites available with CUE, we found inverse relationships between CUE and τ w ( Figure S6). Further analysis also shows a significant inverse relationship between woody allocation and τ w for tropical forests (n = 18, data not shown). This indicates that the forests with a high CUE and W p are usually characterized with a short life time, and therefore, AGB cannot be simply predicted by CUE or W p [Malhi et al., 2015]. The apparent relationship between high CUE and W p with τ w is not a direct result of variation in the productivity (W p ); instead, they may be a result of forest ecosystem trade-offs between growth and defenses against large-scale disturbance such as wind or disease (along edaphic gradients) or perhaps reproduction and persistence (along climatic and edaphic gradients) by the surrogate of W p [Stephenson and van Mantgem, , 2011. This has been partly verified in Amazonian forests, where an edaphic decreasing gradient exists from west to east [Baker et al., 2004]. The eastern Amazonian forests growing on infertile soil (thus, with less subcanopy vegetation), more recourses are spent on forest defense while low productivity and slow growth favor a large τ w ; in contrast, the counterparts in western Amazonian forests on fertile soil (with more subcanopy) exhibit much higher productivity and rapid growth, which may result in a small τ w (Figure 3b). Another research by Quesada et al. [2012] also shows that soil fertility plays an essential role in the dynamics of forest structure and biomass turnover in Amazonian forests by driving the patterns of mortality in the area.

Estimation of Global Forest Woody Residence Time
Based on the collected field data, we mapped the global woody residence time by a RF method (Figure 4). To our knowledge, this is the first τ w map at the global scale, and therefore, this map may help improve model simulation of AGB by parameterizing τ w using this data set (Figures 5 and 6). However, there is still uncertainty in the simulations (Figure 5b), which may originate from absence of plots in certainty areas (such as Siberia), the random forest method, and model structure in the description of forest growth/mortality. The LMM analysis shows the important influence of forest age and CUE on τ w , while these two variables could not be integrated as predictors in the RF method because gridded forest age and CUE data at the global scale were not available yet. Therefore, large-scale estimation of these two variables should be to conducted in future to improve the calculation of τ w in the RF method.
A recent research focusing on boreal and temperate forest residence time was published based on remote sensing data [Thurner et al., 2016]. Thurner et al. [2016] calculated forest residence time for both above ground and below ground, and it is therefore difficult to compare their results with ours directly. Following the methods in Thurner et al. [2016], we downloaded AGB from Thurner et al. [2014] and also downloaded the global MODIS NPP and calculated τ w for the boreal and temperate forests at 0.01°× 0.01°resolution Global Biogeochemical Cycles 10.1002/2016GB005557 ( Figures S7 and 8). Similar patterns were found for both of the two results, such as the relatively large values in Europe, northwestern United States, and Siberia ( Figure 8). On the other hand, τ w from our study is larger than that from Thurner et al. [2016] in most of the global forest areas (Figure 8b). This may be due to the fact that τ w from our study was extrapolated from plot values assuming that forests were near equilibrium and the plots used for random forest calculation are largely undisturbed. Actually, secondary and/or human planted forests are found in North America and China and AGB for these forests were usually below equilibrium values [Liu et al., 2014;Galbraith et al., 2013]. We calculated τ w based on their results by using a simple regression model to estimate aboveground NPP from MODIS NPP ( Figure S8). This calculation processes may introduce uncertainty to the resulting τ w and also make it difficult for the comparison. Nevertheless, with more remote sensing data of AGB and/or NPP becoming available for regional or global scale [e.g., Thurner et al., 2014;Hu et al., 2016], the method in Thurner et al. [2016] could have large potential in estimation of global forest τ w .

Implications for Model Simulation
The τ w can induce high levels of uncertainty in model simulation of AGB, both at site and global levels. This has already been shown in other model simulations as well as in field data-based analysis [Delbart et al., 2010;Castanho et al., 2013;Friend et al., 2014;Malhi et al., 2015]. Based on regional estimation of τ w in Amazonian forests, Castanho et al. [2013] found that the spatially explicit τ w could improve the accuracy of modeled AGB. Similarly, a large difference in estimated AGB was observed by different parameterization schemes of τ w in our results (Figure 7).
DGVMs integrated in Earth system models (ESMs) are useful tools in the projection of future carbon fluxes and stocks at a global scale [e.g., Friedlingstein et al., 2006;Sitch et al., 2008]. Default values of τ w for PFTs in DGVMs should be further calibrated based on additional field data ( Table 2). Eight of 12 DGVMs use a fixed background turnover rate (or τ w ) and fire-induced mortality to simulate the overall τ w . Even though some Global Biogeochemical Cycles 10.1002/2016GB005557 models explicitly simulate mortality, these are still with uncertainty due to the unclear mechanisms in mortality per se [McDowell et al., 2011]. No background τ w is provided in Organizing Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE) and Community Land Model (CLM)-DGVM, and these two models explicitly simulate fire disturbance mortality modified from Lund-Potsdam-Jena (LPJ), which underestimates the mortality rate over a large scale. This results in overestimated τ w and thus carbon stock, especially for tropical forests [Krinner et al., 2005]. The Ecosystem Demography version 2 (ED2) model provides a density independent mortality rate for PFTs during different successional stages without providing an explicit background turnover rate. Only one of 12 models (i.e., Jena diversity (JeDi)-DGVM) investigated actually simulated the background τ w using dynamic values; in this case it was determined by air temperature [Pavlick et al., 2013]. The fixed background τ w varies considerably even for the same PFT for the investigated models. For example, the τ w for boreal coniferous forests ranges from 1.5 to 500 years for all the models; many default values are far from the averaged ones found in our meta-analysis. Furthermore, one third of the 12 models use the same values of τ w for all PFTs in their simulations. This indicates that large variations of AGB may be derived in model simulations even with the same meteorological variables [Sitch et al., 2008;Friend et al., 2014]. Furthermore, most of the current DGVMs only consider background mortality and influences of natural disturbances on τ w and do not explicitly simulate human disturbances. Further uncertainties in AGB may also arise from assumptions of the homogeneity of forest stands (such as equal probability of forest mortality for all species) and the steady state in the DGVM calculations. Therefore, even though our estimated gridded τ w may help in the parameterization of models at this stage, more work should be done to explore the forest mortality mechanisms in the future [McDowell et al., 2011] and thus to improve model structure in the description of recruitment and mortality.
Several studies observed the lessened τ w during the most recent decades [Phillips and Gentry, 1994;Phillips et al., 2004], which is different from the fixed temporal τ w in models ( Figure S2). The change in τ w was, to a large extent, caused by climate change such as the warming temperature and/or the increasing atmospheric CO 2 concentration [Phillips and Gentry, 1994;Phillips et al., 2004]. The ESMs projected enhanced forest productivity in the current century mainly caused by global warming and CO 2 fertilization [e.g., Cramer et al., 2001;Sitch et al., 2008;Piao et al., 2013]. Therefore, an invariant τ w may result in an overestimation of simulated AGB according to equation (1) if τ w decreased in the future. Our simulation results show that assuming a time invariant τ w may result in larger AGB compared with using the projected τ w (Figures 7 and S3). Therefore, a better understanding and prediction on the responses of τ w to climate change and atmospheric CO 2 concentration is needed and should be fully investigated.

Conclusions
The τ w is an important parameter that expresses the balances between mature forest recruitment/growth and mortality. By collecting field data from the literature, this study explores the global forest woody residence time and investigates its influences on model simulations of AGB at the global scale. Parameter τ w is highly related with forest age, T and P, but shows different determinants among various PFTs. The estimated global forest τ w based on the collected field data shows large spatial heterogeneity as well. This heterogeneity plays an important role in model simulation of carbon stocks in DGVMs. Parameter τ w could change the resulted AGB in 10 folds based on a site-level test using the Monte Carlo method. At the global level, different parameterization schemes of the Integrated Biosphere Simulator using estimated τ w resulted in a twofold change in AGB for 2100. The results from the study highlight the influences by various biotic and abiotic variables on forest τ w . The global estimation of τ w may help improve the model simulations and lower the parameter uncertainty over the projection of future AGB in current DGVM or ESM models. A clearer understanding of τ w responses to climate change is also needed for the improvement of model prediction of carbon stock in future studies.