Great uncertainties in modeling grazing impact on carbon sequestration: a multi-model inter-comparison in temperate Eurasian Steppe

The impact of grazing activity on terrestrial carbon (C) sequestration has been noticed and studied worldwide. Recent efforts have been made to incorporate the disturbance into process-based land models. However, the performance of grazing models has not been well investigated at large scales. In this study, we performed a spatially explicit model uncertainty assessment in the world’s largest pasture ecosystem, the temperate Eurasian Steppe. Five grazing models were explicitly incorporated into a single terrestrial biogeochemical model to simulate regional C consumption from grazing activity (Cgraze). First, we summarized the underlying mechanisms and explicitly compared the general functions used to describe the processes in different models. Then, the models (five models with 12 simulations) were run in parallel using the same forcing data and livestock distribution map in 2006. Results indicated that the modeled regional Cgraze varied from 0.1–16.1 gC m−2 for the year. The corresponding ratios of Cgraze to aboveground net primary productivity ANPP and net primary productivity (NPP) ranged from 0.08%–24.6% and 0.028%–11.2%, respectively. Parameter sensitivity was further analyzed. Model outputs are highly sensitive to the intake rate (i.e. feeding rate of livestock per day), half maximum intake rate, and initial livestock weight. Our results indicate that great uncertainty exists in simulating Cgraze. We ascribed the major uncertainty to the different process description and poor parameterization. This study calls for more efforts to the comprehensive synthesis of usable dataset, the foundation of a standard observation system and the observe-based inter-comparison to evaluate models, which would facilitate more accurate assessment of C sequestration by pasture ecosystems and lead to better representation in earth system models.


Introduction
Grazing activity is a disturbance to terrestrial systems worldwide. Livestock-induced carbon (C) emissions account for 5.6-7.5 Gt CO 2 -eq or 11%-18% of global human-made GHG emissions based on different analyses (Gerber et al 2013, Herrero et al 2016, Tubiello et al 2013. Grazing activity alone emits an estimated 1.32 Gt CO 2 -eq, which is about 20% of the total emissions from livestock related production and supply (Gerber et al 2013). As a direct production removal process, the resource use efficiency of extensive grazing systems is the lowest among the livestock sector (de Vries and de Boer 2010, Food and Organization 2010, Herrero et al 2013. The current estimates of grazing-related C emissions largely rely on inventory datasets, mainly from the Food and Agriculture Organization of the United Nations (FAO) dataset as reported by Garnett et al (2017). A general approach in estimating the C emission from grazing is to calculate the difference between the feed demand from livestock and the market feed supply (so called 'grazing gap approach') (Fetzel et al 2017). For example, Krausmann et al (2008) estimated the feed demand according to the correlation between the daily feed demand and livestock products (e.g. milk yield and carcass weight). A similar approach was applied in calculating the contribution of grazing C consumption to human appropriated net primary productivity from grazing land (Haberl et al 2007, Krausmann et al 2013. Herrero et al (2013) estimated the feed demand from livestock using a more comprehensive dataset of feed-composition from multiple sources including FAO and the international livestock research institute.
In order to understand how grazing influences the terrestrial C cycling in a process based manner, grazing models were developed and coupled with vegetation models. Earlier attempts were more conceptual to test assumptions and underlying mechanisms, and to characterize specific processes (Bailey et al 1996, Fryxell 1991, Hamilton et al 1973, Illius and O'Connor 2000, Milchunas et al 1988. Other model developments have been successfully tested and applied in specific pasture systems, including the defoliate rate model (DRM) in Argentina andMongolia (Chen et al 2006, No'am et al 1992), the Pasture Simulation model (PASIM) in western Europe (Riedo et al 2000), the Illius and O'Connor model (IOM) in southern Africa (Illius and O'Connor 2000), and the Shiyomi model (SM) in Inner Mongolia, China (Shiyomi et al 2011). Developed in 1994, the SAVANNA model might be the most comprehensive model to simulate the pasture systems at landscape and regional scales (Coughenour 1993). But its application on large scales is also limited by the detailed inputs and parameters required to force the model.
Increasing recognition of the importance of simulating land processes in predicting global C budget has stimulated more recent efforts to focus on incorporating the grazing factor into the global vegetation models (GVMs, typically part of the land model component in Earth system models). At the current stage, grazing activity has been explicitly represented as a C removal process in several GVMs (Erb et al 2017). Further implementations have been planned to incorporate more detail representation of grazing practices in more GVMs and participate in the future inter-model comparison projects (Pongratz et al 2017). Instead of testing specific processes or mechanisms, the more recent development of grazing models try to quantitatively evaluate how grazing activity influences the terrestrial energy and matter cycling and contributes to the greenhouse gas emissions. Examples are the coupling of PASIM into ORCHIDEE (Chang et al 2013, Vuichard et al 2007b, IOM into LPJ-GUESS (Pachzelt et al 2013), SM into BEPS (Chen et al 2017) and DRM into Biome-BGC (Luo et al 2012). These models have all been applied in regional or global assessments of grazing impact on terrestrial C sequestration (Chang et al 2016a, Chang et al 2016b, Christensen et al 2004, Han et al 2016, Pachzelt et al 2015, Vuichard et al 2007a.
With the increasing model developments and applications at large scales, different model structures and parameterizations are applied. Although the models are validated using various observation approaches (Chang et al 2013, Shiyomi et al 2011, their performances have never been evaluated over a unified region. Uncertainties in estimated influence of grazing on regional C sequestration remain unclear. In this study, we present a grazing model intercomparison on the largest contiguous pasture systems on earth, the temperate Eurasian Steppe (TES). We (1) explicitly incorporated the five grazing models (i.e. IOM, PASIM, SM, DRM. and the SAVANNA model) into a single terrestrial biogeochemical model, boreal ecosystem productivity simulator (BEPS), to assess the uncertainties in predicting the grazing impact on regional C sequestration in TES; and (2) analyzed the parameter sensitivity in each model. The aim of this study is to quantify the magnitude of uncertainty in modeling the grazing impact on C sequestration and try to identify the potential uncertainty sources in predicting C consumption from grazing activity.  (Chen et al 1999, Liu et al 1997. This model was continuously improved and applicable for simulating C, water, energy cycling in various ecosystems (Chen et al 2007, Ju et al 2006, Liu et al 2003. It has been validated and applied in multi-scale studies (Zhang et al 2012, Zhang Fangmin et al 2014. This model was recently revised to improve its simulation ability on grassland ecosystems in the TES (Chen et al 2017, Chen et al 2016).

Summary of the grazing model approaches
Grazing activity influences terrestrial C cycling mainly through three processes: direct intake, excretion, and trampling (Liu et al 2015b, Vuichard et al 2007b. Here we summarized the major approaches and functions used in the current grazing models.

Direct defoliation
Direct intake from livestock is the primary impact of grazing on grassland ecosystems. It directly removes biomass from plants and standing litters, which reduces carbon inputs to litter and soil carbon pools. Approaches to model direct defoliation are categorized into three types, described below.
DRM quantifies the grazing intake by assuming a linear correlation with the available biomass, which is basically derived from the defoliation rate equation (No'am et al 1992). The defoliation rate by grazing from a single livestock unit for a single time step (dB graze /dt) is calculated as: where r in is the daily intake rate or efficiency per unit livestock (ha per livestock unit per time step) and B ava is the available biomass for the time step (kg ha −1 , dry matter). An upper limit of forage intake is set in DRM. SM quantifies the grazing intake mainly from the perspective of grazing demand associated with livestock weight (Shiyomi et al 2011): where W is the livestock weight at time step t, r intake is the intake rate and q is the fraction of composition i (i = live biomass or dead biomass). In SM, the forage availability is used as a limiting scalar to confine the livestock intake (table 1). The third approach is based on the responses of intake rate to the internal (digestion ability) and external (forage availability) factors (Fryxell 1991). They are described by the michaelis-menton (M-M) formulation: where C 1 is the intake rate limited by the forage availability (g per day per livestock unit), C 2 is the intake rate limited by the forage quality (g per day per livestock unit), r maxand e max are the intake mass with maximum consumption (g per day per livestock unit) and unrestricted intake mass with maximum digestible energy content (g per day per livestock unit), respectively. Therefore, the real intake rate is predicted as the minimum value of the two constraints: dB graze dt = min( 1 , 2 ).
Simplified from Fryxell's modeling approach, the three models, PASIM, IOM and the SAVANNA model, only consider the external limitation (i.e. equation 3). In some models, other related factors were used to adjust the real intake rate. So the general equation for the current model could be written as: where f represents other limiting factors such as the animal sex composition in PASIM, age level in IOM, and weight condition in the SAVANNA model. As data are not available for the factors involved in f, many of them cannot be explicitly considered in large scale simulations. Forage digestibility is explicitly considered to influence the intake rate in DRM, IOM and the SAVANNA model. In DRM, forage digestibility is set differently for living and dead biomass. In IOM and the SAVANNA model, forage digestibility of living biomass is estimated as a function of soil water content and a function of leaf nitrogen content, respectively. Forage digestibility of the dead biomass are set as fixed values in the two models.

Excretion
Excretion (urine and feces) from livestock returns C to the soil, and facilitates re-distribution of soil nutrients (Stumpp et al 2005), but also negatively affects aboveground vegetation due to the urine scorching effect (Guthery and Bingham 1996). The C return from excretion is represented as a set percentage from grazing intake in SM and DRM. It is calculated as the difference between animal intake and the sum of other C outputs in PASIM. The scorching effect is generally estimated together with the trampling effect (see below).
In PASIM, DRM, and the SAVANNA model, the CH 4 and N 2 O emissions from livestock are explicitly estimated as a fraction of C and N intakes.

Trampling and scorching effect
Trampling mainly causes changes in soil physical and chemical properties including soil alkalinity, soil water holding capacity, and porosity, and alters the environment where microbial decomposition occurs (Francini et al 2014, Hiltbrunner et al 2012, Li et al 2011, Steffens et al 2008. This process, accompanied by the scorching effect, is described as a flux from aboveground C pools to the litter or soil C pools: where B is the biomass loss from trampling and scorching effect (g per day), n is the livestock density (livestock unit m −2 ) and p is the effect coefficient, equals to 0.008 (per livestock unit) (Vuichard et al 2007b).

Animal energy cycling
Three of the five models (SM, IOM, and the SAVANNA model) have explicit animal energy storage pools (IOM uses an animal fat reservoir instead of a weight pool), which describe the balance between the energy inputs and outputs, and the potential influence on grazing intake. The other two models (PASIM and DRM) specifically describe energy inputs and outputs but have no stated variables for the animal mass dynamics.
In the models, the animal energy dynamics are represented by a balance equation (table 1). The base expenditure or respiration and excretion are generally defined as the major energy outputs, both of which are represented as set percentages from the intake. The SAVANNA model includes additional outputs like milk production and travel that were not considered in this study.

Data used
The IGBP landcover map (Loveland et al 2000) was used to identify the grassland areas in TES. The 'gridded livestock of the world' (GLW Ver. 2.0, http://livestock.geo-wiki.org/) data from the FAO's Animal Production and Health Division were used to define the livestock distribution in the study region. The data were at a spatial resolution of 0.008333 • . It contains the livestock densities of sheep, goat and cattle in 2006. We converted the total livestock number into the sheep unit based on the FAO conversion rate (table S1 available at stacks.iop.org/ERL/13/075005/ mmedia).

Model simulations
In order to compare the grazing models under the same platform, all five grazing models were integrated into BEPS. We coded the grazing models following the equations from the model introduction papers or tech notes (table 1). Some essential simplifications were made for the SAVANNA model simulations as some of the designed processes, including animal travel, lactation and the sex-related factors, are difficult to be simulated at the regional scale due to data limitation. We kept the original parameterizations in each model but did essential modifications to the livestock-type related parameters such as the mature mass (AHFT) in IOM and the half maximum intake rate ( ) in IOM and the SAVANNA model. The parameters were modified to fit the sheep unit (table S2).
The age-related factor was considered in IOM, so we undertook four simulations with different age-based parameterizations (i.e. age = 0-1, 1-2, 2-3 and >3). For DRM, the original set of grazing efficiency (GE) is very high for pasture systems (No'am et al 1992), so we set three simulations with different GE values (i.e. high efficiency, GE = 670 m 2 per day per sheep unit, medium efficiency, GE = 450 m 2 per day per sheep unit, low efficiency, GE = 200 m 2 per day per sheep unit). For the SAVANNA model, the initial weight was not given from its tech-note or the other relevant publications, so we also set three levels of this parameter (i.e. high initial weight, W ini = 80 kg per sheep unit, medium initial weight, W ini= 60 kg, low initial weight, W ini = 40 kg) to perform the simulations. In total, 12 simulations were explicitly run for the BEPS model coupled with grazing models in this study.
In all simulations, the coupled models were run using the same forcing data and livestock distribution map in 2006. It is assumed that the weight updates were all from grazing activity throughout the year.
C graze in this study was calculated as the difference between C intake and C return from the excrete (C ex ):

Sensitivity analysis
A parameter sensitivity analysis was conducted following the method from Williams et al (2012): As a percentage of livestock density As a percentage of animal intake (Coughenour 1993) a Variables in animal biomass intake section: biomass above,ava : available aboveground biomass (g m −2 ); r intake : actual intake rate as a fraction of animal weight; W: current weight of the livestock (g per livestock unit); r d : intake rate when aboveground biomass is deficient; GE ∶ grazing efficiency (m 2 d −1 per livestock unit); I max : maximum livestock intake rate (g per livestock unit per day); q: intake parameter, K: LAI related intake parameter (m 2 /m 2 ); I max : maximum livestock intake rate (g per livestock unit per day), V consume : consumable aboveground biomass for the livestock (g m −2 ); : half maximum intake rate (g m −2 ); F j : the ratio of the fat depot (g), F max : the maximum of the fat storage (g), I: the actual daily intake (g); E D : daily energy expenditure (  where |S| is the response coefficient equivalent to relative change in C graze (i.e. Δ ) caused by the relative change in a parameter (i.e. Δ ). All the parameters in each model were tested in this analysis. For each test, we imposed a 2% increase in the parameter magnitude to calculate the |S| value. A specific parameter was considered to have high impact on model output if the |S| value exceeds 0.2.

The seasonal ensemble and trends
The monthly outputs of C graze generally agreed with the seasonal pattern of grazing activity in reality. All models predicted the high C graze during the growing season. As constrained by the forage availability, the seasonal trend of C graze is in accordance with the ANPP trend ( figure 1(b)). The grazing consumption is mainly concentrated from May to Sept. The mean and median values of monthly C graze increased with the forage availability from January to August and decreased from August-December. The mismatch of months with ANPP and C graze peaks reflects the continuous increase in forage demand (i.e. fast weight increase during the growth season) after the vigorous period of plant growth. In spite of the common predicted seasonal trends, C graze simulated by different models varied significantly. The monthly ranges and inter-quartile range (IQR, the difference between the 25th and 75th percentiles) of the distribution increased with the increase of forage availability during the growth season. As a result, the range and IQR of the predicted C graze in 2006 are from 0.05-16 gC m −2 yr. IOM specifically includes the age-dependent parameterization. With different sets of average livestock age, the predicted C graze values by IOM were from 3.5 gC m −2 yr (IOM_0_1) to 8.67 gC m −2 yr (IOM_3more), representing a range of more than three-fold.
3.2. Regional uncertainties in C graze and C graze as a percentage of ANPP and NPP The regional distribution of standard deviation (SD) and coefficient of variance (CV) of C graze are presented in figure 2. Spatial evaluation revealed that the absolute uncertainties in simulated C graze , indicated by SD, is higher in areas with higher livestock density such as the Russian grasslands, Southern Kazakh Steppe, Northern Mongolia, and Inner Mongolia, China. The CVs, indicating the relative uncertainties, are relatively low for many areas with high livestock density due to their high C graze , while some of the lightly grazed areas stand out, such as the grasslands in Kazakhstan and Mongolia. Only some of the grasslands in Inner Mongolia and the North Kazakh Steppe showed high values of both absolute and relative uncertainties.
C graze as a percentage of ANPP (C graze %ANPP) and NPP (C graze %NPP) showed similar spatial patterns (figure 3). Areas with high SD are located in the Russian Steppe, the South Kazakh Steppe, Inner Mongolia and some grasslands in Mongolia. Areas with low SD are mainly located in grasslands in Kazakhstan and the eastern part of Mongolia. Similar to the CV distribution pattern of C graze , high CVs in C graze % were found in the major grasslands of Kazakhstan and Mongolia. The range of CV is also similar to that of C graze .
As shown in figure 4, predicted C graze ranges from close to zero to more than 14 gC m −2 yr. The difference between the highest and the lowest values is more than 200-fold. The corresponding ratios of C graze %ANPP and C graze %NPP range from 0.08%-24.6% and 0.028%-11.2%, respectively. The IQR of predicted C graze ranges from 1.25 gC m −2 yr to 16.5 gC m −2 yr, which is much narrower compared to the condition with all model results included. The IQRs of C graze %ANPP and C graze %NPP range from 3.06%-25.5% and 1.12%-8.22%, respectively.

Sensitivity analysis
The sensitivity of simulated C graze to model parameters is presented in table 2. The models with explicit representation of the animal energy pool have more parameters than the ones that do not represent the pool. DRM uses a first order differential equation to represent the intake process and the weight update was not explicitly included, so it is not surprising that the C graze output is almost linearly correlated with GE (S = 0.98). SM uses a linear piecewise function to estimate C intake , while the sensitivity of r intake is partly dispersed to the parameters representing the animal weight updates. For the models using the M-M formulation, the half maximum intake rate ( in IOM and the SAVANNA model and equivalent to K q in PASIM) is the common sensitive parameter. The S values for are similar in IOM and the SAVANNA model, while it is much higher for PASIM as it uses an exponential function. I max is a sensitive parameter in both PASIM and the SAVANNA model. It is probably sensitive in IOM as determined by AHFT, which is a sensitive parameter affecting C graze output. The C graze output of the SAVANNA model and IOM are also very sensitive to the weight-related parameters.  , C x : potential intake rate of livestock (kg d −1 per livestock unit); b PASIM, I max : maximum livestock intake rate (g per livestock unit per day); q: intake parameter, K: LAI related intake parameter (m 2 m −2 ); c SM, W ini : initial weight at the beginning of the year (g per livestock unit), r intake : actual intake rate as a fraction of animal weight, W max : maximum weight over a year (g per livestock unit), R: daily respiration as a fraction of animal weight, F lb : digestibility of live biomass as a fraction of live biomass intake, F db : digestibility of dead biomass as a fraction of dead biomass intake. r d : intake rate when aboveground biomass is deficient. d IOM, AHFT: mature mass of livestock(g), m: metabolic coefficient (MJ kg fat −1 ), : half maximum intake rate (g m −2 ), R me : coefficient of metabolizable energy of grass (MJ kg −1 ), M j : ratio of body mass in an age class from mature mass, F j : the ratio of the fat depot (g), F ,max : ratio of maximum fat storage in an age class from mature body mass; F db : digestibility of dead biomass as a fraction of dead biomass intake. e SAVANNA model: I max : maximum livestock intake rate (g per livestock unit per day); : half maximum intake rate (g m −2 ), W ini : initial weight at the beginning of the year (g per livestock unit), W max : maximum weight over a year (g per livestock unit), W min : minimum weight over a year (g per livestock unit), R me : coefficient of metabolizable energy of grass (MJ kg −1 ), R: daily respiration rate.

Discussions
C graze , representing C removal by livestock, is the fundamental indicator of how grazing activity influences vegetation growth, production, and C sequestration. Although it has been widely studied in this region, the results have not been previously synthesized and compared. As presented in this study, great uncertainties exist in modeling C graze from the pasture systems. Different models output diverse results of C graze in both range and IQR analysis ( figure 3). The high estimations of C graze in this study are comparable to the results based on the grazing gap approach, while the other simulations produce much lower values (table S3). The uncertainty of process-based terrestrial modeling could be attributed to three sources: the data inputs, the algorithms to describe the processes, and the corresponding parameterization.
As the same forcing data were used in this study, the potential uncertainty from data input in the model inter-comparison has been excluded in this study. However, the choice of forcing dataset to estimate C graze would cause biases to the outputs from the real values. The livestock distribution map and the datasets used in the BEPS are the two major inputs to predict C graze in this study. Derived from the official statistics, the GLW dataset could well constrain the total amount of livestock in this region. So the main bias from this dataset would be the detail algorithms to predict the livestock distribution in specific areas (Robinson et al 2014). On the other hand, the BEPS predicts lower NPP in TES comparing to the other models (table 3). One important source for the difference should be the model input of LAI. Unstable surface reflectance, LAI saturation and soil effects in sparse vegetation region are the major sources of deviation in LAI datasets. The latter two factors were reported to cause LAI underestimation (Zhao et al 2012). Liu et al (2012) also reported the lower estimations of grassland LAI in the globalLAI dataset (i.e. the LAI product used in this study) due to those factors. Other factors including the choice of meteorological data and different spatial scales used to perform studies would also lead to the different results in this comparison. As summarized in table 1, the major grazing-related processes have been included and simulated in the current models. However, some processes like animal respiration, trampling and excretion and the related CH 4 and N 2 O emissions are only defined in a conceptual style. Neither of these processes have been fully investigated, nor well validated with enough observations. The simulation of these processes change the energy demands and so affect the forage intake estimation. Regarding to the intake process, forage availability is considered as the major factor and energy demand is used as a constrain factor. The other potential factors such as forage quality and environmental constraints on the intake rate are designed in some models such as IOM and the SAVANNA model (Coughenour 1993, Illius andO'Connor 2000). However, they are difficult to be used at large scales due to the lack of enough observations and useable datasets (Galyean 2016, Senft et al 1985. Another process-based modeling gap is how to consider the feedback of vegetation to grazing activity. Current models regard grazing activity as a direct C removal process, but its final effect on pasture ecosystems should largely depend on the relative growth rate of grasslands (Hilbert et al 1981, Oesterheld andMcNaughton 1991a). During a growing season with suitable environmental condition, grazing help clean the older tissue and detritus (Matches 1992, Milchunas and Lauenroth 1993, Oesterheld and McNaughton 1991b, stimulate the reallocation of carbohydrates and nutrient (van Staalduinen and Anten 2005, Volenec et al 1996, Zhao et al 2008, and so facilitate the regrowth of grasslands. In contrast, grazing accelerates the grass senescence during the drought periods, in addition to the C directly removed (Schönbach et al 2011, Wan andSosebee 2002). This neglect of the seasonal variations of compensatory regrowth should cause further deviations in modeling C graze .
As indicated by the parameter sensitivity analysis, the C graze outputs of different models are highly sensitive to the set of intake rate or maximum intake rate. In DRM and the SAVANNA model, the |S| value of the corresponding parameter is close to 1, indicating its primary importance in predicting C graze under the current model framework. Without explicitly considering the animal energy factor, C graze prediction in DRM and PASIM mainly depends on the two or three parameters that are included. Thus, the outputs of the two models largely rely on the observability and predictability of these parameters. In the models based on the M-M formulation, the half maximum intake rate ( ) is another important source of uncertainty. The contrasting set of in PASIM to IOM and the SAVANNA model could be an important reason of its lower C graze prediction among all the simulations (table S2 and figure 1). Moreover, the outputs of IOM reflect the potential importance of the age-related parameterizations. From our results, the range of the regional C graze could be more than three-fold under different average livestock age-levels (figure 1). Livestock age determines the animal weight, energy expenditure, and energy storage ability (Ådnøy et al 2005). In a highly managed system, it also closely relates to the local slaughtering strategy, which connects to the social and economic aspects of the administrate divisions in the region (Murphy et al 2017, Pacheco et al 2015. For example, the typical slaughter age of Kazakh sheep for meat is from 12-60 months, while the corresponding slaughter age in Inner Mongolia is much shorter due to increasing meat demands (Kosilov et al 2015, Qi et al 2017, Shi and Cui 2015, Torgerson et al 2009, Xiaobing et al 2015. However, this factor is not included or simplified by using a maximum weight in most grazing models. It would also lead to a potential issue in multi-year continuous simulations in the models with an explicit animal weight pool. That is, how to connect the weight decrease from the relatively high weight per unit animal at the end of growing season to the relatively low weight per unit animal at the non-growing season due to slaughtering.
Furthermore, proper parameterization would determine the correct prediction over the current model framework. Based on our synthesis, many of the sensitive parameters are pre-set from a priori knowledge or existing regional or sub-regional and inventory or publication-derived data. Unfortunately, this parameterization strategy is likely the best approach in the current grazing modeling as limited observations have been implemented. Some of the sensitive parameters have not been well noticed in previous studies, such as the half maximum intake rate ( ) in the models using the Michaelis-Menton (M-M) formulation, animal metabolic coefficient (m) in IOM and q and K in PASIM and so on. While for those widely observed parameters such as intake rate and maximum intake rate, the outputs have not been well synthesized for large scale studies (Glindemann 2007, Liu et al 2015a, Müller et al 2014, Maestre et al 2012. These data gaps hamper the rigid data-model inter-comparison over the region. In addition, it is also difficult to perform the data-model fusion as the parameters in vegetation or soil processes due to the lack of continuous and standardized observations. Therefore, our study calls for much closer cooperation on these potentially important spots among the animal science, animal husbandry, and ecosystem modeling communities. The current quantitative knowledge of grazing GHG emission is mainly from inventory-based assessments, a top-down grazing gap approach (Herrero et al 2013, Krausmann et al 2008. The estimation of grazing C emission is 1.32 Gt CO 2 -eq per year, which is even higher than the total C sequestration from pasture systems globally (Henderson et al 2015, Lal 2004). If the current estimation of grazing emission is true, then it should be necessary to consider the explicit representation of grazing activity in ESMs. However, our current analysis indicated the large discrepancies in modeling grazing impact on terrestrial C cycling using process-based models. Both the key process descriptions and parameterizations need further tests and calibrations. To facilitate model development, we propose incorporating: (1) a comprehensive synthesis of usable dataset to calibrate the sensitive parameters; and (2) a standard observation system to evaluate and benchmark models.

Conclusions
In this study, we performed a multi-model intercomparison analysis of process-based grazing models. The results suggest that large uncertainty exists in modeling grazing impact to C sequestration in TES. Our study identifies key issues in descriptions of important processes such as animal respiration, trampling and excretion and after-graze response, and parameterizations such as intake rates and age-related parameters. The modeling gaps dominate the current uncertainty in estimating C graze . The findings in this study call for more inter-community cooperation on parameter calibration and model evaluation, including the comprehensive synthesis of usable dataset, the foundation of a standard observation system and the observe-based inter-comparison to evaluate models, which would facilitate more accurate assessment of C sequestration of pasture ecosystems and better representation of C sequestration in earth system models.

Acknowledgments
This work is supported by the National Key R&D Program of China (2016YFA0600202), National Youth Science Fund of China (41701227) and the priority academic program development of Jiangsu higher education institutions (PAPD).
We acknowledge all producers and providers of datasets we used in this study. We also greatly acknowledge Dr Adrian Pachzelt from Goethe-University Frankfurt, Germany for helping prove the model codes and give useful suggestions, Dr Pavel Groisman from NOAA and Dr Pavel Propasin from Georg-August University Göttingen, Göttingen, Germany for giving suggestions to the prototype manuscript and Dr Yiqi Luo from North Arizona University for giving revision suggestions. This study is largely inspired by the topics and discussions during the recent NEESPI/NEFI events. Without this kind help and encouragement, this work could not be finished.