Spatiotemporal Changes and Driver Analysis of Ecosystem Respiration in the Tibetan and Inner Mongolian Grasslands

: Ecosystem respiration (RE) plays a critical role in terrestrial carbon cycles, and quantiﬁ-cation of RE is important for understanding the interaction between climate change and carbon dynamics. We used a multi-level attention network, Geoman, to identify the relative importance of environmental factors and to simulate spatiotemporal changes in RE in northern China’s grasslands during 2001–2015, based on 18 ﬂux sites and multi-source spatial data. Results indicate that Geoman performed well (R 2 = 0.87, RMSE = 0.39 g C m − 2 d − 1 , MAE = 0.28 g C m − 2 d − 1 ), and that grassland type and soil texture are the two most important environmental variables for RE estimation. RE in alpine grasslands showed a decreasing gradient from southeast to northwest, and that of temperate grasslands showed a decreasing gradient from northeast to southwest. This can be explained by the enhanced vegetation index (EVI), and soil factors including soil organic carbon density and soil texture. RE in northern China’s grasslands showed a signiﬁcant increase (1.81 g C m − 2 yr − 1 ) during 2001–2015. The increase rate of RE in alpine grassland (2.36 g C m − 2 yr − 1 ) was greater than that in temperate grassland (1.28 g C m − 2 yr − 1 ). Temperature and EVI contributed to the interannual change of RE in alpine grassland, and precipitation and EVI were the main contributors in temperate grassland. This study provides a key reference for the application of advanced deep learning models in carbon cycle simulation, to reduce uncertainties and improve understanding of the effects of biotic and climatic factors on spatiotemporal changes in RE.


Introduction
Ecosystem respiration (RE) is the second-largest terrestrial carbon flux after photosynthesis, and plays an important role in the interaction between climate change and carbon dynamics [1,2]. Minor fluctuations in RE caused by natural or human activity can result in significant changes in atmospheric CO 2 concentrations [2][3][4][5]. Although studies have described the global RE pattern and associated temporal changes [6][7][8], there remains considerable uncertainty around the spatiotemporal changes in RE at the regional and and biome scales [9], hindering the identification of carbon sinks and source functions in key areas.
Efforts have been made to understand how climate change affects spatiotemporal changes in RE, but complex interactions among the physical, chemical, and biological processes of RE present a considerable challenge for researchers [10,11]. Seasonal changes in grassland RE can be attributed to temperature (Ta) and moisture [2,12,13]. Spatial changes in respiration can be attributed to vegetation status, plant productivity, or moisture [12,14,15]. However, there have been few analyses of long-term interannual changes in RE. Most of these prior studies have focused on climatic and vegetation factors. Soil factors remain relatively unexplored, including soil organic carbon (SOC) and soil texture (SoilTex). This is because their functions remain poorly understood in the context of complex respiration processes [16][17][18]. However, soil factors are important for carbon flux estimation, and soil respiration accounts for a high proportion of ecosystem respiration [19][20][21][22][23][24]. Therefore, the important influence of soil factors in the RE estimation model cannot be excluded.
To accurately quantify RE, many empirical models [2,14,[25][26][27], semi-empirical models [1,28], and process models [16,29] have been developed using data at different temporal and spatial scales. Machine learning (ML) models show promise for solving complex data patterns, owing to their powerful nonlinear regression and classification abilities [30]. ML has been successfully applied in many fields of ecology, including carbon-water flux [2,[31][32][33][34][35], energy flux [6,36,37], and vegetation state [38][39][40]. However, the interpretation of the modeling results is often difficult, as the impact of inputs on the output cannot easily be understood through networks in comparison with process-based models [30,[41][42][43]. In recent years, the attention mechanism has become one of the main focuses in ML research. Its core principle is to assign greater weight to important information to achieve the purpose of selecting the key information, which improves the interpretability of the results [42,44,45]. Combined with powerful neural network models, such as long shortterm memory networks (LSTM), the attention mechanism has achieved great success in time series prediction [44,46,47]. ML models with a high prediction accuracy provide promising tools for monitoring regional RE.
Grasslands play an important role in the terrestrial ecosystem carbon cycle, especially grassland ecosystems in alpine regions and in arid and semi-arid regions, where fragile ecological environments are more likely to be affected by climate change [14,48]. Grassland is the dominant landscape in China, accounting for 40% of the national land area [2], and is mainly distributed in the Tibetan Plateau (TP) and the Inner Mongolian Plateau (IM) (Figure 1). China's grasslands account for approximately 10% of the global grassland area and form 9-16% of the global grassland carbon pool [49]. In this study, we applied the geo-sensory multi-level attention network (Geoman) model developed by Liang et al. [44] to estimate the RE in northern China's grasslands from 2001-2015, at a spatial resolution of 1 km. We aimed to evaluate the performance of the deep learning model based on the attention mechanism, for estimating the RE of grassland ecosystems in northern China. Based on the effects of multiple environmental factors, including soil, meteorology, and biology, we explored long-term spatiotemporal changes of RE in northern China's grasslands, and analyzed the driving mechanisms and effects of these environmental factors on spatial patterns and interannual and seasonal dynamics of RE.

Study Area
Located in an alpine climate zone, the TP is characterized by lengthy sunshine duration, strong solar radiation, low temperature, and low precipitation (PRE). The average altitude of the TP is more than 4000 m, the mean annual air temperature is in the range of −5.75 to 2.57 • C, and the mean annual precipitation varies from 200 to 600 mm. The terrain of the TP slopes from northwest to southeast. Alpine grasslands cover more than 60% of the surface of the TP and are a unique grassland ecosystem among the world's alpine regions [12]. IM is in an arid and semi-arid zone, and the average altitude is approximately 1000 m, with a terrain that is high in the south and low in the north. The mean annual air temperature in the IM ranges between 3 and 6 • C, and the mean annual precipitation ranges from 200 to 350 mm. Temperate grasslands are the main grassland type in the IM, and are a typical vegetation type (VT) in a temperate continental climate.

Flux and Meteorological Data
Flux and meteorological data were collected from the ChinaFLUX [51], Coordinated Observations and Integrated Research over Arid and Semi-arid China (COIRAS) [52], and the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) [53]. A total of 52 site-year observation data points collected from 18 flux sites during 2003-2014 were used in this study ( Figure 1, Table 1). Among them, nine sites were distributed in the TP alpine grasslands and included 34 site-year observation data points. Another nine sites were distributed in the IM temperate grasslands and included 18 site-year observation data points. These flux sites represent the most extensive grassland ecosystem types, with a wide range of spatial, ecological, and climatic conditions. The carbon flux data were processed through triple coordinate rotation, the Webb-Pearman-Leuning (WPL) correction, and abnormal data rejection [51]. To match the eight-day compositing interval of the MODIS products, the half-hour observed RE, Ta and photosynthetically active radiation (PAR) data were averaged and the PRE was summed over eight days. At the site scale, these meteorological data were synchronously observed using the eddy covariance system, and were gap-filled following the method proposed by Schwalm et al. [54]. Finally, we obtained a total of 1923 items of site observation data with an eight-day temporal resolution. At the wider spatial scale, Ta, PRE, and PAR data from 2001 to 2015 with 1 km spatial resolution were produced using ANUSPLIN software to spatialize the observation data from Chinese meteorological stations. The regional-scale Ta, PRE, and PAR data were then cross-validated based on the site observation data. The validation results showed that the interpolation products were consistent with the site observation data [55,56]. In line with the RE data, the regional Ta and PAR data were averaged and the PRE was summed for eight days.
where ρ N IR , ρ red and ρ blue are the surface reflectance for the near-infrared, red, and blue bands, respectively. The coefficients G, C 1 , C 2 and L used in the EVI algorithm were 2.5, 6, 7.5, and 1, respectively [57]. The original spatial resolution of the EVI was 500 m; it was resampled to 1 km for ease of computation. To further reduce the effects of cloud and to capture the seasonality of the EVI, a data smoothing tool providing a double logistic curve fit in the TIMESAT software was used to smooth the original time series [58].
Digital elevation model (DEM) data was obtained from the Shuttle Radar Topographic Mission (https://srtm.csi.cgiar.org/srtmdata/, accessed on 10 November 2021) and the original data were resampled with a spatial resolution of 90 m to 1 km.

Soil Data
The SoilTex of the surface soil layer (0-30 cm) was retrieved from the Harmonized World Soil Database (http://webarchive.iiasa.ac.at/Research/LUC/External-World-soildatabase/, accessed on 16 November 2021), including the clay content (wt%), silt content (wt%), and sand content (wt%). The soil organic carbon density (SOCD, kg C m −2 ) was calculated using the organic carbon content (wt%), gravel content (vol%), layer thickness (m), and bulk density (kg m −3 ), according to the method described by Carvalhais et al. [59]. Due to limited data availability, we only obtained one issue of SOCD and SoilTex data. We used the mean values of the SoilTex, SOCD, and EVI subsets of 3 × 3 km pixels centered on each flux tower, to better represent the eddy covariance footprint area and to reduce the effects of geolocation errors [34].

Deep Learning Model: Geoman
In this study, we used a deep learning model including a multi-level attention mechanism, known as Geoman. It was originally developed to perform geographic sensor time-series predictions, such as for air quality and water quality [44]. The structure was modified from the original Geoman framework and is shown in Figure 2. Geoman combines the advantages of the attention mechanism and LSTM by distinguishing the importance of multiple time-series data to preserve key spatiotemporal information. The attention module in Geoman consists of spatial and temporal attention. The spatial attention includes two sub-modules, namely global spatial attention that expresses the inter-site correlations, and local spatial attention that expresses the intra-site correlations. The temporal attention in the decoder stage captures the time dependencies. LSTM learns interactions at multiple timescales and can be extended to encoder and decoder modules [60]. Limited by the availability of ecological data, and also by differences in the spatial density between ecological observation sites and other geographic monitoring sites, including loop detector systems and urban groundwater monitoring systems, we only used local spatial attention to describe the dynamic correlation between each local feature and the target series. For simplicity, in this paper the term 'Geoman' refers to the model without global spatial attention. The local spatial attention mechanism can be expressed simply by Equations (2) and (3).
where [; ] represents concentration operation. v l , b l ∈ R T , W l ∈ R T×2m and U l ∈ R T×T are learnable parameters and are determined by the back-propagation algorithm during the training process. The superscript T is the time window of input data. h t−1 ∈ R m and s t−1 ∈ R m are the hidden state and cell state at time t − 1, constituting the short-term memory and long-term memory in the LSTM encoder, respectively. The superscript m is the number of hidden neurons in LSTM. N l denotes the number of feature variables. x k is the k-th local feature, namely the environmental variables (VT, Ta, PAR, EVI, PRE, SOCD, or SoilTex) of the flux sites. e k t is the correlation score between the k-th local feature and the target series (i.e., RE). Its size is determined by the local input feature x k and historical state information (h t−1 and s t−1 ). α k t is the attention weight calculated by normalizing the score of the k-th feature using a softmax function, and represents the relative importance of the k-th feature to the target series. After training, α k t was obtained from the trained Geoman prediction process. The external factor fusion module can fuse other relevant features affecting the RE, such as elevation and weather forecasts. The weather forecasts used in this study refer to Ta and PRE at the time stamp of the predicted RE, and provide time-related factors to further improve the prediction accuracy.

Model Training and Evaluation
The hyperparameters of Geoman were optimized and determined using the grid search method [61]. All data from the 18 sites were fed into the model for training, without considering the interactions between sites. The performance of Geoman was evaluated using a sample-based 10-fold cross-validation strategy in which the observation data were randomly assigned 10 folds. Each fold contained approximately 10% of the data and was used consecutively for model training, validation, and testing to ensure independence of the three datasets [60]. According to this partition principle, the entire data were divided into non-overlapping training, validation, and test data at a ratio of 8:1:1. The cross-validation process was repeated 10 times. The estimated RE from all 10 folds were compared with the observed RE.
Three frequently used statistical indicators, i.e., the coefficient of determination (R 2 ), root-mean-square error (RMSE), and mean absolute error (MAE) were selected as evaluation metrics to analyze the performance of the model. These metrics were calculated as follows: 1 Figure 2. The structure of Geoman adapted from Liang et al. [44]. LSTM: long short-term memory networks. Attn: attention. Concat: concatenation layer.ŷ t : the predicted value at time t. C t : the context vectors at time t. h 0 : the initial state of encoder. VT: vegetation type. Ta: temperature. PAR: photosynthetically active radiation. EVI: enhanced vegetation index. PRE: precipitation. SOCD: soil organic carbon density. SoilTex: soil texture.

Model Training and Evaluation
The hyperparameters of Geoman were optimized and determined using the grid search method [61]. All data from the 18 sites were fed into the model for training, without considering the interactions between sites. The performance of Geoman was evaluated using a sample-based 10-fold cross-validation strategy in which the observation data were randomly assigned 10 folds. Each fold contained approximately 10% of the data and was used consecutively for model training, validation, and testing to ensure independence of the three datasets [60]. According to this partition principle, the entire data were divided into non-overlapping training, validation, and test data at a ratio of 8:1:1. The cross-validation process was repeated 10 times. The estimated RE from all 10 folds were compared with the observed RE.
Three frequently used statistical indicators, i.e., the coefficient of determination (R 2 ), root-mean-square error (RMSE), and mean absolute error (MAE) were selected as evaluation metrics to analyze the performance of the model. These metrics were calculated as follows: where y i and y i are the observed and estimated RE, respectively. y and y are the means of y i and y i , and n is the number of observed samples. When all parameters were optimized, the trained Geoman model was applied to predict eight-day RE at all spatial grids, using environmental data with l km spatial resolution and eight-day temporal resolution. The grid numbers for each grassland subtype are shown in Table 2. The DEM, SOCD, and SoilTex data were set unchanged at all time steps due to their limited availability.

. Model Evaluation
According to the results of the 10-fold cross-validation, the predicted RE was in line with the observed RE ( Figure 3). Most of the scattered points occurred approximately at the 1: where and are the observed and estimated RE, respectively. and are the means of and , and n is the number of observed samples. When all parameters were optimized, the trained Geoman model was applied to predict eight-day RE at all spatial grids, using environmental data with l km spatial resolution and eight-day temporal resolution. The grid numbers for each grassland subtype are shown in Table 2. The DEM, SOCD, and SoilTex data were set unchanged at all time steps due to their limited availability.

Model Evaluation
According to the results of the 10-fold cross-validation, the predicted RE was in line with the observed RE ( Figure 3). Most of the scattered points occurred approximately at the 1:1 line (R 2 = 0.87, RMSE = 0.39 g C m −2 d −1 , MAE = 0.28 g C m −2 d −1 ), indicating that Geoman can accurately estimate the RE value for northern China's grasslands. We further evaluated the model simulation for the seasonal dynamics of RE at the site scale. Results showed that Geoman can accurately capture the seasonal dynamics of RE in the seven grassland subtypes ( Figure S1). The best agreement between the predicted and observed RE was for typical steppe and meadow steppe (R 2 = 0.97 and 0.95, respectively), followed by the alpine shrub meadow (R 2 = 0.91), alpine Kobresia meadow (R 2 = 0.89), desert steppe (R 2 = 0.88), alpine swamp meadow (R 2 = 0.86), and alpine meadow steppe (R 2 = 0.86).
The seasonal variation in the RE at all sites formed a unimodal curve, peaking in July or August. The model continued to perform well during the high value period (from July to August), the growing season (from May to September), and the low value period (from December to January) in the non-growing season. Across the seven grassland subtypes, the alpine Kobresia meadow, the alpine shrub meadow, and the meadow steppe had higher RE peaks, reaching 4-6 g C m −2 d −1 (Figure S1a,b,f), while the alpine meadow steppe and desert steppe had relatively low RE peaks of 0.5-2 g C m −2 d −1 (Figure S1d,f). We further evaluated the model simulation for the seasonal dynamics of RE at the site scale. Results showed that Geoman can accurately capture the seasonal dynamics of RE in the seven grassland subtypes ( Figure S1). The best agreement between the predicted and observed RE was for typical steppe and meadow steppe (R 2 = 0.97 and 0.95, respectively), followed by the alpine shrub meadow (R 2 = 0.91), alpine Kobresia meadow (R 2 = 0.89), desert steppe (R 2 = 0.88), alpine swamp meadow (R 2 = 0.86), and alpine meadow steppe (R 2 = 0.86).
The seasonal variation in the RE at all sites formed a unimodal curve, peaking in July or August. The model continued to perform well during the high value period (from July to August), the growing season (from May to September), and the low value period (from December to January) in the non-growing season. Across the seven grassland subtypes, the alpine Kobresia meadow, the alpine shrub meadow, and the meadow steppe had higher RE peaks, reaching 4-6 g C m −2 d −1 (Figure S1a,b,f), while the alpine meadow steppe and desert steppe had relatively low RE peaks of 0.5-2 g C m −2 d −1 (Figure S1d,f).

Relative Importance of Environmental Factors
In this study, we quantified the relative importance of each environmental variable for predicting the RE, using the dynamic attention weights. We analyzed the differences between seven grassland subtypes during the growing season and non-growing season ( Figure 4, Table 3). Results demonstrate that vegetation type and soil texture are the two most important environmental factors for RE estimation in northern China's grasslands. The importance of other environmental variables also represents spatiotemporal characteristics with ecological significance. In the growing season, temperature was found to be more important than precipitation in alpine grasslands, whereas the opposite was true for temperate grasslands. This indicates that restrictive climatic conditions have a considerable effect on accurate estimation of RE. In the non-growing season, SOCD was more important than climatic variables for RE estimation. However, temperature was more important than SOCD in alpine meadows, which may be associated with the low temperature limit of the RE in the area. In this study, we quantified the relative importance of each environmental variable for predicting the RE, using the dynamic attention weights. We analyzed the differences between seven grassland subtypes during the growing season and non-growing season ( Figure 4, Table 3). Results demonstrate that vegetation type and soil texture are the two most important environmental factors for RE estimation in northern China's grasslands. The importance of other environmental variables also represents spatiotemporal characteristics with ecological significance. In the growing season, temperature was found to be more important than precipitation in alpine grasslands, whereas the opposite was true for temperate grasslands. This indicates that restrictive climatic conditions have a considerable effect on accurate estimation of RE. In the non-growing season, SOCD was more important than climatic variables for RE estimation. However, temperature was more important than SOCD in alpine meadows, which may be associated with the low temperature limit of the RE in the area.  Table 3. Relative importance statistics for the environmental variables. All the notations in the grassland subtypes column are consistent with those in Figure 4.  Evaluation of the relative importance of the environmental variables indicates that the importance of SOCD in the non-growing season exceeded that of the climatic factors. Therefore, we further evaluated the accuracy of the RE estimation for the low value period from December to January in the non-growing season. It was found that adding SOCD to the model could improve interpretability by 44% compared with removing this variable ( Figure S2a,b). When the soil texture was removed from the model, the results appeared to have been overestimated during the low value period in the non-growing season ( Figure  S2c). Therefore, it is necessary to integrate soil texture into RE models because this can to a certain extent solve the problem of overestimation of low RE values. Table 3. Relative importance statistics for the environmental variables. All the notations in the grassland subtypes column are consistent with those in Figure 4.  (Figure 5a). Regarding the alpine grasslands on the TP, the mean annual RE was 467.44 g C m −2 yr −1 , exhibiting a clear decreasing gradient from southeast to northwest (Figure 5b). The alpine shrub meadow had the largest RE (573.64 g C m −2 yr −1 ), followed by the alpine Kobresia meadow (515.28 g C m −2 yr −1 ), the alpine swamp meadow (380.94 g C m −2 yr −1 ), and the alpine meadow steppe (129.29 g C m −2 yr −1 ). Regarding the temperate grasslands in IM, the mean annual RE was 270.43 g C m −2 yr −1 , exhibiting a clear decreasing gradient from northeast to southwest (Figure 5b). The highest RE value occurred in the meadow steppe (534.28 g C m −2 yr −1 ) in the northeast, and the lowest value for RE was recorded in the desert steppe in the southwest (88.15 g C m −2 yr −1 ). The mean annual RE in the typical steppe was 283.15 g C m −2 yr −1 . Statistics indicated that the RE for the alpine Kobresia meadow and typical steppe accounted for more than 60% of the total RE in northern China's grasslands (Figure 5c). Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 24

Temporal Variation for RE in Northern China's Grasslands
Seasonal variation in RE at the regional scale formed a unimodal curve peaking during summer (Figure 6a), which is consistent with the site-scale verification results ( Figure  S1). The RE for alpine grasslands in the TP was always higher than that of the temperate grasslands in the IM. The RE peak value for the TP was in the range of 2.5-3.0 g C m −2 day −1 , and its trough value was around 0.5 g C m −2 day −1 . The IM had a lower peak value of 2.0-2.5 g C m −2 day −1 , and its trough value was significantly lower than 0.5 g C m −2 day −1 . We analyzed seasonal dynamics for the mean 15-year seasonal RE in the seven grassland subtypes and found that vegetation heterogeneity had a significant impact on RE ( Figure  6b). The seasonal change curves for RE in the alpine meadow steppe and desert steppe were relatively flat, and the differences in RE between the growing season and non-growing season were relatively small, which may have resulted from hydrothermal condition limitations. Compared with these two grassland subtypes, the other grassland subtypes had a higher RE peak value during the growing season.

Temporal Variation for RE in Northern China's Grasslands
Seasonal variation in RE at the regional scale formed a unimodal curve peaking during summer (Figure 6a), which is consistent with the site-scale verification results ( Figure S1). The RE for alpine grasslands in the TP was always higher than that of the temperate grasslands in the IM. The RE peak value for the TP was in the range of 2.5-3.0 g C m −2 day −1 , and its trough value was around 0.5 g C m −2 day −1 . The IM had a lower peak value of 2.0-2.5 g C m −2 day −1 , and its trough value was significantly lower than 0.5 g C m −2 day −1 . We analyzed seasonal dynamics for the mean 15-year seasonal RE in the seven grassland subtypes and found that vegetation heterogeneity had a significant impact on RE (Figure 6b). The seasonal change curves for RE in the alpine meadow steppe and desert steppe were relatively flat, and the differences in RE between the growing season and non-growing season were relatively small, which may have resulted from hydrothermal condition limitations. Compared with these two grassland subtypes, the other grassland subtypes had a higher RE peak value during the growing season.
From 2001 to 2015, the RE of northern China's grasslands increased at a rate of 1.81 g C m −2 yr −1 (Figure 6c). The alpine grassland RE in TP increased at a rate of 2.36 g C m −2 yr −1 , which was faster than the increase rate of 1.28 g C m −2 yr −1 in the IM temperate grasslands (Figure 6c). We further analyzed the increase rate of RE in the seven grassland subtypes from 2001 to 2015 and found that the maximum rate of increase occurred in the alpine shrub meadow, while the minimum increase rate was in the desert steppe (Figure 6d). The rate of increase for RE in the TP is greater than that of the IM, which shows that the TP has a stronger feedback response under the background of global climate change. Remote Sens. 2022, 14, x FOR PEER REVIEW 11 of 24  (Figure 6c). The alpine grassland RE in TP increased at a rate of 2.36 g C m −2 yr −1 , which was faster than the increase rate of 1.28 g C m −2 yr −1 in the IM temperate grasslands (Figure 6c). We further analyzed the increase rate of RE in the seven grassland subtypes from 2001 to 2015 and found that the maximum rate of increase occurred in the alpine shrub meadow, while the minimum increase rate was in the desert steppe (Figure 6d). The rate of increase for RE in the TP is greater than that of the IM, which shows that the TP has a stronger feedback response under the background of global climate change.

Environment Effects on RE Seasonal Dynamics
The relationship between RE and Ta in northern China's grasslands at an eight-day timescale showed an explicit exponential form (R 2 = 0.97; Figure 7a). The sensitivity of RE to temperature in the TP alpine region was greater than that of the IM temperate region. This temperature-sensitive heterogeneity can be attributed to the regional climatic conditions and the vegetation types. There was a clear positive correlation between RE and PRE throughout northern China's grasslands, and the R 2 between the two was 0.90 ( Figure 7b). We also found that the RE reached a certain saturation point with the increase in PRE. When the PRE reacheds a specific threshold, the RE no longer increased and even tended to decline (Figure 7b). The saturation thresholds for RE to PRE differed in different climatic regions. In the temperate grasslands, RE began to decrease after the PRE exceeded 35 mm, while PRE at this intensity retained a continuous positive stimulus for RE in alpine grasslands. This may be associated with the different moisture requirements of the

Environment Effects on RE Seasonal Dynamics
The relationship between RE and Ta in northern China's grasslands at an eight-day timescale showed an explicit exponential form (R 2 = 0.97; Figure 7a). The sensitivity of RE to temperature in the TP alpine region was greater than that of the IM temperate region. This temperature-sensitive heterogeneity can be attributed to the regional climatic conditions and the vegetation types. There was a clear positive correlation between RE and PRE throughout northern China's grasslands, and the R 2 between the two was 0.90 ( Figure 7b). We also found that the RE reached a certain saturation point with the increase in PRE. When the PRE reacheds a specific threshold, the RE no longer increased and even tended to decline (Figure 7b). The saturation thresholds for RE to PRE differed in different climatic regions. In the temperate grasslands, RE began to decrease after the PRE exceeded 35 mm, while PRE at this intensity retained a continuous positive stimulus for RE in alpine grasslands. This may be associated with the different moisture requirements of the different grassland types. A significant linear positive correlation (R 2 = 0.93) between the RE and the EVI was observed in the northern grasslands of China (Figure 7c). Similar to Ta, PAR and RE showed an exponential positive correlation (R 2 = 0.69, Figure 7d). However, under the same radiation intensity, the sensitivity of RE to PAR was higher in the alpine grasslands than in the temperate grasslands. different grassland types. A significant linear positive correlation (R 2 = 0.93) between the RE and the EVI was observed in the northern grasslands of China (Figure 7c). Similar to Ta, PAR and RE showed an exponential positive correlation (R 2 = 0.69, Figure 7d). However, under the same radiation intensity, the sensitivity of RE to PAR was higher in the alpine grasslands than in the temperate grasslands. We further analyzed the relationship between RE and environmental variables on a seasonal scale in the seven grassland subtypes (Table 4). It was found that Ta could best explain the RE seasonal dynamics compared with other variables, in five of the seven grassland subtypes. Ta was also the secondary environmental factor explaining the RE seasonal dynamics in the remaining two grassland subtypes. To summarize, Ta was found to dominate the seasonal dynamics of RE in northern China's grasslands.  We further analyzed the relationship between RE and environmental variables on a seasonal scale in the seven grassland subtypes (Table 4). It was found that Ta could best explain the RE seasonal dynamics compared with other variables, in five of the seven grassland subtypes. Ta was also the secondary environmental factor explaining the RE seasonal dynamics in the remaining two grassland subtypes. To summarize, Ta was found to dominate the seasonal dynamics of RE in northern China's grasslands.

Environment Effects on RE Interannual Dynamics
With respect to each grassland type, we analyzed the interannual dynamic relationship between RE and environmental variables (Figure 8, Table 5). In the four alpine grassland subtypes, the R 2 between the annual RE and annual average Ta was 0.73-0.91 (Figure 8a). The R 2 between the annual RE and annual average EVI, excluding the alpine meadow steppe, was 0.44-0.63 (Figure 8c). There was no significant correlation between the annual RE and the annual total PRE or PAR (Figure 8b,d). In the three temperate grassland subtypes, the R 2 between the annual RE and the annual total PRE was 0.2-0.64 (Figure 8b). The R 2 between the annual RE and the annual average EVI was 0.28-0.63 (Figure 8c), but there was no significant correlation with the annual average Ta or the annual PAR (Figure 8a,d). In contrast to the positive promotion effect of PAR on RE at the seasonal scale (Figure 7d), the effect of annual total PAR on annual RE was not significant (Figure 8d). We found that when the mean annual PRE exceeded 400 mm, which was the case in all the alpine grasslands and in MS among the temperate grasslands (Table 5), annual average Ta always provided a positive stimulus for annual RE (Figure 8a). This suggests that in grasslands with adequate moisture availability, the interannual dynamics of RE were controlled by factors other than PRE. Compared with the significant positive correlation between annual average RE and annual total PRE in temperate grasslands, the correlation for these in alpine grasslands was negative or close to zero (Figure 8b). When  We found that when the mean annual PRE exceeded 400 mm, which was the case in all the alpine grasslands and in MS among the temperate grasslands (Table 5), annual average Ta always provided a positive stimulus for annual RE (Figure 8a). This suggests that in grasslands with adequate moisture availability, the interannual dynamics of RE were controlled by factors other than PRE. Compared with the significant positive correlation between annual average RE and annual total PRE in temperate grasslands, the correlation for these in alpine grasslands was negative or close to zero (Figure 8b). When temperature was restricted, the annual total PRE had almost no stimulating effect on the annual RE in the same grassland type. In summary, the interannual change in RE in alpine grassland is mainly attributed to the annual average Ta and annual average EVI, and the interannual change in RE in temperate grassland is attributed to the annual total PRE and annual average EVI.

Environment Effects on RE Spatial Variation
The red fitting line in Figure 8 shows that annual RE was correlated with the annual average EVI across the seven grassland subtypes, and the R 2 reached 0.91. This indicates that EVI can best explain the spatial variation in RE in northern China's grasslands, compared with other environmental factors. In alpine grasslands, including AS, SW, KO, and SH, the mean 15-year EVI showed an increasing trend, and RE also showed an increasing trend ( Figure 8c, Table 5). In temperate grasslands, EVI and RE also increased simultaneously from DS and TS to MS. These results highlight that RE and EVI changed synchronously in the same spatial direction in the same climate zone.
In the temperate grasslands, including DS, TS, and MS, the mean annual RE increased synchronously with an increase in the mean annual PRE (Figure 8b). Within AS, SW, KO, and SH in the alpine grasslands, the mean annual PRE and the mean annual RE also increased simultaneously. However, simultaneous variation in the PRE and the RE was only observed in the same climate zone. Across grassland types in the same climate zone, the mean annual PRE controlled the magnitude of the mean annual RE.
There was no significant correlation (R 2 = 0.004, p = 0.52) between annual RE and annual average Ta throughout northern China's grasslands (Figure 8a). In the alpine grasslands, including AS, SW, KO, and SH, the annual RE increased with an increase in the annual average Ta. However, in the temperate grasslands, there no synchronization phenomenon was observed. These results indicate that the dependence of annual RE on the annual average Ta tends to be stronger in regions with more moisture. The influence of Ta on RE in grassland ecosystems is regulated by moisture, and the regulation strength varies with different vegetation types.
Across the seven grassland subtypes, the SOCD, clay, and silt content had positive effects on the mean annual RE, whereas sand content had a negative effect (Figure 9), which is consistent with previous research [27]. Within AS, SW, KO, and SH in the alpine grasslands and DS, TS and MS in the temperate grasslands, RE increased simultaneously with the increase in SOCD and clay content. The R 2 between the RE and the SOCD, clay, silt, and sand content were 0.80, 0.68, 0.45, and 0.81, respectively, implying that soil factors also affect the spatial variation of RE. EVI and soil factors are critical environmental factors that determine the spatial variations in RE in northern China's grasslands. Their dominant roles illustrate the importance of grassland type and soil factors in regulating the spatial pattern of RE.

Model Comparison
Accuracy comparison between Geoman and the previous models in estimating RE in northern China's grasslands is displayed in Table 6. Geoman performed best, with the best explanation and minimum bias, when compared with seven previous models, including two semi-empirical models, one satellite-based model, three traditional machinelearning models, and one deep-learning model. Compared with semi-empirical models, deep-learning models can automatically extract the complex nonlinear dependence relationship between the target variable and the explaining variable, without excessive restriction from predetermined theoretical assumptions about complicated respiration processes. Based on similar methods, carbon flux estimations such as gross primary productivity (GPP) have been used as benchmarks and are widely employed in the evaluation and improvement of process models [62][63][64]. The combination of LSTM and the attention mechanism in the encoder and decoder stages further improves the predictive ability of the Geoman model (see Section 3.1.1). The LSTM in Geoman can adapt to time-varying behaviors and use dynamic data to update the parameters to improve the predictions. The local spatial attention module provides a means of distinguishing the importance of different input variables. The time attention module gains insight into past patterns of encoder outputs for different environment variables, and decoder outputs for the target variables. However, this insight can be difficult to model using traditional ML models. EVI and soil factors are critical environmental factors that determine the spatial variations in RE in northern China's grasslands. Their dominant roles illustrate the importance of grassland type and soil factors in regulating the spatial pattern of RE.

Model Comparison
Accuracy comparison between Geoman and the previous models in estimating RE in northern China's grasslands is displayed in Table 6. Geoman performed best, with the best explanation and minimum bias, when compared with seven previous models, including two semi-empirical models, one satellite-based model, three traditional machine-learning models, and one deep-learning model. Compared with semi-empirical models, deeplearning models can automatically extract the complex nonlinear dependence relationship between the target variable and the explaining variable, without excessive restriction from predetermined theoretical assumptions about complicated respiration processes. Based on similar methods, carbon flux estimations such as gross primary productivity (GPP) have been used as benchmarks and are widely employed in the evaluation and improvement of process models [62][63][64]. The combination of LSTM and the attention mechanism in the encoder and decoder stages further improves the predictive ability of the Geoman model (see Section 3.1.1). The LSTM in Geoman can adapt to time-varying behaviors and use dynamic data to update the parameters to improve the predictions. The local spatial attention module provides a means of distinguishing the importance of different input variables. The time attention module gains insight into past patterns of encoder outputs for different environment variables, and decoder outputs for the target variables. However, this insight can be difficult to model using traditional ML models. Despite the many challenges in adopting DL in ecology, DL can provide predictions with unparalleled precision without prior hypotheses. Geoman has been designed to make time-series predictions by understanding the impact of multiple explanatory variables on the target variables. However, the interpretation of the prediction results from the DL model structure is often insufficient, because the impact of the inputs on the outputs cannot be easily understood through neural networks [30]. In this study, an attempt was made to use the attention mechanism to explain the impact of environmental variables on RE, ecologically meaningful results were generated (see Section 3.1.2). The integration of DL models and ecological-process-based models represents a promising development in the domain of ecological study, and should substantially improve the interpretability of the DL models [65].

Climatic and Biotic Control over RE
Climatic conditions regulate the spatiotemporal patterns and variation in RE [1,66]. Temperature is the main factor affecting the RE rate in the absence of other environmental stresses [67]. This is because respiration requires the catalysis of enzymes, and temperature is a necessary condition for enzyme decomposition [68]. According to our results (see Section 3.3.2), warming promoted an increase in RE in alpine grasslands across the four alpine grassland subtypes as well as within the same subtype. It did not lead to an increase in RE in temperate grasslands, and even led to a decrease in RE in TS and DS (Figure 8a). Warming had no effect on RE, which may be related to the limited precipitation in temperate grasslands ( Figure S3b) that restricts plant production and microbial activity [69]. This also indicates that warming can increase the RE under favorable environmental conditions when moisture availability is not a restrictive factor. PRE provides the necessary water for vegetation growth, affecting canopy structure and plant production [70], and ultimately improves autotrophic respiration. PRE changes the soil microenvironment, affecting the dissolution and diffusion of organic solutes, gases, and enzymes [67,71]. This promotes microbial activity [72], and improves heterotrophic respiration. However, a saturation phenomenon affects RE under conditions of high precipitation (Figure 7b), which has been reported in previous studies [13,27]. This may be because precipitation in excess of the saturated soil moisture reduces the oxygen content in the soil [27,68]. It also inhibits the activity of plant roots and aerobic microorganisms, and reduces RE in grasslands. The different levels of saturated PRE for RE in alpine grasslands and temperate grasslands may reflect differences in the water requirements for different vegetation types.
The increase in the annual total PRE, and not the annual average Ta, had a positive effect on the annual RE within the same temperate grassland subtype (Figure 8a,b). This may be associated with the greater increase in precipitation than in temperature from 2001-2015 ( Figure S3a,b). The annual average Ta had a greater stimulating effect on the increase in annual RE in the same alpine grassland subtype (Figure 8a), which was associated with the increase in the annual average Ta ( Figure S3a). These results illustrate that PRE and Ta were the dominant climatic factors for interannual changes in RE in temperate grasslands and alpine grasslands, respectively. IM temperate grasslands are mainly distributed in arid and semi-arid zones, and limited moisture strongly restricts vegetation growth and microbial activity, which each affect RE changes [34]. Although plants that grow in alpine regions have adapted to the low temperature environment, their production and respiration processes are comparatively restricted. A minor temperature rise can contribute to an increase in RE by prompting the grassland to turn green, and increasing the biomass, the microbial activity, and litter decomposition [73][74][75].
The promoting effect of PAR on RE at the seasonal scale can be attributed to the effect of the seasonal PAR cycle on photosynthesis. Stronger photosynthesis tends to mean stronger respiration. However, at the interannual scale, the influence of annual PAR on annual RE changes was not significant in the alpine grasslands nor in the temperate grasslands (Figure 8d), which may be related to the long-term stability of solar radiation ( Figure S3d).
This study indicates that EVI is the most important environmental factor affecting the spatial pattern of RE in northern China's grasslands (Figure 8c), which is in accordance with previously reported results [12,16]. As a vegetation index commonly used for estimation of vegetation productivity, EVI is strongly associated with GPP in space and time [76,77]. GPP is the main supply for autotrophic respiration, and determines its intensity [24]. Decomposition of organic matter by microorganisms, including heterotrophic soil respiration, is also supported by photosynthesis. As shown in Figure S3c, the rate of increase for EVI in temperate grasslands from 2001 to 2015 was substantially greater than that in the alpine grasslands. This illustrates the faster increase in productivity of temperate grasslands, although the rate of increase for RE in temperate grasslands over 15 years was slower than that in the alpine grasslands (Figure 6c). This may be due to the greater uptrend of carbon sequestration in temperate grasslands [78].
The soil carbon substrate supplies various components of the RE, and determines the distribution of nutrients in aboveground and underground plant tissues, as well as the phenological patterns of vegetation growth [67]. Soil organic carbon (SOC) is one of the main forms of soil carbon that affects spatiotemporal patterns of respiration [1,79]. With a relatively large soil carbon pool and gradually melting permafrost, alpine grasslands on the TP have the potential to release significant amounts of CO 2 [80]. Our results also demonstrate that the RE and its rate of increase on the TP are larger than those in the IM (Figure 6c). SoilTex indirectly influences RE, primarily by controlling water and nutrient availability [81,82]. Clay and silt soils have higher soil water-holding capacities and longer water-holding times [83]. This provides more available water for vegetation growth, microbial decomposition, and respiration. Soils with higher clay and silt content can also store more SOC and nutrients for respiration processes [14,27]. The negative impacts of sand content on respiration could be attributed to larger soil porosity, resulting in lower soil water availability. This may constrain microbial activity and vegetation growth, leading to reduced respiration.
The effects of biotic and climatic factors on RE are strongly interlinked. Adequate hydrothermal conditions increase vegetation growth and autotrophic respiration [84]. These hydrothermal and vegetation conditions further regulate the SOC and affect the soil heterotrophic respiration [85]. Biotic and climatic factors synergistically stimulate RE, but at different spatial and temporal scales, and one of these variables may be particularly significant.

Limitations and Prospects
The carbon, water, and energy flux estimations for the regional scale based on the ML models represent a process of upscaling ecological knowledge of the carbon cycle from the site scale. Therefore, it is essential to evaluate the spatial representativeness of the flux observation sites [86,87]. Most of the current flux observation sites were established gradually over the previous few decades, with their locations often having been determined by expert experience and consideration of accessibility, including traffic, logistics, and topographic factors. Site density and uniformity substantially vary in different regions, owing to the lack of a systematic design. The 18 flux sites used in this study exhibited a nonuniform spatial distribution (Figure 1). In situ observations representing a few surrounding kilometers cannot fully represent the climate, environment, and flux in remote regions [88]. This increases the level of uncertainty in regional-scale carbon-flux estimations [89,90]. It has been reported that the current network of flux towers in northwestern Tibet lacks representativeness, and that adding flux towers could considerably improve the ability of the existing flux network to represent local environmental conditions [91]. It is expected that the ecological environment will receive more public attention, and the rapid expansion of big data in ecological systems engineering will establish an increasing number of flux sites, so their representativeness will become more credible.
The time range of carbon flux data employed in this study was concentrated in the period 2003-2014 (Table 1). This approximately coincides with the long-term regional estimation of RE. However, we obtained only one-year flux data for the AR, NMC, ZF, HLBE, SZWQ, and XLS sites. The small amount of observation data not only means these sites were under-represented, but may also make it impossible for the Geoman model to learn adequately the features of these flux data during model training. The estimated RE in this study is based on time-series environmental data with an eight-day resolution. These environmental variables cannot fully reflect various extreme climatic events that occurred during the eight-day interval [92], such as extreme heat waves or rainfall in summer and extreme cold in winter. The limited availability and considerable uncertainty of SOC and soil texture data restrict their effective application in estimating the RE at the regional scale [17,59,93]. The continuous sharing of ecological data, and the development of data augmentation and digital soil mapping technologies will provide reliable solutions to these problems.

Conclusions
Compared with several previous semi-empirical models and traditional ML models, Geoman performed better (R 2 = 0.87, RMSE = 0.39 g C m −2 yr −1 , MAE = 0.28 g C m −2 d −1 ) when estimating RE in northern China's grasslands. We found that soil factors have a significant effect on the estimation of RE in the low-value period in the non-growing season, thereby suggesting that RE models should not exclude the role of soil, especially in grassland ecosystems. Our results revealed the occurrence of substantial spatiotemporal changes in RE in northern China's grasslands. On the spatial scale, the RE of alpine grassland showed a decreasing trend from southeast to northwest in the TP, and that of temperate grassland showed a decreasing trend from northeast to southwest in the IM. On the temporal scale, the RE of grasslands in northern China showed a significant increasing trend (1.81 g C m −2 yr −1 ) from 2001 to 2015. The rapid increase in RE in alpine grasslands indicates greater response to climate change than in temperate grasslands. Biotic and climatic control over RE varied at the temporal and spatial scales. Ta was the main factor that influenced the seasonal dynamics of RE in grassland ecosystems. Ta and EVI controlled the interannual dynamics of RE in alpine grasslands, whereas PRE and EVI were the main controlling factors in temperate grasslands. EVI and soil factors were found to control the spatial changes in RE in grassland ecosystems. This study highlights the high accuracy of DL in simulating spatiotemporal changes in the carbon cycle, and enhances our understanding of the interactions between RE and climate change.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/rs14153563/s1, Figure S1: Eight-day time series plots for observed and predicted RE at all sites; Figure S2: Predicted RE vs. observed RE in the low value period in the non-growing season; Figure S3: Interannual dynamics of annual average Ta, annual total PRE, annual average EVI, and annual total PAR from 2001-2015 at the regional scale. Data Availability Statement: All data can be accessed through the provided website and references.