A Satellite-Based Model for Simulating Ecosystem Respiration in the Tibetan and Inner Mongolian Grasslands

It is important to accurately evaluate ecosystem respiration (RE) in the alpine grasslands of the Tibetan Plateau and the temperate grasslands of the Inner Mongolian Plateau, as it serves as a sensitivity indicator of regional and global carbon cycles. Here, we combined flux measurements taken between 2003 and 2013 from 16 grassland sites across northern China and the corresponding MODIS land surface temperature (LST), enhanced vegetation index (EVI), and land surface water index (LSWI) to build a satellite-based model to estimate RE at a regional scale. First, the dependencies of both spatial and temporal variations of RE on these biotic and climatic factors were examined explicitly. We found that plant productivity and moisture, but not temperature, can best explain the spatial pattern of RE in northern China’s grasslands; while temperature plays a major role in regulating the temporal variability of RE in the alpine grasslands, and moisture is equally as Remote Sens. 2018, 10, 149; doi:10.3390/rs10010149 www.mdpi.com/journal/remotesensing Remote Sens. 2018, 10, 149 2 of 20 important as temperature in the temperate grasslands. However, the moisture effect on RE and the explicit representation of spatial variation process are often lacking in most of the existing satellite-based RE models. On this basis, we developed a model by comprehensively considering moisture, temperature, and productivity effects on both temporal and spatial processes of RE, and then, we evaluated the model performance. Our results showed that the model well explained the observed RE in both the alpine (R2 = 0.79, RMSE = 0.77 g C m−2 day−1) and temperate grasslands (R2 = 0.75, RMSE = 0.60 g C m−2 day−1). The inclusion of the LSWI as the water-limiting factor substantially improved the model performance in arid and semi-arid ecosystems, and the spatialized basal respiration rate as an indicator for spatial variation largely determined the regional pattern of RE. Finally, the model accurately reproduced the seasonal and inter-annual variations and spatial variability of RE, and it avoided overestimating RE in water-limited regions compared to the popular process-based model. These findings provide a better understanding of the biotic and climatic controls over spatiotemporal patterns of RE for two typical grasslands and a new alternative up-scaling method for large-scale RE evaluation in grassland ecosystems.


Introduction
Grasslands are one of the most widespread vegetation types and store one-fifth of the total global carbon in its vegetation and soil [1,2]; thus, they play an important role in the global carbon cycle.The extensive grasslands in China make up ~10% of the total area of grasslands in the world and have been estimated to store 9-15% of the total carbon held in the grasslands worldwide [3].Northern China's grasslands located in the Tibetan Plateau and Inner Mongolian Plateau, in particular, constitute the majority (more than 70%) of the grasslands in China and represent two significant grassland types worldwide (i.e., alpine and temperate grasslands) [4].Moreover, they are sensitive to climate changes due to their unique plateau topography, the extreme cold, arid and semi-arid ecological environment and the high soil carbon density [5][6][7].Therefore, a better understanding of the carbon cycle in northern China's grasslands is necessary for it to serve as a sensitivity indicator of regional and global carbon cycles.
Ecosystem respiration (RE) is an important outflow component of ecosystem carbon cycle and represents the second largest carbon exchanges between terrestrial ecosystems and the atmosphere [8].However, as the result of the complex interactions among physical, chemical, and biological processes, RE varies greatly at different temporal and spatial scales [9].Extensive studies at regional or global scales suggest that RE estimates still remains a large uncertainty [10].Therefore, it is particularly essential but challenging to accurately evaluate the regional RE in northern China's grasslands for a quantitative assessment of the terrestrial ecosystem carbon budget and its response to future climate changes.
Generally, modeling is the commonly adopted method for inferring RE at a large scale.Remote-sensing models, with relatively simple model structures and high spatiotemporal resolutions, provide promising tools for monitoring regional RE.Currently, with the development of fine-scale remote-sensing technology and long-term regional and global FLUXNET observations, extensive studies have revealed the robust statistical and physiologically meaningful relationships between RE and satellite products, e.g., the land surface temperature (LST) [11][12][13] and vegetation index (VI) [14][15][16].Based on this information, many semi-empirical satellite-based models were developed and well validated at the plot scale [13,[16][17][18][19].However, RE at large scale involves complex spatial and temporal variation processes [20].Few satellite-based models have an explicit representation of the spatial process of RE.Specifically, the basal respiration rate, an important parameter in RE models to describe the location-specific respiration rate, is usually set as a constant parameter, which introduces considerable error in the simulated spatial variability of RE [8,20,21].Therefore, it is critical to identify the controlling mechanism involved in the spatial variations of RE with more comprehensive data, not just in the temporal variations, thus serving as a valuable basis to better understand patterns of RE at regional scale and to scale-up from specific sites to vegetation biomes.In addition, both climatic and biotic controls over RE should be considered comprehensively in RE simulation [22,23].In satellite based models, RE has often been linked to temperature or plant productively solely.Some studies advanced RE estimates by accounting for the combined two effects of temperature and plant productivity [13,15].The moisture limitation on RE is often lacking in most of the current satellite-based respiration models, leading to a potential underestimation of drought effects on respiration (e.g., [19,24]).However, grasslands are one of the most sensitive ecosystems to the alteration of water regimes [25].Water stress was found to have a strong influence on northern China's grasslands, especially in southwestern Inner Mongolia and northwestern Tibet [26].Recently, some satellite-based models gradually began incorporating the moisture effect on RE, but it still depends which satellite product and form are selected to represent the moisture factor and the response process of RE to moisture.For instance, a linear equation of an indirect water indicator (diurnal LST difference, LST diff ) was utilized in the RECO model during the dry season for water-limited biomes, which were mostly concentrated in the Mediterranean climate region without a temperature limitation [27]; thus, its utility in northern China's grasslands still needs to be validated.Root-zone soil moisture derived from the Advanced Microwave Scanning Radiometer-Eos was adopted in a soil respiration (Rs) model in the Midwest USA, and it yielded a lower explanation capacity for seasonal variation of Rs than the surface soil moisture under extreme moisture regimes [28].Currently, increasing attention is being given to the moisture effect on RE [29,30], as most climate change scenario analyses have predicted that the frequency and severity of drought will increase in most biomes across the world [31].It is practical to develop a satellite-based RE model that integrates temperature, plant productivity, and moisture effects together in both temporal and spatial processes, to yield reliable simulations of regional RE for northern China's grasslands.
This study was based on 57 site-year flux-tower data and remote-sensing data from northern China's grasslands.The major objectives of this study were threefold: (i) to quantify the relationships of climatic and biotic factors and the spatial and temporal variability of RE in the Tibetan alpine grasslands and the Inner Mongolian temperate grasslands; (ii) to develop a high-resolution satellite-based model for simulating RE at regional scale on the basis of comprehensive consideration of vegetation productivity, temperature, and moisture; and (iii) to assess the model performance in terms of the temporal dynamics as well as the spatial pattern simulation of RE in northern China's grasslands.

Study Area
This study was conducted in the alpine grasslands on the Tibetan Plateau and the temperate grasslands on the Inner Mongolian Plateau.The Tibetan region is in the alpine climate zone, where temperature appears as the controlling factor for the grassland ecosystems.With an average elevation of ca.4000 m, the mean annual air temperature on the plateau ranges between −5.75 and 2.57 • C, and the mean annual precipitation varies from 200 to 600 mm.The alpine meadow covers more than half of the plateau, representing not only a typical ecosystem in central Asia alpine environment but also a unique ecosystem within the alpine regions of the world [32].Located in the arid and semi-arid zone, the grassland ecosystems in the Inner Mongolian region are mainly limited by moisture.With an average elevation ca.1000 m, the mean annual temperature in this area varies between 3 and 6 • C, and the mean annual precipitation ranges from 200 to 350 mm.The temperate grasslands are a typical vegetation type under the temperate continental climate and they represent an important component of the Eurasian grasslands [33].According to the Atlas of China's Grassland Resources

Flux-Tower Data
A total of 57 site-year data from 16 grassland sites during 2003-2013 were retrieved from the ChinaFLUX and the Coordinated Observations and Integrated Research over Arid and Semi-Arid China (COIAS) databases, including 43 site-years from 11 sites on the Tibetan Plateau and 14 site-years from 5 sites on the Inner Mongolia Plateau (Table 1).The observed eddy-covariance data were processed through a three-dimensional coordinate rotation, WPL (Webb-Pearman-Leuning) correction, and invalid data exclusion [34].Subsequently, the missing nighttime RE and daytime RE data were calculated using the Lloyd-Taylor equation based on net ecosystem exchange (NEE) observations during the nighttime [35].The entire procedure was completed using the ChinaFLUX CO2 data processing system [36].The half-hourly flux data was summed to obtain daily values, and the site-years with more than 30% of the daily RE missing were eliminated.To comply with the temporal scale of the Moderate Resolution Imaging Spectroradiometer (MODIS) 8-day composite imagery, the processed daily RE were averaged within the same periods.

Flux-Tower Data
A total of 57 site-year data from 16 grassland sites during 2003-2013 were retrieved from the ChinaFLUX and the Coordinated Observations and Integrated Research over Arid and Semi-Arid China (COIAS) databases, including 43 site-years from 11 sites on the Tibetan Plateau and 14 site-years from 5 sites on the Inner Mongolia Plateau (Table 1).The observed eddy-covariance data were processed through a three-dimensional coordinate rotation, WPL (Webb-Pearman-Leuning) correction, and invalid data exclusion [34].Subsequently, the missing nighttime RE and daytime RE data were calculated using the Lloyd-Taylor equation based on net ecosystem exchange (NEE) observations during the nighttime [35].The entire procedure was completed using the ChinaFLUX CO 2 data processing system [36].The half-hourly flux data was summed to obtain daily values, and the site-years with more than 30% of the daily RE missing were eliminated.To comply with the temporal scale of the Moderate Resolution Imaging Spectroradiometer (MODIS) 8-day composite imagery, the processed daily RE were averaged within the same periods.

Remote-Sensing Data
At the site scale, the MODIS products of the LST, enhanced vegetation index (EVI), and land surface water index (LSWI) at each site-year (2003-2013) were downloaded from the University of Oklahoma Data Center (http://www.eomf.ou.edu/visualization/manual/), based on the pixels where the flux towers located, to match up with the flux data for model parametrization.
At the spatial scale, the 8-day (best observation in 8 days) MODIS land surface reflectance datasets (MOD09A1, Level 3, Collection 5) during 2001-2010 were downloaded from NASA Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/) to generate the spatial EVI and LSWI data using Equations ( 1) and (2), as model drivers in regional simulation.
The MOD11A2 and MYD11A2 products (Level 3, Collection 5), which are 8-day average values of cloud free observations, were also downloaded from the LP DAAC data pool to obtain the spatial daytime LST and nighttime LST data, respectively.The spatial resolution of EVI and LSWI is 500 m, and that of LST is 1 km.For the regional application in China's grasslands, the extracted EVI and LSWI were resampled to a 1-km resolution.The temporal resolution of these products is 8-day, which ensures a better data quality as selected on the basis of the absence of clouds or aerosol, etc., and is confirmed by extensive RE modeling research to be able to robustly capture the seasonal variation of RE [24,27,28].Additionally, to further reduce the effect of cloud and to capture the seasonality of EVI, the original time series were smoothed using the double logistic curve fit in the TIMESAT software (Lund University, Lund, Sweden) [38].

Model
RE at the regional scale is a coupling result of complex spatial and temporal processes.The process of spatial variability refers to the basic trend of RE changing with the climate and substrate at a specific location, which is especially significant in regional simulations [39]; the process of temporal variability refers to the cyclical changes of RE induced by the seasonal dynamics of environmental factors and plant growth.A multiplication factor is usually used to express the interaction of the above two processes in an RE simulation [20][21][22]27].Based on the remote-sensing model RECO [27] and the traditional climate-driven model T&P [20,22], we put forward a modified model of RE for grasslands by taking the effects of temperature, vegetation productivity, and moisture on both the spatial and temporal variability of RE into account and making full use of the high-resolution satellite data as the model drivers.To account for spatial and temporal variability separately, we partition RE into a site-specific reference respiration rate (Reref ) and the remaining seasonal variation (f (T,P,W)), which is detrended for site characteristics (Equation ( 3)).The structure of the model is shown below (Figure 2).
Remote Sens. 2018, 10, 149 6 of 20 seasonality of EVI, the original time series were smoothed using the double logistic curve fit in the TIMESAT software (Lund University, Lund, Sweden) [38].

Model
RE at the regional scale is a coupling result of complex spatial and temporal processes.The process of spatial variability refers to the basic trend of RE changing with the climate and substrate at a specific location, which is especially significant in regional simulations [39]; the process of temporal variability refers to the cyclical changes of RE induced by the seasonal dynamics of environmental factors and plant growth.A multiplication factor is usually used to express the interaction of the above two processes in an RE simulation [20][21][22]27].Based on the remote-sensing model RECO [27] and the traditional climate-driven model T&P [20,22], we put forward a modified model of RE for grasslands by taking the effects of temperature, vegetation productivity, and moisture on both the spatial and temporal variability of RE into account and making full use of the high-resolution satellite data as the model drivers.To account for spatial and temporal variability separately, we partition RE into a site-specific reference respiration rate (Reref) and the remaining seasonal variation (f(T,P,W)), which is detrended for site characteristics (Equation ( 3)).The structure of the model is shown below (Figure 2).

Representation of Spatial Variability of RE
Reref is the site-specific respiration rate at the reference temperature (Tref).It is used to describe the difference in basal respiration between sites and years, thereby reflecting the spatial variability of RE [8].Numerous studies have indicated that a spatially heterogeneous Reref results in a better estimation of RE at different sites [20][21][22]27].The long-term spatial variation of RE or Reref is found to be mainly driven by the mean annual temperature [40][41][42], mean annual plant productivity [43,44], and mean annual moisture [45][46][47].Thus, the functional description of Reref in this study is developed based on its dependencies on the temperature, productivity, and moisture factor (Figure 3).Since Reref becomes less suitable to capture site characteristics at higher Tref values [27], Tref was set as the mean springtime (DOY 96-144) temperature (i.e., 3 °C for the Tibetan alpine grasslands, and 7 °C for the Inner Mongolian temperate grasslands), and the Reref for each site-year was calculated as the mean respiration rate at the temperature approximate to Tref during the spring time.Finally, Reref is estimated in our model as follows:

Representation of Spatial Variability of RE
Reref is the site-specific respiration rate at the reference temperature (Tref ).It is used to describe the difference in basal respiration between sites and years, thereby reflecting the spatial variability of RE [8].Numerous studies have indicated that a spatially heterogeneous Reref results in a better estimation of RE at different sites [20][21][22]27].The long-term spatial variation of RE or Reref is found to be mainly driven by the mean annual temperature [40][41][42], mean annual plant productivity [43,44], and mean annual moisture [45][46][47].Thus, the functional description of Reref in this study is developed based on its dependencies on the temperature, productivity, and moisture factor (Figure 3).Since Reref becomes less suitable to capture site characteristics at higher Tref values [27], Tref was set as the mean springtime (DOY 96-144) temperature (i.e., 3 • C for the Tibetan alpine grasslands, and 7 • C for the Inner Mongolian temperate grasslands), and the Reref for each site-year was calculated as the mean respiration rate at the temperature approximate to Tref during the spring time.Finally, Reref is estimated in our model as follows: where EVI mean is the mean annual springtime EVI, LST mean is the mean annual daytime LST (LST d ), LSWI mean is the mean annual growing season LSWI (DOY 136-272), and p 1 to p 4 are parameters to be inferred herein.Figure 3 displays the Reref dependencies on EVI mean , LST mean and LSWI mean .

Quantitative Relationships between RE and Biotic and Climatic Factors
The observed RE values were partitioned into Reref and f(T,P,W) to account for the spatial and temporal variability of RE separately.The quantitative relationships between these two components and the climatic and biotic factors provided a clear physiological basis for the development of the model.

Quantitative Relationships between Reref and Biotic and Climatic Factors
Multiple regression analyses between Reref and the remote-sensing data show that Reref has a significant dependence on EVImean, LSTmean and LSWImean, which can be approximated linearly (Figure 3).The variation of Reref in the Tibetan alpine grasslands is primarily attributed to the change in vegetation productivity (EVImean, R 2 = 0.612) and to the water regime (LSWImean, R 2 = 0.334).Temperature plays a negligible role in the Reref (LSTmean, R 2 = 0.128).Similarly, productivity (R 2 = 0.576) and moisture (R 2 = 0.480) tend to be more important than temperature (R 2 = 0.461) for Reref in the temperate grasslands of Inner Mongolia.However, the temperature limitation among the temperate grassland sites is characterized by a negative slope.

Representation of Temporal Variability of RE
The seasonal dynamics of RE are reflected by f (T, P, W), which is mainly controlled by the temporal variability of temperature [48,49], vegetation productivity [43,50], and water [20][21][22]40].Especially in the arid and semi-arid ecosystems, water is the main limiting factor of the seasonality of RE [51,52].The values of f (T,P,W) are obtained by calculating the ratio of the observed RE from the flux tower and the site-specific Reref partitioned from the observed RE. Figure 4 displays the f (T, P, W) dependencies on the temporal variation of EVI, LST and LSWI.The form of f (LST, EV I) scalar was based on the LST-driven Lloyd-Taylor equation [48] and EVI-driven linear equation, while the form of fw scalar was derived from a Michaelis-Menten equation firstly proposed by Raich et al. [53].During the growing season, we calculate the f (T, P, W) as follows: where LST n is nighttime LST, EVI std is the standardized EVI (8-day EVI divided by EVI mean ), p 5 to p 7 are the parameters to be sought, E 0 is the activation energy, T 0 is the minimum temperature for respiration, which is set to 227.13 k (−46.02• C) as in the original Lloyd-Taylor model, k is the half-saturation constant of the hyperbolic relationship, and 0.5 is an empirically given value to maintain the numerator/denominator of fw > 0.
Remote Sens. 2018, 10, 149 9 of 20 3.1.2.Quantitative Relationships between f(T,P,W) and Biotic and Climatic Factors The regression analyses between f(T,P,W) and the remote-sensing indexes (Figure 4) indicate that the temporal variability of RE has a strong exponential dependency on temperature, which could be described using the LSTn-driven Lloyd-Taylor equation.RE still has a significant linear response to the seasonal dynamics of EVI, which is decided by the behind mechanism that growth respiration supplies energy at the cost of consuming some photosynthate, causing the RE to be approximate to a constant fraction of plant productivity [57,58].In terms of the response to temporal variation of LSWI, RE is almost linearly related to LSWI firstly (when LSWI < 0.1) and then it reaches saturation during the final stage in the Inner Mongolian region with a relatively low soil water content, whereas in the Tibetan region with a relatively high soil water content, the response of RE to LSWI follows a gradual saturation.The responses of RE to moisture under the two different water regimes both agree well with the typical Michaelis-Menten curve.The temporal variation of RE in the alpine grasslands of the Tibetan Plateau is mainly correlated to temperature (LSTn, R 2 = 0.752), followed by plant productivity (EVIstd, R 2 = 0.617), while water (LSWI, R 2 = 0.451) plays a minor role in regulating the temporal dynamics of RE in this temperature-limited region (Figure 4).In contrast, water exerts a strong influence (LSWI, R 2 = 0.656) on the temporal variability of RE in the temperate grasslands of Inner Mongolia, which is comparable with the determinant factor (LSTn, R 2 = 0.687), while plant productivity exerts a relatively weak impact (EVIstd, R 2 = 0.485).

Model Parameterization and Validation
Based on the quantitative dependencies of spatiotemporal variability in RE to these biotic and environmental factors, we developed a semi-empirical RE model described in 2.3.The long-term flux tower-observed RE was used for parameter estimation and validation.Detailed parameter values for Reref and f(T,P,W) in the two grasslands are presented in Table 2.These parameters are mostly empirical coefficients varying with the vegetation types and structures, site conditions, and histories.One of the parameters for Reref, for example, p3, has the opposite sign for the temperate grasslands (negative) relative to the alpine grasslands (positive).This phenomenon implies a negative correlation between the spatial variation of RE and mean annual temperature in Inner During the non-growing season, temperature, but not water, becomes the main limiting factor for RE in northern China's grasslands [24].In addition, the surface covered by ice and snow in the non-growing season has a strong absorption of the shortwave infrared band that causes LSWI to increase sharply; thus, it cannot reflect the actual water content in the winter [54].Therefore, the fw scalar is only considered for the growing season (Equation ( 5)).During the non-growing season, we calculate the temporal variability of RE as follows: The different periods of MODIS indices (i.e., the mean springtime of EVI, the mean annual of LST d in Reref, and LST n in f (T, P, W)) used here have already been selected by Jägermeyr et al. [27] as the best index for fitting the relationship with RE.LSWI, selected by us as the model driver, is a more direct proxy for moisture than LST diff , which was used in RECO but had no significant relation with the temporal dynamics of RE in our study (R 2 < 0.2, p > 0.05); in addition, land surface moisture was found to have greater influence on RE than deep soil water [20,28], which is especially instrumental in northern China's grasslands with the large number of roots distributed in the soil surface [55].Therefore, all these drivers appear to be chosen empirically but are based on the physiological dependencies of RE, thereby ensuring reliable estimates from our model.

Model Parameterization and Validation
The two components, i.e.Reref for each site and year and the remaining seasonal variation f (T,P,W) partitioned from the overserved RE, were used to optimize parameters for Equations ( 4) and ( 5) separately to better constrain each component of the model.The nonlinear least-squares curve fitting method was utilized to estimate the model parameters.In addition to joint-sites estimation for the two grasslands in two regions, the leave-one-out cross validation method was employed to test the validity of the joint estimation of parameters (i.e., one subtype was excluded at a time, and the data for all other subtypes were used to estimate the parameters for each grassland type).The coefficient of determination (R 2 ) and root mean square error (RMSE) were chosen to assess the model performance.Regional-scale application is based on PFT-specific parameters retrieved from the entire site-data of alpine grasslands and temperate grasslands to assure most robust parameters.
To further verify the model improvement by incorporating the moisture effect on both spatial and temporal variability, the LSWI mean in Reref and fw in f (T,P,W) were removed, respectively, to test the accuracy of the model.In addition, we also compare some other semi-empirical or process-based models to assess the ability of our model.Specifically, for the RE estimation at a site scale, we chose the two original models that our model was based on, namely T&P [22] and RECO [27]; for the RE estimation at a regional scale, an independent regional RE product simulated by a mainstream dynamic global vegetation model CLM4CN, derived from Piao et al. [56], was applied.

Quantitative Relationships between RE and Biotic and Climatic Factors
The observed RE values were partitioned into Reref and f (T,P,W) to account for the spatial and temporal variability of RE separately.The quantitative relationships between these two components and the climatic and biotic factors provided a clear physiological basis for the development of the model.

Quantitative Relationships between Reref and Biotic and Climatic Factors
Multiple regression analyses between Reref and the remote-sensing data show that Reref has a significant dependence on EVI mean , LST mean and LSWI mean , which can be approximated linearly (Figure 3).The variation of Reref in the Tibetan alpine grasslands is primarily attributed to the change in vegetation productivity (EVI mean , R 2 = 0.612) and to the water regime (LSWI mean , R 2 = 0.334).Temperature plays a negligible role in the Reref (LST mean , R 2 = 0.128).Similarly, productivity (R 2 = 0.576) and moisture (R 2 = 0.480) tend to be more important than temperature (R 2 = 0.461) for Reref in the temperate grasslands of Inner Mongolia.However, the temperature limitation among the temperate grassland sites is characterized by a negative slope.

Quantitative Relationships between f (T, P, W) and Biotic and Climatic Factors
The regression analyses between f (T, P, W) and the remote-sensing indexes (Figure 4) indicate that the temporal variability of RE has a strong exponential dependency on temperature, which could be described using the LST n -driven Lloyd-Taylor equation.RE still has a significant linear response to the seasonal dynamics of EVI, which is decided by the behind mechanism that growth respiration supplies energy at the cost of consuming some photosynthate, causing the RE to be approximate to a constant fraction of plant productivity [57,58].In terms of the response to temporal variation of LSWI, RE is almost linearly related to LSWI firstly (when LSWI < 0.1) and then it reaches saturation during the final stage in the Inner Mongolian region with a relatively low soil water content, whereas in the Tibetan region with a relatively high soil water content, the response of RE to LSWI follows a gradual saturation.The responses of RE to moisture under the two different water regimes both agree well with the typical Michaelis-Menten curve.
The temporal variation of RE in the alpine grasslands of the Tibetan Plateau is mainly correlated to temperature (LST n , R 2 = 0.752), followed by plant productivity (EVI std , R 2 = 0.617), while water (LSWI, R 2 = 0.451) plays a minor role in regulating the temporal dynamics of RE in this temperature-limited region (Figure 4).In contrast, water exerts a strong influence (LSWI, R 2 = 0.656) on the temporal variability of RE in the temperate grasslands of Inner Mongolia, which is comparable with the determinant factor (LST n , R 2 = 0.687), while plant productivity exerts a relatively weak impact (EVI std , R 2 = 0.485).

Model Parameterization and Validation
Based on the quantitative dependencies of spatiotemporal variability in RE to these biotic and environmental factors, we developed a semi-empirical RE model described in 2.3.The long-term flux tower-observed RE was used for parameter estimation and validation.Detailed parameter values for Reref and f (T, P, W) in the two grasslands are presented in Table 2.These parameters are mostly empirical coefficients varying with the vegetation types and structures, site conditions, and histories.One of the parameters for Reref, for example, p 3 , has the opposite sign for the temperate grasslands (negative) relative to the alpine grasslands (positive).This phenomenon implies a negative correlation between the spatial variation of RE and mean annual temperature in Inner Mongolia.Among the parameters for f (T,P,W), the activation energy E 0 and the half saturation constant k are two key parameters that have a special physiological implication.In the alpine grasslands, E 0 has a higher value than the temperate grasslands.This result accurately reflects the increasing temperature sensitivity with the decreasing temperature.Regarding the half saturation constant, k, the value is higher (0.207) in the temperate grasslands than the alpine grasslands (0.021).The drier ecosystem has a higher half saturation point, accurately characterizing the different responses of RE to moisture among the alpine and the arid and semi-arid regions.The model parameter estimation reveals that approximately 77.8% and 74.9% of the variation in the observed RE across the alpine and temperate grasslands, respectively, could be explained by this model, albeit with a small bias (Figure 5; Table 3, joint-sites).Specifically, the model could explain more than 55% of the spatial variability of RE (Reref ) in the two regions.The simulation error in the Tibetan alpine grasslands (RMSE = 0.62 g C m −2 day −1 ) is higher than that in the Inner Mongolian temperate grasslands (RMSE = 0.11 g C m −2 day −1 ).The temporal variability (f (T,P,W)) is generally explained to a greater degree (more than 65%) compared to the spatial variability, while the simulation error in the temperate grasslands (0.91 g C m −2 day −1 ) is higher than that in the alpine grasslands (0.35 g C m −2 day −1 ).The cross-validation among subtypes in the alpine and temperate grasslands shows that approximately 75.6-91.7%and 57.5-80.3%,respectively, of the variation in the observed RE could be predicted by the model, with the values of RMSE ranging from 0.841 to 1.544 g C m −2 day −1 and 0.405-1.021g C m −2 day −1 for the alpine and temperate grasslands, respectively (Table 3).Except for the typical steppe, where the parameters derived from the limited sites in the meadow steppe and desert steppe do not represent for the typical steppe well, the other vegetation types differ slightly in performances from the cross-validation and parameter estimation, which verifies the robustness of the parameterization results.

Model Simulation for Seasonal and Inter-Annual Dynamics of RE at the Site Scale
Using the model parameters in Table 2, we simulate time-series RE of all 57 site-years and then further compared them to the flux tower-observed RE for each vegetation subtype in the alpine and temperate grasslands.In most cases, the temporal variation of the RE is tracked robustly by our model (Figure 6).The agreement between the modeled and observed RE is the best for the desert steppe and meadow steppe (R 2 > 0.9), followed by the alpine shrub meadow (R 2 = 0.872), the alpine meadow steppe (R 2 = 0.805), the alpine Kobresia meadow (R 2 = 0.778), the typical steppe (R 2 = 0.743), and the alpine swamp meadow (R 2 = 0.715).Ground covered by water in the alpine swamp meadow would affect the spectral reflectance and further disturb the model drivers (EVI and LSWI), and this may be the main reason for the relatively low performance in this vegetation type.
The seasonal variation of RE in China's grasslands could be described as a unimodal curve peaking in July or August.Across the seven vegetation types, the mean annual peak respiration rate ranges from 1.10 to 5.84 g C m −2 day −1 .The highest value occurring in the alpine meadow or the meadow steppe is approximately four times higher than that of the alpine meadow steppe or the desert steppe.There are also significant inter-annual fluctuations in the RE at some sites; for instance, the observed RE decreased remarkably as precipitation decreased in 2011 in SZWQ (Figure 6f), 2005 and 2006 in NM (Figure 6e), and 2009 in HBKO (Figure 6a).In addition, our model could accurately reproduce the inter-annual dynamics of RE due to the variation in precipitation across years.

Model Simulation for Seasonal and Inter-Annual Dynamics of RE at the Site Scale
Using the model parameters in Table 2, we simulate time-series RE of all 57 site-years and then further compared them to the flux tower-observed RE for each vegetation subtype in the alpine and temperate grasslands.In most cases, the temporal variation of the RE is tracked robustly by our model (Figure 6).The agreement between the modeled and observed RE is the best for the desert steppe and meadow steppe (R 2 > 0.9), followed by the alpine shrub meadow (R 2 = 0.872), the alpine meadow steppe (R 2 = 0.805), the alpine Kobresia meadow (R 2 = 0.778), the typical steppe (R 2 = 0.743), and the alpine swamp meadow (R 2 = 0.715).Ground covered by water in the alpine swamp meadow would affect the spectral reflectance and further disturb the model drivers (EVI and LSWI), and this may be the main reason for the relatively low performance in this vegetation type.

Model Simulation for Spatial Patterns of RE at the Regional Scale
Based on the parametrized RE model dependent on prescribed PFTs (alpine and temperate grasslands in Figure 1) and the spatial driver data of MODIS EVI, LST, and LSWI, the spatial distribution of the mean annual REs from 2001 to 2010 were obtained (Figure 7a).The predicted annual RE of China's northern grasslands is 258.18 ± 12.07 g C m −2 yr −1 on average.For the alpine grasslands on the Tibetan Plateau, the alpine shrub meadow has the largest RE value (434.69 ± 8.83 g C m −2 yr −1 ), followed by the alpine swamp meadow (413.82 ± 8.21 g C m −2 yr −1 ), the alpine Kobresia meadow (368.64 ± 7.52 g C m −2 yr −1 ), and the alpine steppe meadow (146.57± 3.21 g C m −2 yr −1 ), exhibiting a clear decreasing gradient from the southeast to the northwest with increasing elevation across the Tibetan Plateau.For the Inner Mongolian Plateau, RE also shows a clear decreasing gradient from northeast to southwest.Specifically, the highest value of RE occurs in the meadow steppe (470.27 ± 39.24 g C m −2 yr −1 ) in the northeast, and the lowest value of RE occurs in the desert steppe in the southwest (144.91 ± 21.75 g C m −2 yr −1 ).The mean annual RE in the typical steppe is 193.6 ± 30.84 g C m −2 yr −1 .
Reref appears to be a key component in this model in terms of quantifying the spatial variability of RE across sites.Nevertheless, the spatial patterns of Reref have rarely been investigated due to the lack of spatially explicit algorithms to allow the upscaling of Reref.Here, we took abiotic and biotic factors into consideration to parameterize Reref spatially and revealed that the mean annual Reref has a clear spatial heterogeneity among different vegetation types in northern China's grasslands (Figure 7b), ranging from 0.22 g C m −2 day −1 in the desert steppe to 2.25 g C m −2 day −1 in the meadow steppe.In addition, the spatial distribution of Reref, to a large extent, determines the spatial pattern of the mean annul RE (Figure 7a,b).The seasonal variation of RE in China's grasslands could be described as a unimodal curve peaking in July or August.Across the seven vegetation types, the mean annual peak respiration rate ranges from 1.10 to 5.84 g C m −2 day −1 .The highest value occurring in the alpine meadow or the meadow steppe is approximately four times higher than that of the alpine meadow steppe or the desert steppe.There are also significant inter-annual fluctuations in the RE at some sites; for instance, the observed RE decreased remarkably as precipitation decreased in 2011 in SZWQ (Figure 6f), 2005 and 2006 in NM (Figure 6e), and 2009 in HBKO (Figure 6a).In addition, our model could accurately reproduce the inter-annual dynamics of RE due to the variation in precipitation across years.

Model Simulation for Spatial Patterns of RE at the Regional Scale
Based on the parametrized RE model dependent on prescribed PFTs (alpine and temperate grasslands in Figure 1) and the spatial driver data of MODIS EVI, LST, and LSWI, the spatial distribution of the mean annual REs from 2001 to 2010 were obtained (Figure 7a).The predicted annual RE of China's northern grasslands is 258.18 ± 12.07 g C m −2 yr −1 on average.For the alpine grasslands on the Tibetan Plateau, the alpine shrub meadow has the largest RE value (434.69 ± 8.83 g C m −2 yr −1 ), followed by the alpine swamp meadow (413.82 ± 8.21 g C m −2 yr −1 ), the alpine Kobresia meadow (368.64 ± 7.52 g C m −2 yr −1 ), and the alpine steppe meadow (146.57± 3.21 g C m −2 yr −1 ), exhibiting a clear decreasing gradient from the southeast to the northwest with increasing elevation across the Tibetan Plateau.For the Inner Mongolian Plateau, RE also shows a clear decreasing gradient from northeast to southwest.Specifically, the highest value of RE occurs in the meadow steppe (470.27 ± 39.24 g C m −2 yr −1 ) in the northeast, and the lowest value of RE occurs in the desert

Biotic and Climatic Control over RE
Our analysis reveals that among biotic and climatic factors, plant productivity and moisture exert stronger influences than temperature on regulating the spatial pattern of RE.These two factors together can well explain the spatial patterns of RE, accounting for 54% and 58% of the variation in the temperate grasslands and alpine grasslands, respectively (Table 4).The important role of plant productivity and moisture in regulating the spatial pattern of RE is in good accordance with other studies across the world [20,46,59,60].Most of the variation in RE could be attributed to the difference in plant productivity among sites, with a relatively smaller proportion being further explained by moisture, but the moisture effect differs significantly between the two grasslands and plays a more crucial role in the water-deficit temperate grasslands than the alpine grasslands (Figure 3, Table 4).Although temperature has been noted to be positively correlated with RE across biomes at a global scale [40,41] and an important limitation on RE in alpine ecosystems, the inclusion of LSTmean does not significantly affect the explanation of Reref variation (Table 4).In the alpine grasslands, studies increasingly suggest that plant biomass is the direct determinant of RE distribution [61][62][63], while moisture changes, but not temperature, serve as another indirect determinant via regulating plant growth and distribution in the vast alpine grasslands (4400-4800 m) [64].The above two explanations describe how temperature is less important to governing the spatial pattern of RE in the alpine grasslands.However, the dependence of Reref on LSTmean appears to be stronger and negative in the temperate grasslands (Figure 3, Table 2).In fact, the relation of RE with temperature can be confounded by moisture in arid and semi-arid ecosystems [63,65].Specifically, the gradient of precipitation mainly defines the declining gradient of RE from the northeast to the southwest in the Inner Mongolian Plateau, which is exactly opposite to the gradient of temperature [66], thus leading to the apparently close and negative spatial correlation between temperature and RE.

Biotic and Climatic Control over RE
Our analysis reveals that among biotic and climatic factors, plant productivity and moisture exert stronger influences than temperature on regulating the spatial pattern of RE.These two factors together can well explain the spatial patterns of RE, accounting for 54% and 58% of the variation in the temperate grasslands and alpine grasslands, respectively (Table 4).The important role of plant productivity and moisture in regulating the spatial pattern of RE is in good accordance with other studies across the world [20,46,59,60].Most of the variation in RE could be attributed to the difference in plant productivity among sites, with a relatively smaller proportion being further explained by moisture, but the moisture effect differs significantly between the two grasslands and plays a more crucial role in the water-deficit temperate grasslands than the alpine grasslands (Figure 3, Table 4).Although temperature has been noted to be positively correlated with RE across biomes at a global scale [40,41] and an important limitation on RE in alpine ecosystems, the inclusion of LST mean does not significantly affect the explanation of Reref variation (Table 4).In the alpine grasslands, studies increasingly suggest that plant biomass is the direct determinant of RE distribution [61][62][63], while moisture changes, but not temperature, serve as another indirect determinant via regulating plant growth and distribution in the vast alpine grasslands (4400-4800 m) [64].The above two explanations describe how temperature is less important to governing the spatial pattern of RE in the alpine grasslands.However, the dependence of Reref on LST mean appears to be stronger and negative in the temperate grasslands (Figure 3, Table 2).In fact, the relation of RE with temperature can be confounded by moisture in arid and semi-arid ecosystems [63,65].Specifically, the gradient of precipitation mainly defines the declining gradient of RE from the northeast to the southwest in the Inner Mongolian Plateau, which is exactly opposite to the gradient of temperature [66], thus leading to the apparently close and negative spatial correlation between temperature and RE.
In contrast, temperature becomes the major controlling factor in terms of the temporal pattern of RE (Figure 4).The significant exponential dependence of seasonal RE on temperature that was observed in these two regions is consistent with many other grassland ecosystems [67,68].However, there are obvious differences in the temperature control over RE among different climate zones.For the alpine grasslands, temperature explained 75.2% of the seasonal variability of RE and was the strongest controlling factor.However, the magnitude of the temporal dynamics of RE due to temperature decreased to 68.7% for the temperate grasslands.This result is in accordance with the argument that the seasonal variation of RE may depend much less on temperature when moisture become a limiting factor [23,65,69].As expected, moisture exerts a strong impact on the seasonal changes of RE in the arid and semi-arid temperate grasslands, almost to the same degree as temperature (Figure 4).While the dependence of RE on LSWI is relatively weak, productivity serves as a relatively stronger controlling factor (Figure 4).Although some studies support the close correlation between the seasonal dynamics of RE and plant growth [68,70], this correlation still contradicts the more universal argument that the abiotic factors (temperature and moisture), rather than biotic factors (plant growth), mainly regulate the seasonal patterns of RE [67].The most significant cause of this strong dependency on production in alpine grasslands is the relatively high root biomass density in this grassland community [71]; thus, autotrophic respiration, which is closely related to plant biomass, contributes to a relatively large proportion of the total RE [62].

Model Evaluation
The model demonstrates that flux-tower observed RE can be estimated with an applicable degree of accuracy across grasslands using only remote sensing data.We first evaluate the precision of our model at a site scale by comparing it to other recent semi-empirical approaches, e.g., the two original developed-based models, T&P [22] and RECO [27] that achieved a mean site-level accuracy of R 2 = 0.679 with RMSE =1.227 g C m −2 day −1 and R 2 = 0.625 with RMSE = 0.969 g C m −2 day −1 , respectively, for these 16 grassland sites, while our model performed better with a higher degree of explanation and a lower bias (Table 3), verifying the model improvements with better representation of Reref and further incorporation of the moisture effect.An independent regional RE product from a popular process-based model CLM4CN [56] was employed to further evaluate the model performance at a regional scale.The CLM4CN presented a mean annual RE value (309.18 ± 5.15) within an order of magnitude of this study, but it still resulted in an overestimation.The general spatial pattern of RE appeared similar but there were still distinct differences, especially in the water-deficit area.For example, in the southwest of the desert grasslands which is driest, the CLM4CN failed to produce the lowest values.Moreover, the CLM4CN unreasonably produced a higher RE in the alpine meadow steppe than the surrounding alpine meadow (Figure 7a,c).These distinct differences are likely associated with the different representation of Reref in the models, which has been proven to largely determine the spatial pattern of RE (Figure 7a,b).Reref was set as the same constant for all biomes in the CLM4CN [72], as in the most process-based models [73][74][75].However, more and more studies demonstrate that a large error would be introduced if a constant Reref is used [76].In fact, Reref varies remarkably in space [9,22,60].Therefore, the better representation of Reref makes our model more sensible for capturing the regional characteristics of RE in grassland ecosystems and arid and semi-arid areas, in particular.
The consideration of moisture effects on both the spatial and temporal variabilities of RE is another major factor contributing to the better performance of our model.For the spatial variability of RE (Reref ), even though it has been gradually confirmed to vary systematically with the spatial variation of productivity (GPP or LAI) and temperature [8,[20][21][22]27,50], to our knowledge, this process has not yet been linked to moisture change in satellite-based RE studies at large scales and with a high resolution.Our findings reveal that Reref has a significant linear relation with LSWI mean in both the Tibetan and Inner Mongolian grasslands (Figure 3, p < 0.01), which is also supported by another recent study in northern China [24].Considerable improvement was found in the explanation of Reref (12%, Table 4) in the temperate grasslands by the inclusion of the moisture effect, while the improvement is relatively small (4%) in the alpine grasslands due to the colinearity of plant biomass and soil moisture [62].
Moisture also exerts significant influence on the seasonal dynamics of RE in northern China's grasslands (Figure 4), which further explained a large proportion of RE in the arid and semi-arid vegetation types (Table 5).Obvious gaps would be found in the performance between the model without fw and the model with fw for the arid and semi-arid ecosystems.The simulated RE with fw has a higher correlation coefficient with observed values (R 2 increased by 12-24%) and a smaller bias from the observed values (RMSE reduced by 13-41%) in the desert steppe, typical steppe, and alpine meadow steppe.Although various equations (e.g., linear, quadratic, parabolic, logarithmic, exponential, and hyperbolic form) have been proposed to describe the response process of RE to moisture [69, 77,78], our findings reveal that the Michaelis-Menten equation is the most applicable form for northern China's grasslands, as the moisture content is hardly able to reach beyond field capacity and toward saturation in these water-limited grasslands [62,66].Therefore, the process of the steep decline in RE at very high moisture levels could be ignored, while the remaining two phases, i.e., the strong increase at low moisture levels and the plateau of RE in response to a broad range of near optimum water contents, agree well with the typical Michaelis-Menten-shaped curve.

Model Applications and Limitations
Our model was developed with a robust biological basis with long-term and inter-sites observations, and it performed well in estimating RE for northern China's grasslands solely driven by MODIS products.The main advances introduced by this model are the better representation of the Reref process and the comprehensive consideration of moisture effect combined with temperature and productivity effects on RE.The derived parameters reported in Table 2 may be considered as an optimized parameterization for the application of the model at a regional scale.Our model may facilitate the simulation of RE at larger spatial scales with a high temporal resolution and could also be combined with the satellite-based GPP models or MODIS-retrieved GPP products to conduct a regional-scale NEE simulation.
However, the further application of our model at a larger scale probably has several limitations: (1) The spatial variability (Reref ) is explained to a relatively lower degree compared to the temporal variability (f (T,W,P)).Since RE is also influenced by a number of other biotic and abiotic factors, such as age of the ecosystem [57], nutrient availability [79], and acclimation and disturbance effects, we expect that more comprehensive information could be incorporated in the model to better represent the variations of Reref ; (2) Although the satellite-derived VIs can to some extent reflect the impacts of GPP on respiration [80,81] and have also been widely used to explain the seasonal variation in Rs as VIs are related to soil carbon through litter input [16,82,83], they are still difficult to mirror in terms of the influence of the belowground organic matter on RE; (3) The performance of our model for other vegetation types or in other regions needs to be validated in further study.Without further testing, the model is strictly applicable to similar ecosystems (i.e., temperate and alpine grassland ecosystems).Along with the development of satellite algorithms for soil moisture, plant biomass, nutrient availability, human activity, and so on, we will optimize the model to improve its accuracy.

Conclusions and Implications
A statistical-mechanistic model driven by remote-sensing data was formulated for simulating the spatial and temporal variability of RE.The effect of plant productivity, temperature, and moisture were comprehensively considered using the publicly available MODIS products EVI, LST, and LSWI, thereby yielding a better performance for the RE simulation across seven vegetation types in the Tibetan alpine grasslands and Inner Mongolian temperate grasslands.The temporal variations as well as the spatial patterns of RE in northern China's grasslands were well reproduced by the model.The inclusion of LSWI in the Reref and the seasonal dynamics of RE significantly improved the model performance in arid and semi-arid ecosystems.Our results suggest that a good representation of the spatial process and moisture effect on RE should be considered in the next-generation satellite-based RE models.
Biotic and climatic control over RE is different for temporal and spatial processes.Plant productivity and moisture mainly contribute to the spatial variation of RE in northern China's grasslands.The influence of moisture on the spatial pattern of RE is stronger in the temperate grasslands than in the alpine grasslands.Temperature plays a minor role in regulating the spatial pattern of RE.In contrast, temperature tends to be more important in controlling the seasonal patterns of RE in northern China's grasslands, whereas in the temperate grasslands of Inner Mongolia, moisture exerts strong impact on the seasonal variation of RE, almost equally important as temperature.These conclusions drawn from the Tibetan and Inner Mongolian grasslands, on behalf of the alpine grasslands and Eurasian temperate grasslands across the word, can provide valuable information for large-scale estimates of RE and better understanding the response of RE to climate change in grassland ecosystems.

Figure 2 .
Figure 2. Structure of the ecosystem respiration (RE) model for grasslands.

Figure 2 .
Figure 2. Structure of the ecosystem respiration (RE) model for grasslands.

Figure 4 .
Figure 4. f (T,P,W) response to nighttime LST (LST n ), standardized EVI (EVI std ), and LSWI for the Tibetan alpine grasslands and the Inner Mongolian temperate grasslands.

Remote Sens. 2018, 10 , 149 11 of 20 Figure 5 .
Figure 5.Comparison between the observed and predicted ecosystem respiration (RE) in the Tibetan alpine grasslands (left) and the Inner Mongolian temperate grasslands (right).

1 Figure 5 .
Figure 5.Comparison between the observed and predicted ecosystem respiration (RE) in the Tibetan alpine grasslands (left) and the Inner Mongolian temperate grasslands (right).

Figure 6 .
Figure 6.Time-series plots for observed (black square dot) and modeled (red round dot) ecosystem respiration (RE) at all sites.The text in the abscissa represents the abbreviation of the site name, and the following number represents the year.

Figure 6 .
Figure 6.Time-series plots for observed (black square dot) and modeled (red round dot) ecosystem respiration (RE) at all sites.The text in the abscissa represents the abbreviation of the site name, and the following number represents the year.

Figure 7 .
Figure 7. Predicted spatial distribution of the mean annual ecosystem respiration (RE).(a) and the reference respiration (Reref); (b) by our model and mean annual RE by CLM4CN; (c) in northern China's grasslands.

Figure 7 .
Figure 7. Predicted spatial distribution of the mean annual ecosystem respiration (RE).(a) and the reference respiration (Reref ); (b) by our model and mean annual RE by CLM4CN; (c) in northern China's grasslands.

Table 1 .
Main characteristics of the 16 flux-tower sites over China's grasslands.

Table 2 .
Parameter estimation of Reref and f (T,P,W) for two grassland types in two regions.
* Among the parameters for f (T,P,W), the values inside the parentheses are for the non-growing season, while the others are for the growing season.

Table 3 .
Joint-sites estimation and leave-one-out cross validation for each vegetation type.

Table 4 .
Relationship between Reref and the combined effects of EVImean, LSTmean, and LSWImean.

Table 4 .
Relationship between Reref and the combined effects of EVI mean , LST mean , and LSWI mean .

Table 5 .
Comparison of model performance between the model without fw and the model with fw for each vegetation type.