Vegetation biomass change in China in the 20th century: an assessment based on a combination of multi-model simulations and field observations

Vegetation biomass is a key and active component of the carbon cycle. Though China’s vegetation biomass in recent decades has been widely investigated, only two studies have quantitatively assessed its century-scale changes so far and reported totally opposite trends. This study provided the first multi-model estimates of China’s vegetation biomass change for the 20th century and its responses to historical changes in environmental and anthropogenic factors, based on simulations evaluated with the field observations from 3757 inventory plots in China and bias-corrected using machine learning (Gaussian process regression). A significant decline in vegetation biomass over the 20th century was shown by bias-corrected simulations from the six Dynamic Global Vegetation models (DGVMs) with trends ranging from −32.48 to −11.10 Tg C yr–1 and a mean trend of −17.74 Tg C yr–1. Land use and land cover change (LULCC) was primarily responsible for the simulated downward trend (−50.71 to −24.28 Tg C yr–1), while increasing atmospheric CO2 concentration lead to increased vegetation biomass (+9.27 to + 13.37 Tg C yr–1). Climate change had limited impacts on the long-term trend (−3.75 to + 5.06 Tg C yr–1). This study highlights the importance of LULCC for historical reconstruction and future projection of vegetation biomass over China. It also suggests that the incorrect change in China’s forest area for 1980–2000 in the LULCC dataset used as model input data of many existing and ongoing model intercomparison projects (MIPs) has likely led to inaccurate estimations of historical vegetation biomass changes in China.


Introduction
China is the most populous country and the largest carbon emitter in the world, and has undergone distinct land-use and climate changes during the 20th century (Fang et al 2018). Assessment of vegetation biomass and its change over China is essential for understanding the regional, and even the global, carbon cycle (Kondo et al 2013, Fang et al 2018.
Many studies have estimated China's vegetation biomass in recent decades using various approaches. Estimates of present-day vegetation biomass in China range from 13.34 to 35.23 Pg C based on numerical models (Li et al 2004, Ni 2001, and references therein), and from 13.09 to 34.68 Pg C based on field and satellite observations (Carvalhais et al 2014, Li et al 2015, Xu et al 2018, Fang et al 2018and references therein). Fang et al (2018) and Tang et al (2018) estimated that forests, grasslands, shrublands, and croplands contributed to 80%, 10%, 5%, and 4% of current vegetation biomass, respectively, using extrapolation from site-level field observations based on a land cover map estimate. Xu et al (2018) investigated the spatial changes of 2004-2014 vegetation biomass by collecting data based on literature review, and concluded that mean annual precipitation is an important driver of the spatial changes. Several studies reported that vegetation biomass increased after ∼1980 probably as a result of forest expansion and regrowth (Fang et al 2001, 2007, Tang et al 2018. However, until now, only Wang et al (2007) and Mao et al (2009) have investigated the long-term change of China's vegetation biomass over the 20th century. Both studies were based on a single model, specifically, the Integrated Terrestrial Ecosystem Cbudget model (InTEC) in Wang et al (2007) and a modified version of the Sheffield Dynamic Global Vegetation Model (SDGVM) in Mao et al (2009). The two studies reported opposite trends: China's vegetation biomass for the 20th century was decreased by 2.15 Pg C from Wang et al (2007) while increased by 2.07 Pg C from Mao et al (2009).
Historical changes in vegetation biomass are regulated by climate, atmospheric CO 2 concentration, and land use and land cover change (LULCC) (Kim et al 2017, Fang et al 2018, Xin et al 2018. Dynamic Global Vegetation Models (DGVMs) simulate the response of global terrestrial ecosystems to these drivers, and integrate biogeochemistry, biogeophysics, and vegetation dynamics at the land surface within a physically and chemically consistent modelling framework (e.g. Foley et al 1996;Sitch et al 2003). They can be used to isolate how each driver affects vegetation biomass, providing that the models produce (or can be induced to produce via bias correction) realistic simulations of the modern state (Piao et al 2011, Ichii et al 2013, Chen et al 2016. However, up to now, there has been no attempt to evaluate vegetation biomass simulated by DGVMs or correct their simulation bias before applying them to investigate the historical changes and drivers of global and regional vegetation biomass. The latter is important given that DGVMs aim to model the global simulations with low spatial resolution and generally do not conduct regional optimization. Furthermore, although Mao et al (2009) found that atmospheric CO 2 concentration change was more important than changes in temperature and precipitation in explaining vegetation biomass change of the 20th century, there has been no study to investigate the potential impact of other drivers on the historical change.
In this study, we used simulations by six stateof-the-art, fire-enabled DGVMs which were run following the same protocol and the same forcing datasets. We evaluated these simulations using sitespecific biomass data from 3757 inventory plots in China. Biases in each model simulation were statistically corrected using machine learning (Gaussian process regression, GPR). The bias-corrected simulations provided the first multi-model estimates of vegetation biomass for the 20th century. We then used three sensitivity experiments in which potential drivers of 20th century vegetation changes, specifically, climate, CO 2 concentration and land-use, were constant, respectively, throughout the historical period to analyze the causes of the simulated changes in vegetation biomass.  (table 1). All the DGVMs are not stochastic except for LPJ-GUESS. LPJ-GUESS includes stochastic representation of individual establishment, mortality, and disturbance impacts. This does not affect the interpretability of the results because the simulations used here feature sufficient (50 for LPJ-GS and 25 for LPJ-GSB) replications so that the effect of the stochastic variability is small and the effect of drivers is dominant.

Model simulations
All the DGVM simulations started from a spinup simulation for 1700 conditions to get the initial fields for subsequent 1701-2013 simulations (figure 1, see Rabin et al 2017 for details). For the spin-up runs, models were recycled until carbon values in the slowest soil carbon pool varied by less than 1%   between consecutive 50 year periods for every grid cell. The atmospheric reanalysis data were from CRU-NCEP v5.3.2. Annual global atmospheric CO 2 concentration was derived from ice core and NOAA monitoring station data ( In the no climate change simulation, the climate throughout the run was the recycled time-varying 1901-1920 atmospheric forcing, but all other inputs were allowed to vary as in the simulation. In the constant CO 2 simulation, the atmospheric CO 2 concentration was held constant at 277.33 ppm for the whole of the simulation. In the no LULCC simulation, land cover in 1700 was used throughout the simulation. All model outputs were converted to 1º spatial resolution, based on bilinear interpolation from low to high resolution for CLM4.5 (∼1.9 • latitude × 2.5 • longitude), JSBACH (1.875 • ), and JULES (∼1.2 • latitude × 1.9 • longitude), and on area-weighted averaging for those models with an original spatial resolution of 0.5º.
Not all of the FireMIP models were used in this study. LPJ-GUESS-GlobFIRM and LPJ-GUESS-SIMFIRE-BLAZE had the same vegetation model but used different fire models. Evaluation by Andela et al (2017) and Li et al (2019) showed that the fire model GlobFIRM simulated the global burned area and fire carbon emissions poorly. MC2 was developed for regional applications and applied globally without adequate calibration (  some bugs were identified in temporal variability of carbon fluxes and pools in its new version. Therefore, the simulations of the three models (LPJ-GUESS-GlobFIRM, MC2 and CLASS-CTEM) were excluded from our analyses.

Field observations
The field observations of vegetation biomass used here were collated and derived from Chinese forest survey inventories of Luo (1996) and Pan et al (2004) and from literature review (supplement B.xlsx and supplement B reference.docx). The vegetation biomass is derived as the sum of leaf, stem and branches, and root biomass or the sum of above-ground and below-ground biomass. The dataset from Luo (1996) and Pan et al (2004) focused on forest plots, and provided information on location, forest type, and leaf, branch, stem, and root biomass for trees for all plots. It included 1266 forest plots, sampling tropical and monsoon forests, subtropical evergreen broadleaf and coniferous forests, temperate deciduous broadleaf forests, boreal evergreen/deciduous coniferous forests where the measurements were made on forests aged from 3 to 350 years during 1956-1993. Most of the organ biomass (leaf, branch and stem, root) was obtained by the biomass allometric growth model method. Missing values of tree roots were calculated from the ratio of above-ground to below-ground biomass of the same species in nearby stands (Luo 1996).
Furthermore, we compiled an additional dataset of above-/below-ground vegetation biomass for 1112 forest plots and 1379 grassland plots measured during 2000-2013 from published literature. We selected plots with observations of both above-and belowground biomass and a clear description of the functional types or species. Observations in the literature were mainly obtained by the clear-cutting method, average standard tree method, and/or biomass allometric growth method. A coefficient of 0.45 was used to convert vegetation biomass density to C density (kg C m −2 ) (e.g. Xu et al 2018).

Bias correction for vegetation biomass
To derive bias-corrected outputs from each of the DGVMs, we built statistical models between the default (uncorrected) vegetation biomass simulations and observations of vegetation biomass (total vegetation biomass) as follows. First, according to the vegetation-type or species records in the field observations and the plant functional type (PFT) information in the model outputs, we sorted the observations and simulations into evergreen broadleaf, deciduous broad-leaf, evergreen needle-leaf, and deciduous needle-leaf tree PFTs and grass (not separating C3 and C4) PFT. Second, simulations were temporally averaged to match the representative periods of observations. Third, we built statistical models based on simulation-observation pairs that belonged to the same PFT over all grid cells.
Three methods for bias correction were tested: first order linear regression: cubic regression: and machine learning (Gaussian process regression, GPR). Here, a i,j,n (n = 0, .., 3) was regression coefficient for the ith functional type of the jth model, B i,j,obs and B i,j,uncorrected were the spatial arrays of field observations and the default (uncorrected) model outputs, respectively. For the machine learning GPR, we used the GPR package of Matlab version 2018b. We carried out batch training based on the GPR model to optimize the parameters of kernel function (ardsquared exponential function). Hyperparameters were then determined by automatic hyperparameter optimization based on minimizing five-fold crossvalidation loss. We adopted the commonly-used leave-one-out cross-validation framework (LOOCV, Wilks 1995) to build and evaluate the statistical models based on the three bias-correction methods for preventing over-fitting. In LOOCV, one observation-simulation pair was left out as test set, and the statistical model was built by using remaining observation-simulation pairs (n − 1) as training set. Then, the bias-corrected simulation value was predicted with the statistical model and the independent left-out simulation (uncorrected) value. The process was repeated n times to obtain n bias-corrected simulation values. The correlation coefficient (R) and root mean square error (RMSE) between simulations and field observations were used as skill scores (skill is higher with higher R and lower RMSE). We selected the method with highest cross-validation skill as the correction method for all models and functional types which was denoted as correction function f. Analyses of cross-validation skill (table S1 (available online at stacks.iop.org/ERL/15/094026/mmedia)) showed that machine learning GPR outperformed the first order linear and cubic regressions for all functional types and models except deciduous forest in JULES.
We derived the bias-corrected biomass simulations (B i,j,k,corrected ) for functional type i of model j in year k during 1950-2013 as: Then, the total vegetation biomass ratio between with and without bias-correction averaged during 1950-2000 was used as bias-correct coefficient to modify the mean total vegetation biomass in 2001-2010, as well as trend in vegetation biomass change in the 20th century. This method assumes that the statistical relationship between the multi-year average of DGVM simulations and field observations is constant.
Before bias correction, all the DGVMs showed a systematic bias in simulated vegetation biomass for tree, and grass PFTs compared to observations (figures S1-2). Specifically, CLM4.5, JSBACH, and JULES overestimated the vegetation biomass for evergreen tree PFTs, while ORCHIDEE underestimated it. LPJ-GS and ORCHIDEE overestimated the vegetation biomass for deciduous tree PFTs (figure S1). All models tended to underestimate grass vegetation biomass ( figure S2). Moreover, the spatial heterogeneity of vegetation biomass was considerably underestimated by JULES and LPJ-GSB for tree PFTs (figure S1) and by all models except for CLM4.5 for grass PFTs (figure S2).
The bias-corrected results showed a significant improvement compared to the default simulations (figure 2). The correlation coefficients between corrected simulations and observations were much higher than uncorrected ones, and the RMSE was also lower in the bias-corrected simulations. All correlation coefficients for bias-corrected simulations passed the Student's t-test at the 0.001 significance level, indicating that the bias-corrected simulations can reproduce China's vegetation biomass reasonably well. In addition, the bias-corrected simulations provided an estimate of China's vegetation biomass averaged for the period 2001-2010 ranging from 10.57 ± 0.05 to 15.47 ± 0.13 Pg C with a MMEM of 12.48 ± 0.06 Pg C (table 3)

Historical changes over the 20th century
All the six DGVMs exhibited an overall decline in vegetation biomass over the 20th century for China   for the 20th century. The different change of ORCH-IDEE after the 1960 s was caused by the stronger response to rising [CO 2 ] (figure 4(b)) partly due to no modeling of nitrogen limitation (table 1) and the weaker response to LULCC since the 1980 s (figure 4(c)) likely due to weaker decline in forest area (figure 6(a)) than other models. The MMEM of all DGVMs also exhibited a significant downward trend (-17.74 Tg C yr −1 ).
All the DGVMs and MMEM agreed that rising [CO 2 ] between 1901 and 2000 tended to increase vegetation biomass. The trend ranged from + 9.27 Tg C yr −1 to + 13.37 Tg C yr −1 , with a trend of around + 10.04 Tg C yr −1 ( figure 4(b)). Rising [CO 2 ] increased vegetation biomass over most regions of China (figure S4) probably because the models were constructed to respond to changes in [CO 2 ] by changing light absorption and light-use efficiency (Norby et al 2005, Mao et al 2009 and also changing forest water-use efficiency (Keenan et al 2013).
The model simulations for all the DGVMs and MMEM showed that LULCC decreased vegetation biomass over the 20th century ( figure 4(c)). The effect of LULCC ranged from a strong decrease of −50.71 Tg C yr −1 for LPJ-GS to a moderate decrease of −24.28 Tg C yr −1 for JULES with an estimate of −39.47 Tg C yr −1 for the MMEM. LULCC impacts were the dominant driver for long-term changes in vegetation biomass over China for MMEM and all models ( figure 4). Spatially, the largest impacts of LULCC on vegetation biomass were in the northeastern and southwestern forested regions ( figure S5).
The impacts of climate change on the long-term trends of vegetation biomass were generally limited (figure 4(d)), with trends ranging from −3.75 Tg C yr −1 (CLM4.5) to + 5.06 Tg C yr −1 (JSBACH). There was no significant trend (at the 0.05 level) in the MMEM (-0.02 Tg C yr −1 ). This was probably because the impacts of climate change over the 20th century are generally weak and have large spatial heterogeneity (figure S6).
The DGVM simulations without bias-correction exhibited similar trends in vegetation biomass and similar responses to the various drivers, but the downward trend and responses were stronger than those in the bias-corrected simulations (Figs. S7-11), partly because the magnitude of vegetation biomass was overestimated in the uncorrected simulations (table 1, figure S3). For example, the decline in the vegetation biomass over the 20th century in the uncorrected simulations ranged from −71.70 to −14.36 Tg C yr −1 (-31.81 Tg C yr −1 for MMEM) (figure S7), compared to trends of −32.48 to −11.10 Tg C yr −1 for bias-corrected simulations and −17.74 Tg C yr −1 for the MMEM (figure 4(a)).
The bias-corrected MMEM showed a significant decline in vegetation biomass over most of China for the 20th century, especially in the southwest ( figure  5(a)). Only limited areas in northern China showed a slight increase in vegetation biomass. LULCC was the main reason for changes in vegetation biomass over 72.42% of China's land area ( figure 5(b)), and largely driven by a decline in forest area. From 1901 to 2000, forest area in the DGVMs declined from −1.38% in JULES to −9.43% in ORCHIDEE ( figure  5(a)). The areas dominated by the impact of rising [CO 2 ] were small, while climate change was the dominant driver for the increase in vegetation biomass only in part of northwestern China ( figure 5(b)). Again, these conclusions were still robust without bias-correction, where LULCC was dominant over approximately 74.94% of China (figure S11).

Discussion and conclusions
Because long-term vegetation biomass observations over large regions are unavailable, DGVMs are a good choice for quantifying the trends in regional and global vegetation biomass and isolating the impacts of each driver on the scale of a century (Piao et al 2011). Before trend and attribution analyses, evaluation should be done and statistical or dynamical post-processing are often needed for the regional use of DGVM simulations because many parameters are global constants and some schemes are not optimized to fit the region of focus (Oleson et al 2013).
Our evaluation results showed that the DGVMs often overestimated the amount of vegetation biomass in China and underestimated the spatial heterogeneity of vegetation biomass for tree and grass PFTs. The overestimation in China's vegetation biomass is partly due to underestimated turnover rate and overestimated photosynthesis for the DGVMs (table S3, supplement C) as well as overestimated forest area by 12.2% for LPJ-GS. The underestimation in spatial heterogeneity of vegetation biomass for a PFT is likely because the DGVMs generally do not consider the diversity in plant traits and individual growth within a PFT (Scheiter et al 2013).
To reduce the bias in DGVM simulations, machine learning GPR was used as bias-correction method. The bias-corrected DGVM simulations were much more skillful, and could simulate China's vegetation biomass reasonably well and reduced the uncertainty in DGVM simulations (smaller inter-model spread). This shows how statistical bias-correction can be used as the post-processing of model outputs to improve simulation skill and increase the reliability of subsequent quantification of temporal and spatial changes in terrestrial ecosystems. This method of bias-correction could also be applied to simulations of future changes in the vegetation biomass of China. Providing there are sufficient observations, it would be possible to use GPR bias-correction in other regions or for other simulated variables.
Our results also showed that China's vegetation biomass decreased over the 20th century with significant trend of −17.74 (−32.48 to −11.10) Tg C yr -1 . This is consistent with the results from Wang et al (2007)  between different drivers to be identified. Four additional sensitivity experiments, including no change in all drivers, no CO 2 and land cover change, no climate and CO 2 change, no CO 2 and climate change, would facilitate a more comprehensive attribution analysis in the next phase of FireMIP (i.e. total n 2 −1 simulations, n is the number of drivers).
However, our bias correction still has limitations. First, the machine learning GPR is a statistical method and has the common problem that assumes the statistical relationship between dependent variables and independent variables (multi-year average of DGVM simulations and field observations in this study) to be constant. Second, although the field data that we collected and used had been quality controlled, they still included uncertainty and possibly errors due to sampling design, measurement techniques, and observer bias. Such uncertainty and errors were not taken into account in our bias-correction, and may affect our estimation of trends in China's vegetation biomass. In addition, the field data used in this study represented multi-year averages and did not provide information on temporal change, so the coefficients of bias-correction models were based on multi-year averages. That is, this study just corrected the multiyear amounts and spatial pattern, which would affect our trend estimates. Long-term field measures which could provide information on temporal change are needed to improve bias-correction further and our quantitative understanding of large-scale and longterm changes in vegetation carbon pool and their drivers.
Our estimated decline in vegetation biomass after ∼1980 is likely incorrect. It is not consistent with earlier studies based on Chinese forest inventory data (Fang et al 2001(Fang et al , 2007 which reported an increase in vegetation biomass from 1981 to 2000. This opposite trend is mainly because the opposite change in China's forest area used by these studies. The Chinese forest inventory data, also reported in Food and Agriculture Organization of the United Nations (FAO) national-level forest resource assessment, shows a sharp rise in China's forest area due to reforestation policy since ∼1980 s (figure 6(b)) (Li and Li 1996, Ma et al 1997, State Administration of Forestry and Grassland 2017. However, in this study, the change in China's forest area in the DGVMs, except for ORCHIDEE whose forest area change was constrained by Houghton et al (2012), was characterized by a decrease after 1980 (figure 6(a)). The FireMIP models used the cropland and/or pasture area changes from Hurtt et al (2011), and prescribed or modeled the changes in tree, grass, and or shrub PFTs areas which shared the remaining grid area after assigning areas for water/ice/urban (constant) and cropland/pasture (prescribed), an approach widely used in international MIPs, e.g. TRENDY (Sitch et al 2015),  (Lawrence et al 2016), that is the LULCC input data for many ongoing MIPs, e.g., CMIP6 as a basis for the next Intergovernmental Panel on Climate Change Assessment Report (IPCC AR6) (Eyring et al 2016). Given the importance of forest area change in simulated vegetation biomass change verified in the present study, the use of incorrect LULCC data will result in inaccurate simulation of the carbon sink and source in China for the last two decades of the 20th century in these MIPs. Therefore, the Chinese forest inventory data, and satellite-based land cover change data from Song et al (2018) are recommended to LUH group for improving the accuracy of LULCC and thus simulations of various MIPs and IPCC reports.