1 Introduction

Arid and semi-arid lands have low vegetation density and meager rainfall covering large areas (30–40%) on the Earth`s land surface (Bastin et al. 2017).These ecosystems are under the influence of extreme environmental conditions and the risk of global warming, while approximately 38% of the world’s population lives in them (Huang et al. 2017). These territories cause a remarkable amount of anthropogenic CO2 emissions and are expected to be affected more than climate change than humid regions (Berg et al. 2016; Bastin et al. 2017). The role of these ecosystems in the global carbon cycle is increasingly important and they are highly sensitive for the diversity of climatic factors (Ahlström et al. 2015; Biederman et al. 2018). In this context, many studies have been carried out using different methods in different ecosystems around the world to determine carbon exchange in dryland ecosystems and to understand the factors controlling these dynamics (Hastings and Oechel 2005; Jasoni et al. 2005; Wohlfahrt et al. 2008; Fu et al. 2009; Jia et al. 2016; Baldocchi et al. 2018; Kong et al. 2022; Yang et al. 2022). Some researchers have studied carbon dynamics in the shrub ecosystem and the controlling environmental factors (Biederman et al. 2018; Jia et al. 2018; Flores-Rentería et al. 2023). However, the temporal and spatial scales used in these studies to determine environmental factors controlling net ecosystem exchange (NEE) were different from Tsogt Ovoo in the Gobi Desert of Mongolia.

As mentioned by many researchers (Xie et al. 2015; Biederman et al. 2018), although drylands dominate the temporal variation of the global terrestrial carbon cycle, the role of these ecosystems in the global carbon cycle is still unclear. According to FAO (2023), shrubs cover approximately 12.5% of the world (1.9 billion hectares) and are the most widespread vegetation type worldwide. For these reasons, more studies are needed to determine the carbon dynamics of shrub-steppe ecosystems and the controlling factors of these dynamics.

The Gobi Desert, located in the south of Mongolia and the north of China, is a vast landscape with mountains and dunes. It is the 6th largest desert of the world with a size of 1.295 million km2. Moreover, it is also the largest desert of the Asian continent. It is not densely populated; the traditionally dominating nomadic lifestyle of pastoralists was mostly replaced by extensive farming in the last 70 years. Shrub ecosystems cover about 2.85 million ha area in Mongolia (FAO 2023). In Mongolia, some studies have been carried out to determine their carbon exchanges and the factors affecting them (Li et al. 2005; Shao et al. 2017; Nakano and Shinoda 2018; Wang et al. 2023). Shao et al. (2017) determined the total NEE of the shrub ecosystem in Mongolia during two growing seasons between June and September (without rain) to be 32 to 53 gC m−2. They found that the seasonal variations in net ecosystem carbon exchange were controlled mainly by soil water content. Furthermore, Nakano and Shinoda (2018) highlighted that the total NEE of grassland varied between -2 and -57.1 gC m−2 in Bayan-Unjuul during the growing seasons and interannual NEE related to precipitation, soil water content and vapor pressure deficit. These studies, however, were either on different ecosystems such as grassland, different locations, or environments in Mongolia on a limited time scale. Although deserts and drylands cover large areas of Mongolia, research on CO2 fluxes in Mongolia is limited and mostly carried out in the grassland ecosystems. While quite a lot of studies addressed similar ecosystems in Inner Mongolia (an autonomous region of China), like Butterbach-Bahl et al. ( 2011); Wang et al. ( 2017) no research was found in which the CO2 exchanges between the surface and the atmosphere were determined by the Eddy Covariance method over the shrub desert steppe ecosystem in multi-time scale in Tsogt-Ovoo district of Mongolia.

The aim of this study was to investigate the carbon dynamics using the Eddy Covariance method over the shrubland in the Gobi Desert, Mongolia (in the Tsogt-Ovoo district) during four growing seasons of 2014, 2015, 2018 and 2019 under the influence of the temporal distribution of different rain regimes in order to fill that gap in our knowledge to better understand this ecosystem. Another aim was to understand the factors controlling carbon dynamics of the shrub ecosystem using the machine learning method. This study is crucial in terms of being the first research conducted at this location in this ecosystem. We looked for answers to the following questions: (i) How do CO2 fluxes change over the shrub steppe ecosystem during the different rainy growing seasons (May–September) in the Gobi Desert, Tsogt-Ovoo district in Mongolia; (ii) Which factors control the carbon exchange (NEE) in desert ecosystems, to what extent and in which direction, as indicated by the SHAP values generated from the predictions?

2 Materials and methods

2.1 Site description

The measurements were conducted at the Tsogt-Ovoo district (44°23′04′′ N, 105° 16′ 59′′ E, 1232 m asl) in the Gobi Desert, Mongolia. The land was flat and homogeneous, and fetch length was 4 km in all directions (Kimura et al. 2016). Tsogt-Ovoo locates in the desert-steppe zone, which occupies 21.9% of the total land area of Mongolia (Jamsran et al. 2018). The long-term (1970–2019) annual mean temperature at the Tsogt-Ovoo was about 5.1 °C and annual average rainfall amount was about 104.4 mm. In general, most part of the rain (approx. 70%) falls between June and August. The growing season starts in May and lasts until mid-September (Buyantogtokh et al. 2022). In Tsogt-Ovoo, the surface is covered with snow from mid-October to the end of April. Other months are dry and cold compared to growing season. The soil of the observation station contains 11% clay, 12% silt and 77% sand. The vegetation community is dominated by shrub (Reaumuria and Salsola) community (Ishizuka et al. 2012).

2.2 Data collection

In this study, the flux measurements were carried out using the Eddy Covariance (EC) system established by the Arid Land Research Center (ALRC) of Tottori University to analyze the carbon exchange in the desert shrub ecosystem during the growing season from 1st May to 30th September in 2014, 2015, 2018, and 2019. CO2 flux measurements were made by using the Eddy Covariance system, which consists of a CO2/H2O open path infrared gas analyzer (LI-COR 7500A) and a 3D-sonic anemometer (Wind master Pro, Gill Instruments) and is installed at 3.06 m above the soil surface. In the EC system, data is stored in the data logger (CR1000, Campbell Sci.) in an interval of 10 Hz. At the same time, another station is installed to measure other meteorological variables. It consisted of air temperature and relative humidity (HMP-155D-10, Vaisala, at 1.96 m), precipitation (CTKF-1, OTA, at 0.65 m), air pressure (PTB210, Vaisala, at 1.5 m), wind speed and direction (YG-5103, Young, at 3 m), radiation (CNR4 radiometer, Kipp & Zonen, at 1.54 m) soil water content (ML2x ADR, Delta-T, at 1, 2.5, 5, 10, 15, 20, 30 and 50 cm depths) and soil temperature (C-PTG-10, Climatec, at 1, 2.5, 5, 10, 15, 20, 30 and 50 cm depths) sensors and the averaged data were stored in the datalogger (CR1000, Campbell Sci.) in an interval of 30-min. Sensors were maintained regularly due to dust transport in the desert area. In cases where the precipitation sensor was unable to measure due to the dust transport because of strong wind and excessive soil erosion, the missing precipitation data were completed from the station which was established by the Hydrology and Environment Institute of the Mongolian Meteorology Service (5 km away). Furthermore, due to data loss in the 2016 and 2017 growing periods, these years were not considered. Additionally, flux measurements could not be made between 1st and 10th May in 2014, and between 13th and 30th September in 2019 due to technical reasons. In some periods, SWC values at 30 and 50 cm depths could not be measured due to sensor failure. For this reason, SWC data up to 20 cm depth were used in the analyses. In the study, NDVI data of MOD13Q1 (Terra) and MYD13Q1 (Aqua) with a resolution of 250 m were used to monitor the development in shrubs during the measurement periods.

2.3 Data processing and gap filling

In this research, the recorded 10 Hz raw data are processed by using the software (EddyPro 7.0.9, LICOR Inc. USA, www.licor.com/eddypro). In these processes, double coordinate rotation, detrending (block average), detection and removal the spikes, rotation, spectral, cross wind, frequency response, density (Webb-Pearman-Leuning (WPL)) corrections were applied, and 30 min CO2 (NEE), LE and H fluxes data were computed (Webb et al. 1980; Lloyd and Taylor 1994; Falge et al. 2001; Reichstein et al. 2005). Angles attacks were corrected for Gill sonic anemometer. In addition to the above given steps, the half-hourly CO2 flux data were further removed from the dataset by considering precipitation, extreme fluxes, low quality data (QC = 2) and friction velocity. Subsequently, we used the R-based package REddyProc developed by the Max Planck Institute for Biogeochemistry (Reichstein et al. 2005; Wutzler et al. 2018) to fill the flux gaps in 30 min processed data and partitioning of net ecosystem exchange (NEE) to ecosystem respiration (Reco) and gross primary production (GPP). GPP was calculated as GPP = Reco – NEE. GPP = 0 at night and GPP ≥ 0 during daytime, whereas Reco ≥ 0 for both night and daytime. In this study, negative NEE fluxes measured at nighttime were considered as gap. Friction velocity (u*) threshold was applied for the gap filling, and it was determined as 0.11, 0.20, 0.10 and 0.11 m s−1 for 2014, 2015, 2018 and 2019, successively. Therefore, we assumed Reco = NEE for missing NEE at nighttime, (Reichstein et al. 2005).

Negative NEE flux indicates the flux towards the surface (i.e., uptake); in opposite case from the surface to the atmosphere (release). In energy fluxes, the positive values of the net radiation show that it is from the atmosphere to the surface, while the negative values ​​demonstrate that it is from the surface to the atmosphere. Additionally, positive values in LE and H fluxes indicate that there is a flux from the surface to the atmosphere, and in the opposite case from the atmosphere to the ground. 37% of the gap in the half-hourly NEE, 20% of LE and 7.4% of H data were rejected and filled during the data quality control and gap filling procedure.

2.4 Data analysis

In this research, the relationship between NEE and various environmental variables was predicted using the XGBOOST method. Three categories of input data were incorporated: meteorological factors, energy fluxes, and plant growth indicators, with NDVI being a key measure for plant growth. The SHAP algorithm was employed to identify the most influential input variables for predicting NEE. This approach allowed for a detailed exploration into how individual input variables contribute to the variations in NEE, providing a nuanced understanding of their impact. For the generation of SHAP plots, data from the XGBOOST models were utilized. This method enabled a more focused and in-depth analysis, leveraging the strengths of XGBOOST and SHAP in predictive modeling and interpretability, respectively. The use of these advanced techniques facilitated valuable insights into the complex interplay of environmental factors affecting NEE.

2.4.1 eXtreme Gradient Boosting (XGBOOST)

The XGBoost algorithm, developed initially by Friedman (2002) and later refined by Chen and Guestrin (2016), is an advanced machine learning technique using gradient boosting for enhanced accuracy in classification and regression problems. It involves sequentially training decision trees on a dataset, where each new tree is integrated using parallelization (Potdar et al. 2021). This method focuses on improving model accuracy by minimizing an objective function composed of a loss function, calculated through second-order Taylor expansion, and a regularization term to reduce model complexity and prevent overfitting (Koc et al. 2021). Distinguished from other boosting methods, XGBoost uniquely divides data into quantiles for finer interval determination, enhancing predictive accuracy but increasing computational load. It employs a deep-search-first approach for efficient search, starting from a random node and prioritizing deeper nodes, ensuring fast and reliable computation with diverse performance indicators (Tyralis et al. 2021). In this study, the hyperparameters of the XGBoost algorithm were optimized using grid search.

2.4.2 Shapley Additive Explanations (SHAP)

The SHAP algorithm effectively addresses the task of attributing the contributions of independent variables to the ultimate predictions of the dependent variable(s). The brilliance of the SHAP algorithm lies in its ability to offer a transparent and coherent ranking of predictors, achieved through the computation of Shapley values. Consequently, the comprehensive assessment of the prediction scheme hinges on the summation of Shapley values attributed to the input variables, utilizing principles from coalition game theory (Lundberg and Lee 2017). A noteworthy feature of the SHAP technique is its model-agnostic nature, allowing seamless integration into various machine learning algorithms, regardless of their inherent operational principles (Ekmekcioğlu and Koc 2022). The SHAP algorithm excels at revealing intricate relationships between predictors and targets, especially highlighting the impact of attribute value variations, particularly in cases involving extreme values. This facilitates a clear and straightforward representation of how these variations influence the occurrence of output samples (Ransom et al. 2022). Worth mentioning is the recognition within existing literature of the SHAP algorithm’s superiority over alternative techniques that aim to establish feature importance, such as permutation importance, Gini importance, and various correlation analyses like Pearson and Spearman correlations.

2.4.3 Performance evaluation criteria

In this study, statistical performance of the machine learning methods to predict NEE were calculated and compared by using root-mean-squared error (RMSE), coefficient of determination (R2) and Nash–Sutcliffe efficiency (NSE) index. For RMSE, the smaller values show higher accuracy. Also, NSE can take the values between -1 and 1, in which -1 indicates that the model is not able to predict outcomes and 1 denotes perfect predictions. The equations corresponding to the utilized criteria are provided Eq. 1, 2 and 3.

$$RMSE=\sqrt{\frac{1}{n}\sum\nolimits_{i=1}^{n}{({y}_{i}-{\widehat{y}}_{i})}^{2}}$$
(1)
$$NSE=1-\frac{{\sum }_{i=1}^{n}{\left({y}_{i}-{\widehat{y}}_{i}\right)}^{2}}{{\sum }_{i=1}^{n}{\left({y}_{i}-\overline{y }\right)}^{2}}$$
(2)
$${R}^{2}={\left(\frac{\sum_{i=1}^{n}({y}_{i}-\overline{y })({\widehat{y}}_{i}-\widehat{y})}{\sqrt{\sum_{i=1}^{n}{({y}_{i}-\overline{y })}^{2}}\sqrt{\sum_{i=1}^{n}{({\widehat{y}}_{i}-\overline{\widehat{y} })}^{2}}}\right)}^{2}$$
(3)

Variables given in the equations are; \({y}_{i}\)= actual values, \({\widehat{y}}_{i}\)= estimated values, \(\overline{y }\) = average of actual values, \(\overline{\widehat{y} }\) = average of estimated values, n = number of observations.

3 Results and discussion

3.1 Meteorological factors during growing periods

Temporal variation of meteorological factors in four growing seasons is presented in Figs. 1 and 2. A summary of descriptive statistics of daily meteorological data is given in Table S1 (in supplementary material). Total amount of rain in the growing seasons of 2014 and 2015 was 61 and 49 mm, respectively. The highest total rainfall was determined in the growing season of 2018 with 120 mm, whereas the lowest was in 2019 with 24 mm. Total precipitation in the seasons (except 2018) was less than its long-term average (approx. 89 mm). For this reason, the growing season of 2018 was wet, however others were dry. In the growing season of 2015 (May 10), 20.3 mm, and in 2019, approx. 100 mm (May 12) and 10 mm (May 13) snowfall occurred. When all measurement periods are compared, the maximum monthly total rainfall was determined in July (69.5 mm) and in August (38 mm) in 2018. The average soil water content (SWC) of all measurement periods from 1 to 10 cm (SWC(0–10)) and from 1 to 20 cm (SWC(0–20)) were calculated as 0.114 and 0.103 m3 m−3. Additionally, the average SWC(0–20) values of the growing periods of 2014, 2015, 2018 and 2019 were 0.082, 0.098, 0.105 and 0.116, respectively. In 2018, the highest month of both SWC(0–10) and SWC(0–20) was determined as August, with 0.193 m3 m−3. Among the measurement periods, the year with the highest average SWC(0–20) (0.151 m3 m−3) at the beginning of the development period (May) was in 2019. This is an indication of the excess of precipitation (rain and snow) before the growing season.

Fig. 1
figure 1

Time series of meteorological variables during growing seasons

Full size image
Fig. 2
figure 2

Time series of variables during growing seasons

Full size image

The average soil temperature at the surface (1 cm depth, Ts1) varied between 0.75 and 35.07 °C (Table S1). On the other hand, while the average soil temperature at 5 cm depth (Ts5) was determined as 22.6 °C, the lowest determined as 4.97 and the highest as 32.41 °C. In the 2019 season, the average Ts5 was higher than in other seasons. The average soil temperature for all measurements from the surface to 5 cm depth (1, 2.5 and 5 cm, Ts(0–5)) was 22.64 oC. The daily average air temperature (Ta) varied from 0.20 to 29.06 °C during the measurement periods. While the highest Ta was in 2019 as 20.1 °C, the lowest value was measured in 2015 with 18.3 °C. The average air temperature and relative humidity (RH) of all measurement periods were determined as 19.3 °C and 36.9%, individually. While the average RH was the lowest in 2019 (35.3%), the highest RH was found in 2018 with 39.7%.

In addition to these, there was no significant difference between the average global solar radiation (Rg) values in 2014 (268.8 W m−2) and 2015 (267.7 W m−2), however, it was 270.1 and 273.1 W m−2 in 2018 and 2019, respectively. Similarly, the average of daily-reflected short wavelength radiation (SWd) in the period of 2019 was approximately 10% higher than in other measurement periods. One of the reasons for this is the falling snow in May 2019. The average wind speed (u) of all periods was 4.21 m s−1 and this value changed between 0.32 and 11.98. Additionally, the dominant wind direction was south in all seasons. Moreover, the average vapor pressure deficit (VPD) was determined as 13.8 hPa in the whole measurement periods. The daily average albedo of the ecosystem was 0.241 of all measuring periods with the effect of shrub growth and environmental factors, whereas the highest albedo was 0.75 in May 2019. Additionally, in the period of 2014 and 2015, the mean NDVI values were very close to each other and were 0.10 and 0.11, respectively. The highest growing season average NDVI value was determined as 0.18 in 2018, whereas it was 0.13 in 2019. At the beginning of the development period, only in May 2019, the monthly mean NDVI value was approximately 30% higher than in other periods. While the average monthly NDVI value in May, June, and July in 2014, 2015 and 2018 was between 0.10 and 0.11, this value was between 0.13 and 0.14 in the relevant months in 2019. In August 2018, it reached the highest monthly average NDVI value of 0.32 compared to other years. Figure 3a,b,c,d present the status of surface area at the beginning of growing seasons (May) and Fig. 3d,e,f,g show the days when the NDVI was maximum. Furthermore, in our study, 15, 7, 17 and 5 rainy days (P ≥ 1 mm) occurred in 2014, 2015, 2018, 2019, respectively.

Fig. 3
figure 3

The state of the vegetation at the beginning of the growing seasons a) May 2014 b) May 2015 c) May 2018 d) May 2019 and on the days when the NDVI was maximum e) 2014 f) 2015 g) 2018 h) 2019

Full size image

The averages of the difference between the air- and the soil temperatures measured at 1 cm (Ta-Ts1), 5 cm (Ta-Ts5) and 10 cm (Ta-Ts10) depths were calculated as -4.29, -3.30 and -2.82 °C, respectively. The highest temperature difference was determined in 2019.

Daily average net radiation (Rn) varied -44.9 and 170.8 with an average value of 102.4 W m−2 among all measurement periods. The monthly average Rn reached the highest value with 108.7 W m−2 in 2018. In 2014, 2015 and 2019, the average Rn values were calculated as about 100.2, 100.9 and 99.9 W m−2, respectively. Additionally, among all periods, the daily average latent heat flux (LE) value was extremely low as expected, with an average of about 23.3 W m−2 (about 23% of Rn). The lowest average LE was approx. 15.3 W m−2 in 2014 and the highest was 34.1 W m−2 in 2018 period. The average LE determined in the 2015 season was approx. 24.4 W m−2. Although the precipitation in 2014 was higher than in 2019, the average LE in 2014 was less than in 2019. The reason may be the snowfall in May and the high soil water content caused by the precipitation that occurred before the growing period (May 1) of 2019. Furthermore, it has been determined that the daily average sensible heat flux (H) varied between approx. -18.6 and 122.5 W m−2 in all growing seasons. The highest monthly average H was found to be 65.1 W m−2 in 2014 and the lowest H was at 58.8 W m−2 in 2018. On the other hand, the daily mean H value was quite high (about 61.7 W m−2) due to dry conditions and constitutes about 60% of the daily mean Rn.

3.2 Variation of carbon dynamics

The variation of daily carbon dynamics (CO2 fluxes such as NEE, Reco, GPP) for different shrub growing seasons are presented in Fig. 4. Significant differences in CO2 fluxes were detected among four growing seasons. In all periods, the daily average net carbon uptake was determined as about 0.011 gC m−2 and NEE varied from -3.06 to 1.71. Moreover, mean values of NEE/Reco, NEE/GPP and Reco/GPP for all growing seasons were determined as -0.01, -0.01 and 0.99, respectively.

Fig. 4
figure 4

Time series of daily CO2 fluxes

Full size image

In the growing season of 2014, the daily average NEE, GPP and Reco were calculated as 0.55, 0.42 and 0.97 gC m−2, respectively, whereas these values were 0.41, 0.70 and 1.11 gC m−2, in 2015. The higher GPP value in 2015 compared to 2014 may be due to the 21 to 71% higher soil water content in July, August and September, higher evapotranspiration, and more intense precipitation in 2015. The reason for the increase in GPP may be due to the fact that a total of 49 mm of precipitation fell in just 7 days during the development period of 2015, the snowfall and melting of snow in May.

In 2018 season, much more precipitation than other seasons increased shrub growth (which can be seen in NDVI values) and caused to perform more photosynthesis. The frequency and high intensity and frequency of rain events in late July and August of 2018 (more than the long-term total average) and the resulting increased soil water content caused the shrubland ecosystem to rise in plant growth, especially in August, and eventually the carbon uptake of the shrub ecosystem increased significantly. As a result, the daily average net carbon sequestration value reached 0.72 gC m−2 in the growing season of 2018. This situation shows us that the amount and temporal distribution of precipitation can move carbon exchange from release to uptake in the desert ecosystem. This indicates that, as a result of future climate change in the relevant region, changes in the amount and frequency of precipitation may turn these desert surfaces into a source of carbon sinks. In the same season, Reco and GPP were determined as 0.71 and 1.43 gC m−2. In 2019 period, as in 2018, the shrubland ecosystem was determined as a carbon absorption source. In 2019, the daily average NEE, GPP and Reco were determined as -0.28, 0.66 and 0.37 g C m−2, respectively. However, the amount of carbon removed from the atmosphere in this period is less than in 2018. The reason for this may be the distribution and insufficiency of precipitation, less soil water content in July and August compared to 2018. In the 2018 period, the highest (approx. 45% of the total) carbon sink occurred in August with 49.62 gC m−2 alone. Moreover, in August 2018, when the total carbon storage was maximum, the monthly average value of LE was determined as the highest. However, 28% of the total carbon uptake occurred in July 2019.

Finally, the total amount of NEE for the measurement periods of 2014, 2015, 2018 and 2019 were calculated as 79.7, 63.3, -110.4 and -38.7 gC m−2, respectively. In our study, the shrub ecosystem was a carbon release in the growing seasons of 2014 and 2015, whereas it was a carbon sink in 2018 and 2019. The variation of carbon fluxes in the desert shrub ecosystem was very high from year to year. These results are compatible with the findings of Shao et al. (2017) in a different location of Mongolia. However, our overall carbon release findings differ from those of Shao et al. (2017). This might be due to precipitation amount, intensity, frequency, length of measurement period (1 month less in Shao et al. 2017), soil water content, plant growth status, and other environmental conditions. The occurrence of the highest carbon storage in 2018 can be explained by the excess of precipitation in this period. However, although the precipitation measured in 2019 was the lowest of all seasons, it was concluded that the shrub ecosystem’s role as a sink in this season is not related to the amount of precipitation solely in the measurement period. For this reason, in 2019, unlike the other two dry years (2014 and 2015), carbon absorption occurred in the shrub ecosystem. One of the reasons for this may be the increased soil water content because of autumn and winter precipitation pre-growing season. Another reason may be that the plant, which was well developed in the previous period, was protected by autumn and winter precipitation, until the 2019 growing season. Therefore, the average SWC(0–20) in May and June 2019 were higher than in other years. Furthermore, the highest monthly mean LE values were determined in May 2019. Additionally, the NDVI value was higher in 2019 than other dry seasons.

Our study indicates that carrying of the high level of the soil water content from the previous months (pre-growing season) over winter to the beginning of next growing period has a positive effect on shrub growth and carbon sequestration in spring season. Our results are consistent with the findings highlighted by Jia et al. (2016), Biederman et al. (2018) and Flores-Rentería et al. (2023) that precipitation (pre-growing season) in this context soil water content was crucial for the carbon balance of the desert shrubs. Additionally, our results showed that, similar to the results of Kimura (2011), NDVI, an indicator of plant development, started to increase suddenly when the daily rainfall exceeded 15 mm. It has been revealed that each amount of precipitation in the desert ecosystem does not affect soil water content and eventually photosynthetic activity.

In addition to all this, it should not be forgotten that even if the measurement system is regularly checked and maintained, extreme environmental conditions and dust erosion in the Gobi desert ecosystem might negatively affect measurements. In this study, only four growing seasons were analyzed. It would be useful to compare these results with studies to be conducted for a longer period. Furthermore, the differences in NEE values between growing seasons might also result from the grazing of sheep, camels and goats reared in the relevant research area as highlighted by Wang et al. (2017). Since grazing could not be monitored in this study, its effect on NEE could not be evaluated.

3.3 Factors controlling carbon dynamics

At the stage of determining the factors controlling carbon dynamics (net ecosystem exchange), daily data were analyzed separately for each development period of shrub and to cover four development periods. In this study, the relationships between CO2 fluxes (NEE) and meteorological factors, NDVI and energy fluxes were investigated by non-linear method such as XGBOOST, and were further analyzed with the SHAP algorithm. In the implementation of SHAP, Shapley values are determined to ascertain the impact of each predictor within a model on an individual case. This provides insights into how each attribute affects the model’s forecast for a specific data point. The process involves establishing a grand coalition that encompasses all the features in the dataset, in line with the principles of cooperative game theory, which the SHAP algorithm utilizes. The process then involves computing the incremental impact of each feature. This is done by evaluating all possible combinations of the features used and calculating the change in the model’s predictions when a particular feature is either included in or omitted from the set of instances. Ultimately, the Shapley values for a specific feature are derived by averaging these incremental impacts across all possible arrangements of the features. In the study, the grid search algorithm was utilized to optimize the hyperparameters of the XGBOOST algorithm. The hyperparameters of the XGBOOST algorithm were explored within the following ranges: maximum number of boosting iterations [1,100], L2 regularization term on weights [0,1], L1 regularization term on weights [0,1], and the learning rate [0,1]. During the modelling phase with XGBOOST, 70% of the data was used for training, while the remaining 30% was utilized as the test set.

In this study, different combinations of inputs were evaluated for the estimation of NEE values. SHAP values were utilized to determine the variables to be selected for NEE estimation. A total of 5 different estimation periods were considered, encompassing both the evaluation of measurements across all years together and the assessment of each individual year (2014, 2015, 2018 and 2019) separately. In order to avoid model complexity and prevent an increase in computation times, the estimation performance was primarily influenced by 6 input variables. Following the variable importance analysis conducted with the entire dataset, the influential parameters for NEE estimation were determined to be NDVI, Albedo, Ta, H, and Ta-Ts5, as illustrated in Table 1.

Table 1 Feature importance rank and mean SHAP values of all years
Full size table

The estimation performances of models utilizing all variables and the six selected variables through the feature selection method are presented in Table 2. In the model utilizing a total of 18 inputs, including all variables measured at the station and the remotely sensed NDVI values from the grid where the station is located, the XGBOOST algorithms yielded estimations with NSE values of 0.817. When considering only the NDVI variable as an input, the NSE value was calculated as 0.7. Each variable believed to enhance estimation performance was incorporated into the models, leading to an observed improvement in estimation performance up to six variables. The established performance achieved using the identified six variables closely approached that obtained by utilizing nearly 18 variables.

Table 2 NEE estimation results of all daily variable
Full size table

Scatter plots for some of the generated models are presented in Fig. 5. When Fig. 5 is examined, estimations made using all inputs exhibit a consistent distribution along the 1:1 line. However, models constructed with only a single input do not show a consistent distribution along the 1:1 line. Following feature selection, models constructed using the identified six variables also display a well-distributed estimation pattern.

Fig. 5
figure 5

Scatter plots of observed and estimated variables for all years

Full size image

Using XGBOOST models created with data from all periods, a SHAP summary plot is provided in Fig. 6. This graph is generated by arranging the most influential parameters for NEE estimation from top to bottom. The effects of independent variables’ low or high values on the dependent variable can be quickly discerned through this plot. Each point in the figure represents a single prediction, with the color scale indicating feature values. The color scale on the right side of the figure illustrates the high or low values of the respective features. Points on the right side of the figure (positive Shapley values) indicate prediction in favor of a high NEE value, while those on the left side (negative Shapley values) represent the opposite. The SHAP method also provides the significance of features. The higher a variable is positioned in the variable SHAP summary graph, the more important it is in terms of predictive logic. Upon inspecting Fig. 6, it is evident that higher NDVI values exert a decreasing influence on NEE values. A decrease in NEE values indicates an increase in carbon sequestration in the shrub ecosystem. Conversely, lower NDVI values result in a slight increase in NEE values, albeit not as significant as the decrease. This increased carbon release as a result of lower photosynthetic activity. When examining the albedo values, very low albedo values lead to a minor decrease NEE, whereas moderate albedo values tend to increase NEE. However, the influence of albedo on NEE is not as significant as the effect of NDVI. Changes in the other independent variables appear to have relatively less impact on NEE values.

Fig. 6
figure 6

Interpretability of the XGBOOST model using SHAP algorithm for all years

Full size image

The six input variables to be used were determined based on SHAP values for NEE prediction models created using 2014 data (Fig. S1, in supplementary material). These values are indicated in blue on the Fig. S1, while other variables with diminishing impact on estimation performance, highlighted in red, were not included in the model. After examining the mean SHAP values, the variables thought to be influential in NEE prediction are, in order, albedo, Ta-Ts1, Ta-Ts5, RH, VPD, and H. In the estimation models established using the data from the year 2014, when solely the albedo variable was defined as an input, the NSE value was calculated as 0.227. The outcome attained with the six identified variables closely matched the results achieved using around 18 variables. A noteworthy observation in the estimation models constructed using the 2014 data is the notable increase in estimation performance upon the inclusion of the H variable (Table 3).

Table 3 NEE estimation results of 2014 variable
Full size table

Subsequent to the process of feature selection, models constructed employing the identified six variables also reveal a well-distributed predictive pattern. Despite the favorable distribution of models developed using the selected six variables, the overall performance of the models in terms of estimation remains statistically below an acceptable average level (Fig. S2). This observation could offer researchers insights into the challenging nature of NEE value estimation during dry years.

From the analysis of Fig. S3, it is evident that lower albedo values have a positive impact on the increase of NEE values, while higher albedo values tend to decrease NEE values. In the case of Ta-Ts1 values, they exhibit an inverse relationship compared to albedo. Higher Ta-Ts5 values contribute to the enhancement of NEE values, while lower Ta-Ts5 values lead to a reduction in NEE. A decrease in RH values causes a slight increase in NEE, whereas an increase in RH values significantly decreases NEE. Lower VPD values marginally reduce NEE, but moderate VPD values are associated with an increase in NEE. Concerning the H value, both moderate and low values contribute to the elevation of NEE, whereas high values result in a decrease in NEE. However, as mentioned above, the correlation of these six variables with NEE is very low in 2014.

For the 2015 data, effective variables like Ts(0–5), Ta, NDVI, u, SWC(0–20), and VPD were identified for NEE estimation using SHAP values, as depicted in Fig. S4. Table 4 outlines the results for 2015 models using all variables. The XGBOOST algorithm with all 18 inputs reached an NSE of 0.500. When limited to Ts(0–5), the NSE was 0.189. Interestingly, models with only four variables surpassed the performance of those with 18 inputs. However, as observed in Fig. S5, only Model4 displayed a well-distributed scatter pattern, while others were less satisfactory. Fig. S6’s SHAP summary plot for 2015 data models illustrates varying impacts of Ts(0–5), Ta, NDVI, and other variables on NEE. From this figure, it’s observed that higher Ts(0–5) values positively influence the increase of NEE values, whereas lower Ts(0–5) values correspond to a reduction in NEE. Similar to the Ts(0–5) values, Ta values have a parallel effect. An increase in NDVI values leads to a decrease in NEE, while a decrease in NDVI enhances NEE. While rising u values slightly increase NEE, a drop in u values notably decreases NEE. Lower SWC(0–20) values cause an increase in NEE, but moderate to high SWC(0–20) values result in a decrease. Additionally, it appears that other independent variables don’t significantly affect NEE. In this dry growing period, the relationships between NEE and its controlling factors are higher than in the 2014 period.

Table 4 NEE estimation results of 2015 variable
Full size table

In 2018, as shown in Fig. S7, impactful variables for NEE forecasting included NDVI, H, SWC(0–10), Rg, albedo, and LE. Table 5 details the 2018 models’ performance, where the XGBOOST algorithm with all 18 inputs achieved an NSE of 0.922, and NDVI alone resulted in 0.615. Model5, using the XGBOOST algorithm, attained an NSE of 0.929.

Table 5 NEE estimation results of 2018 variable
Full size table

The 2018 period, characterized by increased moisture levels due to rainfall events, significantly contributed to making NEE estimations more feasible. During this period, notable variations in NDVI values were instrumental in enhancing the predictability of NEE values. Additionally, a positive impact of soil water content (SWC) values in the 0–10 cm soil depth range on the performance of the models has been noted (Fig. S8).

The SHAP summary plot depicting models generated utilizing the dataset from the year 2018 is displayed in Fig. S9. Upon close scrutiny of Fig. S9, it becomes apparent that diminished NDVI values contribute positively to the augmentation of NEE values, whereas elevated NDVI values lead to a marked reduction in NEE values. In the case of H values, it can be observed that lower values of H amplify NEE, whereas higher values of H diminish it. The remaining independent variables appear to lack a noteworthy influence on NEE.

After conducting the SHAP analysis, the variables identified as significantly impacting the prediction of NEE, in descending order of importance, are H, LE, Ta-Ts1, Rn, SWC(0–10), and RH, as shown in Fig. S10. The estimated performance of models using all variables for the year 2019 is detailed in Table 6. Focusing on the model that utilizes all 18 inputs, the XGBOOST algorithm yielded an estimation with an NSE value of 0.510. When only the H variable was used as an input, the corresponding NSE value for XGBOOST was 0.345. Similar to other models, the models using 6 independent variables from the year 2019 have generated results that closely align with those from models using 18 independent variables. However, a detailed examination of the NSE values indicates that not all models achieved satisfactory outcomes. This could be due to the complex nature of the year 2019, which experienced a transition from a wet to a dry season. The dynamic nature of this environmental shift might have affected the models’ ability to accurately capture the underlying patterns and dynamics of the system.

Table 6 NEE estimation results of 2019 variable
Full size table

It’s clear that Model0, which includes all variables, demonstrates the most effective scatter pattern for estimating NEE in 2019, as opposed to models with only six variables, which exhibit less satisfactory scattering, as seen in Fig. S11. The SHAP summary plot for models based on 2019 data, shown in Fig. S12, reveals that lower H values positively influence the increase in NEE values, while higher H values contribute to their decrease. Analyzing LE values, it’s noted that while lower LE values slightly reduce NEE, higher LE values substantially increase it. A decrease in Ta-Ts1 values is associated with a slight reduction in NEE, and an increase in these values significantly boosts NEE. Similar patterns are observed for SWC(0–10) and RH values.

4 Conclusion

In this study, we analyzed the dynamics of net ecosystem change determined by the Eddy Covariance method on a shrub-steppe ecosystem, as well as the controlling variables of this change over four growing seasons using a machine learning and SHAP methods. Until now, no such study was available for the Gobi Desert (Tsogt-Ovoo) of Mongolia. Therefore, this study shows the first detailed results of carbon dynamics for the shrub ecosystem in different growing seasons (wet and dry).

The results of our research provide additional information on variables regulating carbon dynamics (NEE) in the shrub ecosystem in Mongolia’s Gobi Desert. The CO2 flux measurements in shrubland ecosystem showed great temporal variability from season to season. Our results indicate that multiple environmental factors may influence the net ecosystem exchange in the desert ecosystem. However, it is not clear whether this ecosystem is a carbon sink or a source during dry seasons. The shrub ecosystem was a net carbon sink when there was sufficient water in the soil in dry seasons, too. This study holds valuable insights for the quest to identify the most influential factors affecting NEE and to discern whether these factors exert amplifying or diminishing effects on NEE. Notably, the research highlights the viability of employing non-parametric methods to evaluate the selection of factors and their effects specifically within desert ecosystems. In the context of desert ecosystems, the assessment of different factors and their effects is crucial due to the unique ecological dynamics of such environments. The use of parametric methods might be limited by assumptions that do not hold true in these settings, making non-parametric approaches like SHAP an invaluable tool. This study substantiates the efficacy of utilizing non-parametric methodologies for assessing factor selection and impacts, enhancing our comprehension of the intricate relationships within desert ecosystems.

The machine learning method showed that NDVI, an indicator of vegetation dynamics, is the most important variable affecting NEE, especially during the wet season. On the other hand, factors controlling NEE of shrub-steppe desert ecosystem are very variable in dry growing seasons. The common variables controlling NEE were determined as H, LE and SWC(0–10) in the growing periods when the soil water content was high. While soil water content was determined as one of the important variables controlling NEE in some growing seasons, a significant relationship between precipitation and NEE was found in neither wet nor dry seasons. Although it was found a significant relationship between the carbon uptake of the shrub ecosystem and environmental factors in the wet season, we could not find it for the dry seasons. In addition to these, our results showed that the well-developed shrub in the wet season can affect the carbon exchange in the next dry season. The results of this research are of great importance in understanding the carbon dynamics of the desert shrub steppe ecosystem and the factors that control it. However, the absolute size of the flux densities (being very small), the interannual variability (very high), the high sensitivity to environmental controls (like precipitation) or to farming practice (shrubs may only be present, when grazing is reduced), the natural landscape heterogeneity and probably also the vulnerability to extreme events (flood and storms) limit an easy extrapolation of our results. Instead, climate change and eventually population growth threaten these ecosystems along with their contribution to stabilize the carbon dioxide concentration of the atmosphere.