Data-driven quantification of nitrogen enrichment impact on Northern Hemisphere plant biomass

The production of anthropogenic reactive nitrogen (N) has grown so much in the last century that quantifying the effect of N enrichment on plant growth has become a central question for carbon (C) cycle research. Numerous field experiments generally found that N enrichment increased site-scale plant biomass, although the magnitude of the response and sign varied across experiments. We quantified the response of terrestrial natural vegetation biomass to N enrichment in the Northern Hemisphere (>30° N) by scaling up data from 773 field observations (142 sites) of the response of biomass to N enrichment using machine-learning algorithms. N enrichment had a significant and nonlinear effect on aboveground biomass (AGB), but a marginal effect on belowground biomass. The most influential variables on the AGB response were the amount of N applied, mean biomass before the experiment, the treatment duration and soil phosphorus availability. From the machine learning models, we found that N enrichment due to increased atmospheric N deposition during 1993–2010 has enhanced total biomass by 1.1 ± 0.3 Pg C, in absence of losses from harvest and disturbances. The largest effect of N enrichment on plant growth occurred in northeastern Asia, where N deposition markedly increased. These estimates were similar to the range of values provided by state-of-the-art C–N ecosystem process models. This work provides data-driven insights into hemisphere-scale N enrichment effect on plant biomass growth, which allows to constrain the terrestrial ecosystem process model used to predict future terrestrial C storage.


Introduction
In the last century, humans have drastically altered the global nitrogen (N) cycle by producing reactive N and spreading it over ecosystems. Reactive N inputs come from fertilizers synthesis by the Haber-Bosch process, N oxides produced by fossil fuels and biofuels, and the cultivation of N 2 -fixing crops (Vitousek et al 1997, Galloway et al 2004, Canfield et al 2010, Peñuelas et al 2020.
Here we focus on the impact of increasing N deposition (Dentener et al 2006, Ackerman et al 2019, on terrestrial ecosystems (Hietz et al 2011, Fowler et al 2013. Previous studies found that most terrestrial ecosystems were limited by N availability, particularly in the Northern Hemisphere at mid-and high latitudes (Elser et al 2007, LeBauer and Treseder 2008, Peñuelas et al 2013, Craine et al 2018, Du et al 2020. Clarifying the effect of atmospheric N input on the growth of terrestrial plants is thus critical to understand terrestrial carbon (C) storage, i.e. how much of the current land C sink is caused by atmospheric deposition (Gruber and Galloway 2008, Reay et al 2008, Schulte-Uebbing et al 2022.
Field N-enrichment experiments have been conducted in various terrestrial ecosystems during the last four decades ( figure 1(a)). The data from such site-scale field experiments were and continue to be used to explore the effects of N enrichment on ecosystem C cycling (LeBauer and Treseder 2008, Janssens et al 2010, Song et al 2019, Du et al 2020. For instance, meta-analyses of N-enrichment experiments showed that the mean effect of N enrichment on site-scale biomass was positive (LeBauer and Treseder 2008, Yue et al 2017, Schulte-Uebbing and de Vries 2018, Song et al 2019. However, the effect of N enrichment on biomass varies drastically across experiments, due to local conditions such as climate, vegetation, background soil fertility, N-enrichment intensity and duration, and experimental design (Xia and Wan 2008, Stewart 2010, Greaver et al 2016. This large variability poses a substantial challenge to the data-driven quantification of regional-or globalscale responses of terrestrial biomass to elevated N deposition. We quantified the effects of climate, soil characteristics, and N-enrichment intensity on the response of Northern Hemisphere vegetation biomass to N enrichment using two machine-learning algorithms. Given that there may be divergence in plant aboveand below-ground adjustment strategies under resource stress (Freschet et al 2018, Tumber-Dávila et al 2022, we investigated the dominant source of variation in N enrichment effects in both aboveground biomass (AGB) and belowground biomass (BGB). The machine-learning algorithms were trained using data from peer-reviewed N-enrichment experimental studies (see Methods). In total, we compiled 597 observations of the response of AGB to N addition from 100 sites, and 176 observations of the response of BGB from 42 sites (table S1, supplementary data). All the N-enrichment experiments were paired with a control and a treatment at the same location. The AGB and BGB responses cumulate effects of N addition throughout the experimental period.
The data on AGB and BGB responses covered a range of vegetation types (figure 1(a) and table S1) and intensity of N enrichment from 0.2 to 64 g N m −2 y −1 at AGB sites (figure 1(b)) and from 1 to 56 g N m −2 y −1 at BGB sites (figure 1(c)). The atmospheric N deposition change during 1993-2010 reaches at highest 1.3 g N m −2 y −1 (figure 1(a)). Thus, N enrichment intensity in the experiments somewhat covered the change in atmospheric N deposition, albeit with a significant bias towards the higher end and beyond. In this study, in addition to metaanalyses approach, we used boosted regression trees (BRT) (Elith et al 2008) and random forest (RF) models (Breiman 2001) to predict the responses of AGB and BGB to the intensity of N-enrichment. The conclusions from the two machine-learning algorithms were consistent, so we present the results obtained with BRT in the main text (and the RF results in supplementary information).

Data set of N-enrichment experiments
We collected data for AGB and BGB in N-enrichment experiments from four meta-analyses: Song et al (2019), Schulte-Uebbing and de Vries (2018), Yue et al (2017), andTian et al (2016). We obtained the data in Song et al (2019) from the authors and the data in the other three studies that were not included in Song et al (2019) from references provided therein. We used data from N-enrichment experiments in natural terrestrial ecosystems between 30-90 • N. The median experimental duration was three years (figure S1). We collected a total of 597 records (including replicates and years in each site) for the response of AGB to N addition from 100 sites and a total of 176 records (including replicates and years in each site) for the response of BGB from 42 sites (table S1, figure 1(a), supplementary data).

Meta-analysis of observed effects of N enrichment on above-and belowground biomass
The effect of N enrichment on AGB (or BGB) likely varies across the N-enrichment experiments due to the spatial heterogeneity in climatic, soil, and experimental characteristics. Thus, in the meta-analysis, we used random-effects models assuming that the effects being estimated in the different studies are not identical, but follow some distribution representing the between-study variability (Gurevitch et al 2018). We conducted the meta-analysis using the 'escalc' and 'rma.uni' functions in the 'metafor' package of R software (Viechtbauer 2010). Specifically, the effects of N on AGB and BGB were measured by estimating the mean response ratio RR = ln , where X t and X c are mean biomasses in the N-enrichment and control treatments, respectively (Hedges et al 1999, Lajeunesse 2011. This was performed by setting the parameter 'measure' as 'ROM' in the 'escalc' function in the 'metafor' package of R software (Viechtbauer 2010). The weighted response ratio (RR w ) was calculated as the weighted average of RR using the weights ω i = 1/(ν i + τ 2 ), where ν i is the variance of the effect size within the ith study and τ 2 is the between-study estimated by a restricted maximum-likelihood estimator (Viechtbauer et al 2015). Parameter estimation was performed by setting the parameter 'method' as 'REML' in the 'rma.uni' function in the 'metafor' package of R software (Viechtbauer 2010). The percent changes of AGB and BGB due to N enrichment were calculated as [exp (RR w ) − 1] × 100%. The effects of N enrichment on AGB and BGB were considered to differ significantly between the N-enrichment and control treatments when the 95% confidence intervals of ∆AGB and ∆BGB did not include zero.

Observation based sensitivities of AGB and BGB to N enrichment
The sensitivities of AGB and BGB to N enrichment (S AGB and S BGB ) were calculated as: where S AGB is the relative response of AGB to N enrichment (% [g N m −2 y −1 ] −1 ), ∆AGB is the N enrichment-induced relative change in AGB over the experimental period (%), and N add is the amount of N added in the treatment plots (g N m −2 y −1 ): where S BGB is the relative response of BGB to N enrichment (% [g N m −2 y −1 ] −1 ), ∆BGB is the N enrichment-induced relative change in BGB over the experimental period (%), and N add is the amount of N added in the treatment plots (g N m −2 y −1 ).

Relative influence of climatic, soil, and experimental characteristics on S AGB and S BGB
The spatial variations of S AGB and S BGB were examined using BRT (Elith et al 2008) and RF models (Breiman 2001). We conducted the BRT analyses using the 'gbm.step' function in the 'gbm' package of R software, with the parameters of 'tree.complexity' as 5 and 'learning.rate' as 0.005. The RF analyses were conducted using the 'randomForest' function in the 'randomForest' package of R software, with the parameters of 'nodesize' as 5 and 'ntree' as 500.
The BRT and RF models were trained using 16 predictor variables: climatic variables (mean annual temperature (MAT) and mean annual precipitation (MAP)), woodiness (woody or nonwoody), foliar N content, soil characteristics (C:N ratio, bulk density, pH (measured in water), cation exchange capacity (CEC), and the contents of organic C, clay, organic phosphorus (P), labile P, and water), and experimental characteristics (AGB and BGB in the control plots, intensity of N addition, and treatment duration). We obtained data for MAT and MAP from the WorldClim2 database (Fick and Hijmans 2017).
To ensure the comparability of N deposition data between BRT and RF models and Multi-scale Synthesis and Terrestrial Model Intercomparison Project (MsTMIP) models, we systematically used the same data set for N deposition (Wei et al 2014a(Wei et al , 2014b. We extracted data for the soil C:N ratio, bulk density, pH, CEC, and the contents of organic C and clay from the gridded Global Soil Dataset for use in Earth System Models (GSDE) (Shangguan et al 2014), and from the WISE30sec database (ISRICWISE) (Batjes 2016). Data for soil-water content were extracted from GSDE (Shangguan et al 2014). Data for the contents of soil organic P and labile P were extracted from Global Gridded Soil Phosphorus Distribution Maps at resolutions of 0.5 • (Yang et al 2014). Data for foliar N content were extracted from global maps of the distributions of plant traits (Butler et al 2017). Data for woodiness (woody or nonwoody) were extracted from the Global Mosaics of the standard Moderate Resolution Imaging Spectroradiometer (MODIS) land-cover type data product (MCD12Q1) in the International Geosphere-Biosphere Programme (IGBP) land cover type classification (Friedl et al 2010). The forests, shrublands, and woody savannas were defined as 'woody' , and the other vegetation types were defined as 'nonwoody' . We used the S AGB and S BGB samples with all 16 variables in the BRT and RF analyses. The data sets for climate, woodiness, N deposition, and soil characteristics were also used in the spatial extrapolation of S AGB and S BGB (figure S2, see below). The machine learning analysis was performed 100 times to examine the relative influence of each predictor of the 16 predictors on the S AGB and S BGB . Here, relative influence of a predictor in BRT analysis is relative contribution of the variable for a BRT model, which was 'contributions' parameter outputted by 'gbm.step' function in the 'gbm' package of R software. Partial-dependence plots for the variables in BRT models were produced using 'gbm.plot' function. Variable importance in RF analysis was assessed using the total decrease in residual sum of squares from splitting regression tree on the variable, which was 'IncNodePurity' parameter in 'importance' object outputted by 'randomForest' function in the 'randomForest' package of R software. Partialdependence plots for the variables in RF models were produced using 'partialPlot' function. The 16 predictors were ranked by the value of their influence on the S AGB and S BGB from high to low. Then, a series of machine learning models including 2-16 predictors were established to examine the performance of simpler models. For each machine learning model, 10-fold cross-validation was used to test the proportion of variance of S AGB (or S BGB ) explained by S AGB (or S BGB ) predicted by the models (R 2 ). The crossvalidation were performed 100 times with the average results shown in the figures. The machine learning model with highest R 2 was used in spatial extrapolation of S AGB and S BGB in the following section.

Spatial extrapolation of S AGB and S BGB
We calculated the spatial distributions of S AGB and S BGB at mid-and high latitudes (30-90 • N) of the Northern Hemisphere using both the BRT and RF models trained by site data and of the gridded climatic and soil variables, with the treatment duration set as 17 years from 1993 to 2010, the intensity of N enrichment as the average change in N enrichment during 1993-2010 relative to 1993, and AGB and BGB in 1993. In the spatial extrapolation analysis, the first year of the treatment duration was set as 1993, because 1993 was the first year of the dataset for global AGB (Liu et al 2015) used in this study. The last year of the treatment duration was set as 2010, because 2010 was the last year of the duration of the MsTMIP models' simulations (Wei et al 2014a) used in this study. BGB was calculated as AGB multiplied by the BGB:AGB ratio (R B2A ) (Liu et al 2015), with the source noted as 'Liu' in the figures. We also used R B2A and total biomass reported by Carvalhais et al (2014) to calculate global AGB and BGB, with the source noted as 'Carvalhais' in the figures. The GEO-CARBON global forest AGB (Santoro et al 2015, Avitabile et al 2016) was also used for the spatial extrapolation of S AGB and S BGB , with BGB calculated as AGB multiplied by R B2A , with the source noted as 'GEO-CARBON in the figures. The relative changes in terrestrial AGB and BGB caused by N deposition during 1993-2010 (∆AGB and ∆BGB) were calculated as S AGB and S BGB multiplied by the average change in annual N deposition during 1993-2010 using driver data (N deposition) of MsTMIP (Wei et al 2014b).
2.6. Total change in terrestrial biomass due to enhanced atmospheric N deposition during 1993-2010 for each grid point BRT-and RF-based change in total biomass (∆TB) were calculated for each grid point using ∆AGB and ∆BGB with AGB and BGB as weights: where ∆TB is the percent change in total biomass due to N enrichment during 1993-2010, ∆AGB is the percent change in AGB due to N enrichment during 1993-2010, ∆BGB is the percent change in BGB due to N enrichment during 1993-2010, AGB 1993 is AGB in 1993 (g C m −2 ), and BGB 1993 is BGB in 1993 (g C m −2 ).

∆TB in Northern Hemisphere
∆TB in Northern Hemisphere (30-90 • N) was calculated as: where ∆TB is the percent change in total biomass due to N enrichment during 1993-2010, i indicates grid cell i (0.5 • × 0.5 • ), n indicates the number of grid cells, ∆TB i is the percent change in total biomass due to N enrichment during 1993-2010 at grid cell i, TB 1993i is total biomass in 1993 at grid cell i (g C m −2 ), Area i is the area of grid i (m 2 ), and TB 1993NH is total biomass in 1993 (g C), calculated as

Terrestrial ecosystem process model simulations
We used total biomass from six terrestrial ecosystem process models with C-N interactions: where ∆TB is the percent change in total biomass due to N enrichment during 1993-2010, TB 1993 is total biomass in 1993 (g C m −2 ) induced by N deposition, and TB 2010 is total biomass in 2010 (g C m −2 ) induced by N deposition. The relative change in Northern Hemisphere plant biomass due to N enrichment the MsTMIP models was calculated using equation (4) as the change in total biomass induced by N deposition during 1993-2010.

Results
Meta-analyses of our dataset revealed that N enrichment on average increased both AGB and BGB in field experiments. AGB was higher by 30 33 27 % (mean and 95% confidence interval) (figure 1(d)) and BGB by 11 14 7 % (figure 1(e)) as compared to each control experiment. The dominant factors influencing the sensitivities of AGB and BGB to N enrichment (S AGB and S BGB , the percent of biomass change over the experimental period due to N enrichment, in % [g N m −2 y −1 ] −1 ) were deduced from the machine-learning algorithms. 16 predictor variables were considered (see Methods). Based on these predictors, the BRT models were able to explain 56%-57% of the variance in S AGB and ∼20% of the variance in S BGB based on our leave-one-out crossvalidation (figure 2). The lower performance for BGB was probably due to the lower amount of data available for this variable. The ranges of climatic conditions and soil properties at the experimental sites of N addition covered those observed in terrestrial ecosystems at Northern Hemisphere mid-and high latitudes, indicating the representativeness of the climatic and soil conditions at the experimental sites (table S2).
The intensity of N enrichment had the largest influence on S AGB based on the BRT models (figures 2(a) and (b)). S AGB decreased with N enrichment (figures 3(a), (b), S3 and S4) from ∼15% [g N m −2 y −1 ] −1 ) at a N input of Figure 2. Relative influences of climatic, soil, and experimental characteristics on SAGB and SBGB. (a) and (b) Relative influence on SAGB using BRT, with soil characteristics from the GSDE and ISRICWISE databases, respectively. The relative influence is shown as the mean of the output from 100 BRT analyses. R 2 represents the proportion of variance of SAGB (or SBGB) explained by SAGB (or SBGB) predicted by BRT based on 2-16 predictors. For each BRT model, the cross-validation were performed 100 times with the averaged R 2 shown in the subfigure. The red vertical line shows the BRT model with the largest value of R 2 . (c), (d) Same as (a) and (b), but for SBGB. 0.2 g N m −2 y −1 , to a value close to zero at a very high N input of 64 g N m −2 y −1 . This was consistent with the results of a regression analysis, which further revealed that S AGB remains close to zero beyond a critical N enrichment intensity of approximately 10 g N m −2 y −1 (figures S5 and S6). N enrichment intensity was also identified as the largest most influential factor on S AGB by the RF models (figures S7-S9). However, the second most influential factor differed between the two approaches: AGB in the control plot in BRT (figures 2(a), (b), S3 and S4) and soil labile P content (or soil organic P content) in RFs (figures S7-S9).
Unlike S AGB , the dominant source of variation in S BGB was BGB in the control plot (figures 2(c) and (d)), while the intensity of N enrichment was ranked second and also had a strong impact on S BGB (figure S5(c)). S BGB nonlinearly decreased as a function of background BGB (figures 3(c), (d), S5(d), S10, and S11). This result was also evident in the raw observation data ( figure S6). The results of the BRT models were generally consistent with those obtained with the RF models results (figures S5, S7, S12 and S13).
In addition to the analysis using all 16 variables, we also performed the analysis with fewer (2-15) variables to establish simpler machine-learning models. BRT models still explained the highest proportion of S AGB variance when using seven most influential variables including soil data from GSDE data set (figure 2(a)) or 13 most influential variables including soil data from ISRICWISE data set ( figure 2(b)), and explained S BGB with only four influential variables (figures 2(c) and (d)). This was Figure 3. Partial-dependence plots of the predicted SAGB and SBGB with independent predictor variables. (a) and (b) SAGB response to climatic, soil, and experimental variables variation using BRT with soil characteristics from the GSDE and ISRICWISE databases, respectively. The range of predictor variables was normalized to 0-1. (c), (d) Same as a and b, but for SBGB. Figures S3, S4, S10 and S11 show the partial-dependence plots with independent predictor variables expressed as absolute value.
The extrapolation analysis shows that AGB generally increased from atmospheric N deposition across the Northern Hemisphere during 1993-2010 (figures 4(a) and S14(a)-(f)), as constrained by observations from the N-enrichment experiments ( figure 1(d)). N enrichment had minimal impacts on BGB, but a strong positive effect on AGB changes, that is ∆AGB (figures 4(a), (b) and S14(a)-(l)). This is likely because N enrichment intensity was the most influential factor of the variation of S AGB , but not of S BGB , in the machine-learning models used for upscaling extrapolation. We further analyzed the spatial distribution of the effect of increased N deposition on the change in total biomass (∆TB) during 1993-2010 by combining the responses of ∆AGB and ∆BGB (Methods). Total biomass was generally enhanced by changes in N deposition across Northern Hemisphere, reflecting the changes in AGB (figures 4(c) and S14(m)-(r)). Biomass increased the most in northeastern Asia, where the largest enhancement of atmospheric N deposition occurred during 1993-2010. The spatial patterns of ∆AGB, ∆BGB, and ∆TB were consistent between the BRT and RF models (figures S14 and S15).
We further investigated ∆TB caused by N enrichment during 1993-2010 using six terrestrial ecosystem process models that include C-N interactions: CLM4, CLM4VIC, DLEM, ISAM, TEM6, and TRIPLEX-GHG from the MsTMIP (Mao et al 2015, Huntzinger et al 2013, Wei et al 2014a. All these models outputted the response of total biomass to N enrichment but did not distinguish between AGB and BGB. None of the models produced the same spatial patterns of ∆TB as those from the machinelearning approaches (figures 4(d) and S16), although the MsTMIP models and the BRT and RF models were all driven by the same N-deposition data set (see methods). The spatial pattern of ∆TB varied greatly across the six MsTMIP models (figure S16). N enrichment generally had a minimal effect on total biomass in ISAM. Total biomass in the other five models responded positively to N enrichment, albeit to different degrees. The multimodel mean indicated that N enrichment increased total biomass in eastern North America, Europe, and eastern Asia (figures 4(d) and (e)), but with large spread across the terrestrial ecosystem process models (figure S16).

Discussion and conclusion
Our study based on machine-learning algorithms indicated that the AGB increase responded nonlinearly to the intensity of N addition, consistent with previous meta-analyses (Arens et al 2008, Bradford et al 2008, Ochoa-Hueso 2016, Tian et al 2016, Prager et al 2017, Xu et al 2018. However, the N enrichment applied in field experiments was typically much higher than the background atmospheric N deposition (figures 1(a)-(c) and S22). Therefore, the overall mean N response estimated by meta-analyses may not accurately represent the larger-scale mean effect of increased N deposition when the strongly nonlinear responses to N addition are not accounted for. In contrast to the meta-analyses, our integrated analyses based on machine-learning approaches did consider the stronger effect of N-enrichment at low doses. The intensity of N enrichment was the dominant cause of terrestrial S AGB variations, when this sensitivity was derived from observations in N-enrichment experiments (figures 2 and S7), with S AGB decreasing from 15% [g N m −2 y −1 ] −1 to nearly 0% [g N m −2 y −1 ] −1 as N enrichment increased from 0.2 to 64 g N m −2 y −1 , (figures S3 and S4). The apparent difference in S AGB between our study and previous meta-analyses therefore indicates that considering realistic N enrichment intensity is recommended in future field experiments and metaanalyses studies focusing on N effect on terrestrial C cycling.
Previous meta-analyses have generally suggested positive effects of N addition on AGB (Yue et al 2016, You et al 2017, Schulte-Uebbing and de Vries 2018 and BGB (Li et al 2015, Yue et al 2016 in N-enrichment experiments. However, when setting N addition as a change in N deposition during the same duration of 1993-2010 and the identical experimental duration in machine-learning approaches, we found that N enrichment only had a minor effect on Northern Hemisphere terrestrial BGB (figure 4(b)), indicating that N enrichment generally decreased the BGB:AGB ratio in Northern Hemisphere terrestrial ecosystems. This finding is consistent with changes in C allocation due to N enrichment in N-limited terrestrial ecosystems, favoring the allocation of C to AGB instead of roots (Chapin 1980, Müller et al 2000, Makela et al 2008, Cambui et al 2011, Yue et al 2021, Peng et al 2022. The positive effect of N enrichment on total biomass was dominated by the response of AGB rather than BGB in the Northern Hemisphere terrestrial ecosystems, due to the different responses of AGB and BGB to N enrichment (figures 4(a)-(c)).
N enrichment-induced decrease of the BGB:AGB ratio maybe related to plant adaptation strategies in changing nutrient conditions. Plant biomass allocation changes due to N enrichment were usually assessed using the two classic mechanisms based on the optimal partitioning hypothesis and the isometric allocation hypothesis, respectively. Under the optimal partitioning hypothesis, biomass is preferentially allocated to the organs that could acquire the most limited resource for plant growth (Bloom et al 1985, Kobe et al 2010. For instance, plants preferentially allocate more biomass to root under N starvation but allocate more biomass to shoot under N enrichment (Mardanov et al 1998, Kobe et al 2010, Chen et al 2013. Under the isometric allocation hypothesis, the biomass allocation is allometric among plant organs but is isometric across various environmental conditions, plant species or vegetation types (Niklas 2004(Niklas , 2005. Plants allocate biomass to each organ following scaling exponents based on individual plant size (Cheng and Niklas 2007). The integrated analysis of N enrichment experiments observations showed that N enrichment decreased plant root:shoot ratio but did not apparently change the allometric relationships among plant organs when the whole set of data from various ecosystems were considered (Peng and Yang 2016, Yue et al 2021, Peng et al 2022. Nevertheless, as shown in figure 5 of Peng et al (2022), there was large uncertainty in the allometric scaling exponents among plant organs with 95% confidence interval ranging from ∼0 to ∼2 under both control and N enrichment conditions. Stronger evidence is still needed to clarify whether the allometric relationships are independent of nutrient conditions and vegetation types. Particularly, a higher number of paired data for AGB and BGB response to N enrichment in each vegetation type is warranted for identifying the mechanism that can most accurately explain the BGB:AGB ratio decrease caused by N enrichment.
The Northern Hemisphere terrestrial ∆TB caused by N enrichment during 1993-2010 varied among the process-based models simulating terrestrial C cycles with C-N interactions (figure 5). ∆TB was under-or overestimated by most models relative to the estimates from the machine-learning approaches. Different representations of the framework of N cycles in C-N models likely leads to great uncertainty in modeling C cycles (Niu et al 2016, Du et al 2018).
Several key mechanisms of C-N cycles remain to be improved in state-of-the-art terrestrial ecosystem process models, such as community composition, contents of labile C and N, allocation and turnover of C and N pools, biological N fixation, and losses of N from the ecosystem via leaching or gaseous emissions (Thomas et al 2015). Our study found that the mismatch between C-N model simulations and observations was widely distributed across Northern Hemisphere terrestrial ecosystems, unlike previous site-scale evaluations of the performance of C-N models (Zaehle et al 2014, Dybzinski et al 2019. Modeling the effect of N enrichment on Northern Hemisphere terrestrial C cycles in state-of-the-art C-N models should therefore be carefully considered. To more clearly evaluate C-N models' performance, it is essential to conduct in-depth comparative research of observations and model simulations that focuses on N enrichment effect on plant C turnover processes such as growth and mortality of leaf, stem and root. The variables reflecting the characteristics of C turnover across plant organs, therefore, are recommended to be included in the standard output variables of model intercomparison projects. In summary, our study provided new insights into the quantification of N enrichment impact on Northern Hemisphere plant biomass. N enrichment intensity was the main cause of the S AGB spatial pattern in Northern Hemisphere. N enrichment had a minor effect on Northern Hemisphere terrestrial BGB, indicating that the BGB:AGB ratio decreased as N increased, unlike S AGB . It is worth noting that the machine-learning models do not explain BGB response well (figures 2 and S7), likely due to the lower amount of observational data to constrain the BGB response to N enrichment. This may lead to considerable uncertainty in the ∆TB extrapolation in Northern Hemisphere. Thus, more effort on observing belowground C-N cycling is recommended in future N-enrichment experiments. Given that there is apparent spatial pattern of N limitation of plant growth across global natural terrestrial ecosystems (Du et al 2020), the effect of N enrichment on plant growth likely varies across different vegetation types. To accurately quantify the difference in N enrichment effect on AGB and BGB across various ecosystems, it will be useful to conduct more N-enrichment experiments in the ecosystems for which few data are available, such as shrublands, wetland and tundra (table S1). Particularly, long-term N-enrichment experiments are still insufficient (figure S1) but worth carrying out to explore the responses of plant C turnover processes such as mortality (Pregitzer et al 2008). Moreover, the components of atmospheric N deposition changed in recent years by increasing reduced N in the United States (Li et al 2016) and oxidized N in China (Yu et al 2019). The impact on biomass of such changes in the ammonium:nitrate ratio of N deposition remains to be studied in future field experiments.
All data that support the findings of this study are included within the article (and any supplementary files).