Aboveground biomass in Inner Mongolian temperate grasslands decreases under climate warming

Grassland aboveground biomass (AGB) is a critical component of the global carbon cycle and reflects ecosystem productivity. Although it is widely acknowledged that dynamics of grassland biomass are significantly regulated by climate 15 change, in situ evidence at large spatiotemporal scales is limited. Here, we combine biomass measurements from six long-term (> 30 years) experiments and data in existing literatures to explore the spatiotemporal changes in AGB in Inner Mongolian temperate grasslands. We show that, on average, annual AGB over the past four decades is 2,561 ka ha-1, 1,496 kg ha-1 and 835 kg ha-1, respectively, in meadow steppe, typical steppe and desert steppe in Inner Mongolia. The spatiotemporal changes of AGB are regulated by interactions of climatic attributes, edaphic properties, grazing intensity and grassland type. Using a 20 machine learning-based approach, we map annual AGB (from 1981 to 2100) across the Inner Mongolian grassland at the spatial resolution of 1 km. We find that on the regional scale, meadow steppe has the highest annual AGB, followed by typical and desert steppe. During 1981-2019, the average annual AGB generally exhibited a declining trend across all the three types of grassland. Under future climate warming, AGB in the study region could continue to decrease. On average, compared with the historical AGB (i.e., average of 1981-2019), the AGB at the end of this century (i.e., average of 2080-2100) would decrease 25 by 14% under RCP4.5 and 28% under RCP8.5, respectively. The decreases in AGB under warming show large disparities across different grassland types and future climate change scenarios. Our results demonstrate the accuracy of predictions in AGB using a machine learning-based approach driven by several readily obtainable environmental variables; and provide new data at large scale and fine resolution extrapolated from field measurements.


Introduction 30
Grassland occupies ~40% of the world land and is an essential component of global terrestrial ecosystems (Hufkens et al., 2016). Grassland provides plenty of ecosystem services such as suppling food to livestock therefore meat and milk to humans (Sattari et al., 2016) and accumulating carbon from atmosphere thus mitigating global warming (O'Mara, 2012). All of these functions are more or less directly dependent on grassland biomass, which has been recognized significantly influenced by environmental changes and anthropogenic activities (Hovenden et al., 2019). Thus, quantifying the dynamics of grassland 35 biomass and revealing underlying mechanisms, particularly at large extents of space and time, are of fundamental importance (Andresen et al., 2018).
Dynamics of grassland biomass are driven by complex interactions among a series of environmental attributes, among which climate is one of the most predominant drivers (De Boeck et al., 2008;Wang et al., 2020). The magnitudes and directions of climate change effects on AGB can vary across different local environments as well. In high altitude or latitude regions, climate 40 warming can avail to biomass accumulation because of the reduced constraints of low temperature on plant growth (Park et al., 2019;Gonsamo et al., 2018). For example, using a data-model integration approach, Zeng et al. (2019) found that annual plant biomass was significantly and positively correlated with temperature in Tibetan Plateau. In the arid and semi-arid systems, however, warming can have harmful effects on biomass formation because warming can aggravate the negative effects of water limitation on plant growth (Fan et al., 2009;Hu et al., 2007). 45 Mean annual climate attributes (e.g., temperature and precipitation) have widely been treated as the potential drivers on spatiotemporal changes in grassland biomass (Fan et al., 2009;Ma et al., 2008). However, growing evidences have demonstrated the importance of seasonality and intra-annual variability of climate, rather than the mean annual climate attributes, in regulating the grassland biomass dynamics (Grant et al., 2014;Godde et al., 2020). These effects of seasonality and intra-annual variability in climates have seldom been considered in studies focusing on large spatial extents, e.g., Inner 50 Mongolia, where China's largest temperate grassland locates in. It is noted that Inner Mongolian grassland accounts for more than half of China's northern temperate grassland area (Department of Animal Husbandry and Veterinary, 1996) and has the nation's largest grassland biomass carbon stock (Piao et al., 2004). The annual productivity of Inner Mongolian grassland, however, tends to vary in response to climate change (Bai et al., 2008). Since the start of 1980s, warming took place in many parts of Inner Mongolia (Wang et al., 2019). Under this warming, the spatiotemporal changes in grassland AGB, however, is 55 still unclear. Although efforts have been taken to quantify AGB variations at the regional scale, these quantifications using mainly remote-sensing approaches generally show large disparities (Guo et al., 2016;Ma et al., 2010a;Long et al., 2010).
Evidences from the datasets independent of remote sensing products may help to disentangle the mysterious spatiotemporal dynamics of AGB at the regional scale. Moreover, existing studies have seldom considered the possible co-regulating effects of soil properties (Jia et al., 2011;Bhandari and Zhang, 2019), grassland types and grazing intensity  Baquerizo, 2017) on AGB, which might lead to large uncertainties and biases in these estimations. In addition, little is known about the consequences of the likely future climate change on AGB dynamics across time and space, which is a substantial concern in the context of global warming.
In this study, we collate to date the most comprehensive dataset of in situ measurements on plant biomass and climatic records in Inner Mongolian grassland from six long-term (more than 30 year) experiments and those data in the region from existing literatures. We calibrate and validate a machine learning-based model for predicting the aboveground biomass in the study region, by treating tens of environmental covariates (climates, soils, grazing intensity, and grassland type) as predicting variables. Then, we map the annual aboveground biomass at a spatial resolution of 1 km in Inner Mongolian grassland over the periods of 1981-2019 (using historical climatic dataset) and 2020-2100 [using climate projections driven by two representative concentration pathways (e.g., RCP4.5 and RCP8.5)]. Based on the mapping products, finally, we analyse the 70 spatiotemporal variations in AGB across the study area over the study periods.

Datasets of grassland aboveground biomass
In this study, we acquired two datasets of in situ aboveground biomass (AGB) in Inner Mongolian grassland. First, we obtained the AGB at six long-term (i.e., more than 30 years) experimental sites across the study region (Fig. 1a, Data S1). These six 75 sites were established by Inner Mongolia Meteorological Bureau of China at early 1980s, measurements of AGB at each site has been carried out year by year since then. At each site, four fenced plots (i.e., four replicates) were set up to collect plant biomass data during plant growing seasons (e.g., from May to September). For each measurement replicate, the plants within a one square meter area were clipped and collected in a cloth bag. The samples were further air-dried to constant weights (weighted once every three days until the percent change in two consecutive weights are less than 2%). It is noted that plant 80 growth rate could peak at different periods across time and space. Following Scurlock et al. (2002), we determined the annual plant biomass as the largest observed monthly biomass during a year (normally at the end of August at Ergun and at the end of September at other three sites). Apart from these six long-term field experiments, we also retrieved a dataset regarding grassland AGB from Xu et al. (2018), who recently conducted a thorough literature synthesis and obtained a comprehensive dataset of plant biomass in the grasslands of northern China. For the dataset constructed by Xu et al. (2018), in this study, we 85 use only the observations conducted in Inner Mongolian grassland and with investigation time and coordinates clearly reported ( Fig. 1a). In general, the Inner Mongolian grassland AGB derived from these two datasets (i.e., long-term experiments and literature synthesis) are comparable (Fig. S1). In total, we obtained 511 individual measurements across 247 locations in Inner Mongolian temperate grasslands (Fig. 1a, Data S1).

Environmental covariates 90
Environmental covariates including climate, soil, grassland type and grazing intensity were retrieved for both AGB driver assessment and model fitting. For climatic covariates, we first obtained the daily climatic records of 120 climatic stations established in Inner Mongolia (Fig. 1b) from the National Meteorological Information Centre (NMIC) of China. The daily climatic attributes such as minimum, average and maximum temperature and precipitation were transformed into monthly time series data using the daily2monthly function in the R package hydroTSM. Based on these monthly data, we calculated 23 95 bioclimatic variables (Table S1) with an annual time step over the period of 1981-2019 by using the biovars function in the R https://doi.org/10.5194/acp-2020-1088 Preprint. Discussion started: 30 October 2020 c Author(s) 2020. CC BY 4.0 License. package dismo. By doing so, we aim to comprehensively consider the possible effects of seasonality, intra-and inter-annual variability of climates (Fick and Hijmans, 2017). By further applying a interpolation algorithm (Thornton et al., 1997) to these 23 bioclimatic variables at the 120 stations, we created the raster layers of the climatic attributes with a spatial resolution of 1 km year by year. For the edaphic covariates, we directly extracted 10 raster soil layers representing key soil physical and 100 chemical properties (Table S1) at a 1 km spatial resolution in the study region from ISRIC-WISE soil profile database (Batjes, 2016).
The grazing intensity in this study was represented by the quantity of three key animals (i.e., cattle, sheep and goats; Table S1) because they are the majority in Inner Mongolian grasslands (National Bureau of Statistics of China, 1981-2019). Here, we first derived the regional distribution data for cattle ( Fig region from Gilbert et al. (2018). Then, we obtained the yearly total amount of each livestock in the study region (Fig. S2 d) from National Bureau of Statistics of China (1981-2019). By assuming a similar spatial distribution of livestock over time, we generated raster layers of each of the three animals year by year over the past four decades using the above-mentioned two datasets. In addition, a spatial layer of grassland type (i.e., meadow steppe, typical steppe and desert steppe; Fig. 1a and Table   S1) at 1 km resolution was derived from the Vegetation Map of China (Zhang, 2007), the digital version of which is publicly 110 obtainable (http://data.casearth.cn/sdo/detail/5c19a5680600cf2a3c557b6b).

Machine learning models to predict grassland biomass
To predict grassland aboveground biomass (AGB) across the region, we generated a suite of machine learning-based predictive models for AGB, treating edaphic, climatic, grassland type and grazing intensity (Table S1) as candidate predictors. Here, data from the 511 measurements ( Fig. 1a and Data S1) were used to fit the models. For the spatial layers of soil properties and 115 grassland type, which were assumed to be constant over time, we retrieved the associated covariates using the geographical coordinates of the 511 measurements. For those variables varying over time (e.g., climatic variables and grazing intensities), we extracted the associated attributes using both the locations and investigating year of the 511 measurements. In fitting the models, AGB is treated as a dependent variable and the environmental covariates (Table S1) are treated as independent variables. Before fitting the models, we converted the categorical variables (i.e., grassland type) to dummy variables. Then, 120 the function findCorrelation in R package caret was used to exclude the environmental covariates with high multicollinearities.
The remaining attributes were further adopted in model training (80% stratified samples) and validation (the remaining 20% stratified samples). We used a 10-fold cross-validation to train a suite of machine learning models using three algorithms [i.e., random forest (RF), Cubist and support vector machines (SVM)], which are implemented in the R package caret. The amount of variance in AGB explained by each model was assessed by the coefficient of determination (R 2 ). The root mean square 125 error (RMSE) was also calculated to compare the model simulations and observations. Apart from the three individual models, we also derived an ensemble model by adopting a principal component analysis (PCA) approach based on the predictions of https://doi.org/10.5194/acp-2020-1088 Preprint. Discussion started: 30 October 2020 c Author(s) 2020. CC BY 4.0 License. the three algorithms. In brief, the smaller an individual model's RMSE is, the more the model's output contributes to final predictions.

Assessment of drivers on AGB 130
The machine learning models themselves provide assessments of the relative importance (RI) of each independent variable in predicting the dependent variable (e.g., grassland AGB in this study). In general, the greater the RI of a variable is, the larger its influence on AGB is. We also adopted the Mantel test (Mantel, 1967) to assess the relationship between similarity of different grassland types and the similarity of environmental covariates. Here, the standardized Mantel's r (ranges from 0 to 1) is used to represent the strength of this relationship (the higher the Mantel's r is, the stronger the correlation is) and the 135 associated significance is indicated by the p value determined from 999 randomization (Legendre and Fortin, 1989). Here, the R package vegan is used to perform the Mantel test.

Regional mapping and uncertainty analysis
Using the fitted machine learning-based ensemble model, we mapped AGB in Inner Mongolian grassland (at a spatial resolution of 1 km) on an annual time step in the history and future. In mapping the historical AGB (i.e., during 1981-2019), 140 the model is run using environmental covariates extracted from the regional data layers (see Environmental covariates).
Prediction uncertainty was quantified using a Monte Carlo analysis to develop the probability density functions (PDF) for each edaphic, climatic and livestock variable within the ranges of mean ±10%. The ensemble machine learning model was then run for 200 times in each grid with each of independent variables assigned from the PDF. The average and coefficient variation (CV, calculated as the standard deviation divided by the average) were then determined in each grid using the 200 model 145 outputs to represent the predicted AGB and the associated uncertainty, respectively.
For predictions of AGB in the future (i.e., 2020-2100), we include the climatic datasets projected by one typical CMIP5 global circulation model (GCM) to save computing resources. In this study, we use the projections output by CESM1-BGC, which was run by National Center for Atmospheric Research (NCAR). Here, we directly obtained the processed climatic products constructed by Karger et al. (2020), who recently generated a downscaled and bias-corrected temperature and precipitation 150 datasets. Specifically, these future climatic datasets were driven by two scenarios of representative concentration pathways (RCP4.5 and RCP8.5) at monthly step in this century. According to the model projections, mean annual temperature (MAT) under both RCPs will continue to increase in the following decades (Fig. S3). The extent of climate warming is generally higher under RCP8.5 than that under RCP4.5 (Fig. S3). After obtaining the future climate datasets, we also use the biovars function in R environment (see Environmental covariates) to calculate the 23 interested bioclimatic attributes (Table S1) for 155 both scenarios of RCPs year by year from 2020 to 2100. In projecting the future AGB dynamics using the ensemble machine learning model, we assume that the soil properties will not significantly change over time (i.e., the same soil inputs used in historical AGB predictions) and current grazing intensity will keep relatively stable (i.e., the average of the recent five years).
In addition, the uncertainty analysis for future AGB predictions were performed using the same approach as that adopted in mapping the historical AGB. All statistical analyses and graphical productions in this study were performed in R v3.6.3 (R 160 Development Core Team, 2020).
The fitted three individual machine learning algorithms (i.e., RF, Cubist and SVM) can explain overall 32%-48% of the variance in observed AGB (Fig. 3a, b and c). The ensemble model of the three algorithms can better simulate the observations than any of those individual models (Fig. 3). On average, 52% of the variance in the observations can be explained by the 170 ensemble model (Fig. 3d). Although the variable importance differed among the three algorithms, climatic and livestock variables seem to substantially affect the AGB dynamics (Fig. S4). After excluding the covariates with high multilinearities, the remaining 10 climatic attributes, 5 edaphic variables and three livestock predictors generally show small autocorrelations (Fig. 4). Mantel test suggests that, compared to the edaphic and livestock attributes, the climatic variables are in general stronger correlators of AGB in the three grassland types (Fig. 4). 175 The regional mapping results of grassland AGB during 1981-2019 show large spatial variations (Fig. 5a). On average, the regional AGB during the past four decades is 1,438 kg ha -1 , the corresponding lower and upper limits of the 95% CI is 479 kg ha -1 , and 2,284 kg ha -1 , respectively (Fig. 5a). Across grassland types, meadow steppe has the highest average AGB (2,194 Mg ha -1 ranging from 1,153 Mg ha -1 to 2,631 Mg ha -1 ), followed by typical steppe (1,552 Mg ha -1 ranging from 539 Mg ha -1 to 2,200 Mg ha -1 ) and desert steppe (893 Mg ha -1 ranging from 405 Mg ha -1 to 1,341 Mg ha -1 , Fig. 5a). Spatially, the average 180 coefficient of variation (CV) in the predictions is lowest in meadow steppe (10.5%), followed by desert steppe (14.6%) and typical steppe (21.8%, Fig. 5d). Over 1981-2019, the regional average AGB displayed a decreasing trend in the first three decades and then kept relatively stable in the recent decade (Fig. 6a). Among the three grassland types, the historical changes in AGB (Fig. 6b, c and d) are in general consistent with that of the total Inner Mongolian grassland AGB (Fig. 6a).
Our predicting results show that future AGB in general decreases under both scenarios of RCPs (i.e., RCP4.5 and RCP8.5,Fig. 185 5 and Table 1). Compared with the historical AGB (i.e., average AGB during 1981-2019, hereafter the same), on average, AGB at the end of this century (i.e., average of 2080-2100, hereafter the same) would decrease by 14% under RCP4.5 (Fig.   5b) and 28% under RCP8.5, respectively (Table 1). The decreases in AGB under climate change show large disparities across https://doi.org/10.5194/acp-2020-1088 Preprint. Discussion started: 30 October 2020 c Author(s) 2020. CC BY 4.0 License. different grassland types and climate change scenarios. Compared with the historical average AGB, AGB at the end of this century under RCP4.5 is estimated to decrease by a smaller extent (i.e.,10%) in meadow steep than those in typical (16%) and 190 desert steep (21% , Table 1). In general, AGB under RCP8.5 would reduce by larger extents compared with those under RCP4.5.
Under RCP8.5, the average AGB at the end of this century is estimated to experience a 24% (in meadow steep), 30% (in typical steep) and 25% (in desert steep) reduction, compared with the averages over 1981-2019 (Table 1). The magnitudes and spatial patterns of CV in the simulations under both RCP4.5 (Fig. 5e) and RCP8.5 (Fig. 5f) are comparable with those during the period of 1981-2019 (Fig. 5d). 195

Discussion
Our results, based on AGB observations derived from six long-term field experiments and literature synthesis, indicate the large spatial disparities in aboveground biomass across different grassland types (Fig. 2). This gradient spatial pattern in AGB, i.e., largest in meadow steep followed by that in typical and desert steep, is comparable with Ma et al. (2008), who carried out a comprehensive field measurements and investigated 113 locations in Inner Mongolian temperate grassland during  2005. On the regional scale, we mapped grassland AGB at high spatial and temporal resolutions, which shows that AGB generally decreases from north-eastern to south-western areas in the study region (Fig. 5a). Such a spatial pattern is also consistent with the maps generated from remote sensing derivations (Fig. S5). This demonstrates the accuracy of our machine learning model's predictions. It should be noted that existing mapping products of grassland AGB use mainly remote sensing approaches requiring inputs from satellite-based datasets (Guo et al., 2016;Ma et al., 2010a;Jiao et al., 2019). Our fitted 205 machine learning model, however, uses only several readily obtainable environmental covariates ( Fig. 4 and Table S1). Our results demonstrate the ability of machine learning approach to effectively extrapolate grassland AGB to much larger spatiotemporal extents (e.g., Fig. 5 and 6).
Our simulation results show that, under the climate warming over the past four decades (Fig. S3), the average AGB generally experienced a declining trend across all the three grassland types in Inner Mongolia (Fig. 6). This demonstrates the possible 210 negative effect of temperature rising on AGB that has been widely reported (De Boeck et al., 2008;Wang et al., 2020), particularly in the arid and semi-arid ecosystems (Ma et al., 2010b). This harmful influence of warming on AGB is explainable.
For example, in a system restrained by water availability (e.g., temperate grassland), warming can not only inhibit plant photosynthesis (Xu and Zhou, 2005) but also enhance evaporation and further intensify water stress (De Boeck et al., 2006) thereby decreasing grassland biomass. Precipitation has generally been recognized to have positive effects on AGB in the 215 temperate grassland (Hovenden et al., 2019;Ma et al., 2010a), which supports our findings in this study. For example, the simulated average AGB is relatively higher in the years with higher MAP (e.g., 1998 and 2012) than those in other years (Fig.   6a). The importance of precipitation on AGB can be more reflected by the spatial patterns of these two attributes, e.g., AGB is much lower in the more arid regions (Fig. 5a) where soils are suffering severer water deficiencies. Apart from climatic https://doi.org/10.5194/acp-2020-1088 Preprint. Discussion started: 30 October 2020 c Author(s) 2020. CC BY 4.0 License. factors, our results also demonstrate the regulating effects of soil conditions and livestock on the dynamics of grassland AGB 220 ( Fig. 4 and S4). This is consistent with several findings highlighting the importance of soil physical and chemical characteristics (Yang et al., 2009;Griffiths et al., 2012) and grazing intensity (Eldridge and Delgado-Baquerizo, 2017) in controlling grassland biomass changes.
We notice that our model predictions show larger interannual variations in AGB (Fig. 6a) than those in the estimations based on remote sensing approaches (Fig. S5). In fact, the remote sensing derived AGB has also been bias-corrected by the field field experiments (Fig. 1a), which can potentially reduce the possible biases in model predictions. In addition, we find that the overall decrease in Inner Mongolian grassland biomass are contributed greater by the decline during the first three decades and the declining trend in AGB seems to be alleviated in the recent decade (Fig. 6). This could be related to the overall slowing climate warming over the recent decade (Fig. S3). In the future, a faster warming (e.g., RCP8.5) climate will lead to a larger reduction in grassland AGB (Fig. 5b and c). It is noteworthy that the accuracy of our predictions on future grassland AGB 235 relies substantially on the robustness of future climate change projections simulated by the GCMs (e.g., CESM1-BGC).
However, although CESM1-BGC (like all other CMIP5 models) can reasonably well simulate changes in temperature, it may not well predict precipitation, particularly for Eastern China where is strongly affected by large-scale atmospheric circulations (Huang et al., 2013). Consequently, it should be cautious in interpreting the likely decrease of AGB under future temperature rising (Fig. 6), because the large uncertainties in projected precipitation may lead to biased predictions in AGB. 240

Conclusions
Our results demonstrate that the aboveground biomass in Inner Mongolian grasslands shows large spatial and temporal variations during the past four decades, which is driven by a series of environmental covariates. Particularly, current climate change characterized mainly by warming has negative effects on AGB across all types of grassland. In addition, our results demonstrate that adopting a machine learning model approach with only a few readily obtainable environmental predictors 245 can accurately capture AGB dynamics, which enables extrapolations of AGB across larger spatiotemporal extents. Moreover, our study provides new data on annual AGB in the study region at fine spatial (1km) and temporal (yearly) resolutions for both historical (1981-2019) and future (2020-2100) periods under different climate change scenarios.
Data availability. The data that support the findings of this study (Data S1) are openly available at: 10.6084/m9.figshare.13108430.