Untangling the impacts of socioeconomic and climatic changes on vegetation greenness and productivity in Kazakhstan

Studies examining the joint interactions and impacts of social-environmental system (SES) drivers on vegetation dynamics in Central Asia are scarce. We investigated seasonal trends and anomalies in drivers and their impacts on ecosystem structure and function (ESF). We explored the response of net primary production, evapotranspiration and normalized difference vegetation index (NDVI) to various SES drivers—climate, human influence, heat stress, water storage, and water content—and their latent relationships in Kazakhstan. We employed 13 predictor drivers from 2000 to 2016 to identify the interactions and impacts on ESF variables that reflect vegetation growth and productivity. We developed 12 models with different predictor–response variable combinations and separated them into two approaches. First, we considered the winter percent snow cover (SNOWc) and spring rainfall (P_MAM) as drivers and then as moderators in a structural equation model (SEM). SNOWc variability (SNOWcSD) as an SEM moderator exhibited superior model accuracy and explained the interactions between various predictor–response combinations. Winter SNOWcSD did not have a strong direct positive influence on summer vegetation growth and productivity; however, it was an important moderator between human influence and the ESF variables. Spring rainfall had a stronger impact on ESF variability than summer rainfall. We also found strong positive feedback between soil moisture (SM) and NDVI, as well as a strong positive influence of vegetation optical depth (VOD) and terrestrial water storage (TWS) on ESF. Livestock density (LSKD) exhibited a strong negative influence on ESF. Our results also showed a strong positive influence of socioeconomic drivers, including crop yield per hectare (CROPh), gross domestic product per capita (GDPca), and population density (POPD) on vegetation productivity. Finally, we found that vegetation dynamics were more sensitive to SM, VOD, LSKD and POPD than climatic drivers, suggesting that water content and human influence drivers were more critical in Kazakhstan.


Introduction
Asian drylands, which account for 30% of global drylands, are in developing countries where livelihoods rely on land and ecosystem services. These countries have experienced socio-ecological, environmental, and institutional shifts causing grassland degradation, livestock mortality and conflicts over water resources (Gutman et al 2020. Rapid population growth and ongoing climate change in these drylands are causing additional pressure on these fragile grassland ecosystems (Hao et al 2018). Consequently, Asian drylands are identified as landuse and climate change hotspots that are vulnerable to ecological and environmental degradation (de Beurs et al 2018). However, little is known about long-term vegetation changes, grassland degradation and drivers associated with human-environmental interactions (Abel et al 2021, Venkatesh et al 2022. Investigating vegetation dynamics and the underlying drivers is crucial for preventing further degradation and restoring degraded grassland ecosystems (Meyfroidt et al 2016).
Kazakhstan is an important land-locked dryland Asian country owing to its large size and economy (Schierhorn et al 2020). Previous studies have found that cropland conversions, land abandonments, and livestock density trends in Kazakhstan were some of the significant factors causing land degradation (de Beurs et al 2015, Rolinski et al 2021. Kazakhstan lost a lot of cropland cover due to socio-political and structural changes in the agricultural sector during the disintegration of the Soviet Union (Frühauf et al 2020). For e.g. crop yields dropped from 23.4 to 10.7 million tons and heads of livestock decreased from 48.6 to 14.5 million heads between 1990 and 2000 in Kazakhstan (Kraemer et al 2015). Livestock management in Kazakhstan is diverse with soviet-style feedlot-based ranching systems being the predominant type. However, private ranch owners cannot afford to move livestock to prime pastures in remote areas, which could alleviate grazing pressures. This is partially explained by the lack of infrastructure and funds to rebuild abandoned watering systems as well as rural outmigration (Dara et al 2020). On the other hand, livestock herders in eastern Kazakhstan have been engaging in transhumance-moving livestock to different ecological zones and from lower to higher elevations-to capture the benefits of seasonally accessible forage resources (Hankerson et al 2019). The availability of these forage resources depends on precipitation that occurs outside of the peak growing season; thus, the variation in seasonal precipitation greatly affects vegetation production (Tomaszewska and Henebry 2020). Furthermore, increasing temperatures have triggered snow cover and glacier mass loss which is expected to accelerate further (Luo et al 2019). Clearly, examining the variability in seasonal rainfall and snow cover in these grassland ecosystems is increasingly crucial, as they impact seasonal and long-term water storage, vegetation phenology and pasture productivity (Qi et al 2017, Petersky et al 2019. The vegetation-precipitation relationship (VPR) in drylands has been the primary focus of research for many decades (John et al 2013. Substantial changes in ecosystem structure and function of vegetation (ESF; list of abbreviations in appendix) has led researchers to evaluate the dynamics of greenness (normalized difference vegetation index-NDVI), water fluxes (evapotranspiration-ET) and ecosystem production (net primary productivity-NPP) in the context of trends in surface temperature, solar radiation, relative humidity, and vapor pressure deficit (Kong et al 2017. These studies suggest that precipitation variability is significantly responsible for alterations in ESF variables (NDVI, ET, and NPP). In addition, studies that investigated the impacts of spring drought on ecosystem carbon dynamics in semi-arid areas also found that this phenomenon had significant socio-ecological impacts (e.g. Zhang et al 2012. Snow cover, in addition to rainfall, is a key driver in Central Asia. Snowfall in Central Asia exceeds rainfall and is the prime contributor to the early summer soil moisture (SM) that drives seasonal biomass in pastoral lands (Apel et al 2018). Recent studies have explored the effects of snow cover change on vegetation over Central Asian regions (Venkatesh et al 2022). Most of these studies, conducted using remote sensing techniques and products, concluded that snow cover dynamics alter NPP (Wang et al 2018, Qiao andWang 2019). Hence, we sought to investigate whether winter percent snow cover and spring season precipitation affect peak season greenness in Kazakhstan using a causal model.
Previous research has focused on structural (e.g. measured by vegetation optical depth-VOD) and water retention (e.g. terrestrial water storage-TWS and SM) driver impacts on vegetation changes in different ecosystems (Deng et al 2020, Ugbaje andBishop 2020). Studies in dryland ecosystems have found a strong relationship between VOD and anthropogenic effects in the context of global climate change (Andela et al 2013. The SM, TWS-vegetation relationships are also crucial as intensive water withdrawals for irrigation often alter the water cycle components, particularly in waterlimited ecosystems. Previous studies have examined the individual effect of seasonal variations in moisture and water content drivers on vegetation greenness in Central Asia (Xie et al 2019). However, there has not been a comprehensive examination of the combined impact of moisture factors on greenness in Kazakhstan.
Scientific studies concerning the interactive effects of social-environmental system (SES) drivers on ESF changes are scarce, with only a few assessing either the effect of hydrological drivers (Xie et al 2016, Zheng et al 2019 or land cover/use change (John et al 2018, Dong et al 2020 on vegetation. Some studies have focused on coupled natural and human systems to investigate the complex ecosystem processes in semi-arid regions (Groisman et al 2009, Chen et al 2015a, 2015b. They suggest that socioeconomic changes or anthropogenic disturbances have produced more drastic impacts than climate change in recent years , Dong et al 2021. It is therefore necessary to examine the implications of these drivers on ESF dynamics in Kazakhstan, a country that has experienced significant land degradation and agricultural land abandonment (Prishchepov et al 2013, Hu et al 2020. A detailed analysis of the direct and indirect causal relationships between SES (climate, structural, water retention and socioeconomic) drivers and ESF response is needed and notably lacking over Central Asian areas (Tomaszewska et al 2020. In this study, we propose the following hypotheses to test whether there is a relationship between peak season greenness, spring precipitation and preceding winter snow cover for Kazakhstan and how this relationship is associated with land cover/use change (LCLUC): (H1) ESFs are positively and directly associated with spring precipitation and preceding winter snow cover; (H2) SM and structural water content (vegetation optical depth) have a strong influence on ESF dynamics; and (H3) Human influence has a more substantial and direct impact than climatic drivers on ESF changes.
We employed structural equation modeling (SEM) to help identify the combined interactive influences between SES drivers and the seasonal effect of these drivers on response variables (NDVI, NPP and ET). Our conceptual foundation of this integrated study of biophysical, climatic, and socioeconomic indicators driven by seasonal variation and change seeks to explain interconnectivity within coupled human and natural systems. We ask the following questions: (1) Does snow cover from the preceding winter and spring rainfall contribute jointly to peak season greenness and productivity? (2) Do socioeconomic indicators have a more substantial influence than precipitation in a water-limited ecosystem like Kazakhstan? (3) Does SM have a stronger impact than TWS on vegetation growth in this dryland region? and (4) Is vegetation affected by changes in anthropogenic water withdrawals (TWS) and water stress (VOD) in Kazakhstan? The study area description, as well as the dataset source and resolutions are listed in figure 1, table 1 and the appendix (sections 'Study area' and 'Data sources').

Methods
We implemented Mann-Kendall trend (MK) and Sen's slope estimator (SS) to identify monotonic upward or downward trends in the SES variables. We calculated standardized anomalies (Z) of all input datasets to maintain consistency. These standardized anomalies were used as input in SEM. More details regarding MK, SS and standardized anomalies are represented in the appendix (sections 'Mann-Kendall trends and standardized anomalies').
We carried out SEM analysis in two different phases by analyzing the interrelationships of ESF variables with various drivers at the provincial level in Kazakhstan (Fan et al 2016).
Phase 1: The first phase considers the two key drivers-precipitation and percent snow cover-as latent variables (description in appendix section 'Structural equation model (SEM)').
Phase 2: The second phase considers precipitation and percent snow cover as SEM moderators. We suggest that these two potential drivers alter the relationship between ESF and other water retention, structural and socioeconomic drivers. Though many previous studies identified that precipitation and percent snow cover would affect ESF, to our knowledge, its role as SEM moderator (appendix section 'Structural equation model (SEM)') has not been explored.
In the first phase, we developed eight SEM models, considering greenness (NDVI) as the response variable. Model-1 consists of four latent constructs . The driving variables were tested for influences using SEM and were removed if they were found to be non-significant or reduced model fit. The eight models with NDVI as a response are represented in table 3. Altogether, 24 SEM models were tested in the first phase of the analysis as we tested the interactions between drivers and three response (NDVI, NPP and ET) variables. In our second phase, we developed four SEM models by considering percent snow cover and seasonal precipitation as SEM moderators. The mean and variance of spring precipitation (P_MAM and P_MAM SD ) and snow cover (SNOWc and SNOWc SD ) were computed across years and tested as SEM moderators between SES drivers and ESF variables. In addition to the 24 models in the first phase, 12 models were tested in the second phase of the analysis to find the interactions between driving variables, SEM moderators and ESF variables. All these 36 models from two phases were developed and tested using SEM with a lavaan v0.6-10 package in R software (Rosseel 2012). Details regarding the SEM are represented in the appendix (section 'Structural equation model (SEM)') and the structural model implemented is illustrated in the path diagram (figure 2). More details regarding the fit statistics and SEM can be found in Fan et al (2016).

Standardized anomalies and spatial trends of response and forcing variables
Standardized anomalies help us understand longterm spatiotemporal dynamics across various provinces in Kazakhstan. The negative anomalies of P_JJA and SM showed strong associations across various provinces in Kazakhstan (supplementary figure  S1). VOD showed similar trends of NDVI, NPP and ET across Kazakhstan's eastern, western, and southern regions. Further, the negative P_MAM, VOD and NDVI anomalies over Kazakhstan's northern and southern provinces illustrated declining P_MAM and VOD effects on NDVI (supplementary figure S1). SNOWc and TWS showed similar trends with a significant increase in northern and north-eastern provinces indicating the contribution of winter snow melt water to summer ground water storage (supplementary figure S1). There was a general increase in LSK D (positive anomaly) across various provinces in  figure S2). The increasing and decreasing trends were depicted in green and brown, respectively, and the significant pixels (p < 0.05) were represented with dots (figure 4 and supplementary figure S2). Kazakhstan has experienced a significant decrease in NDVI, NPP, ET, SM, VOD, spring and summer precipitation and a significant increase of T air and LST in the western provinces of Kazakhstan. Significant greening (increasing trends of NDVI, NPP and ET) was observed in the eastern and northern provinces, along with decreasing temperatures (figure 4 and supplementary figure S2). Socioeconomic drivers, namely, POP D , LSK D and GDPca significantly increased in most Kazakhstan provinces. However, crop yields showed a non-significant trend in most provinces, with only a few provinces having significantly increased CROPh trends (table 2 and supplementary figure S3). A detailed explanation on standardized anomalies, spatial MK trends and Sen's slope estimates is provided in appendix (sections 'Standardized anomalies of response and forcing variables' and 'Spatial trends of response and forcing variables').

Joint interactions between ESF and driving variables 3.2.1. Phase 1: precipitation and snow cover as latent variables
The SEM model-1 indicated that P_MAM and LSK D (HINF-1; abbreviations in appendix) had a strong positive and strong negative influence with NDVI (standardized path coefficients (SPC) of 0.68 and −0.51)(supplementary figure S4(a)). On the other hand, joint interactive influence of POP D , CROPh and GDPca under HINF-2 showed a significant positive influence on NDVI (SPC of 0.44). Finally, winter SNOWc showed a weaker positive influence on summer NDVI in model-1 (SPC of 0.06).
Similarly, multiple models (model-2 to model-8; table 3) were tested by introducing P_JJA (model-2), VOD (model-3), SM (model-4), TWS (model-5), T air (model-6), and LST (model-7) as driving variables in SEM. These details were provided in appendix section 'Precipitation and snow cover as latent variables (model 2-model 7)' . The proximity variable (D CITY ), added in model-8 (table 3, supplementary figure S4(h)), did not result in either a significant positive or negative influence on NDVI, though it did improve model fit (table 3, RMSEA from 0.16 to 0.11). Model-8 exhibited better model fit statistics when compared to the other models in the first phase of the analysis (table 3). Therefore, the latent variables   (figure 5). The statistical fit of model-12 with SNOWc SD as SEM moderator exhibited the overall best fit of all the 12 models tested with NDVI as a response (table 3). The two-phase analysis described in section 2 was repeated to test the interactions between predictor variables, NPP and ET. We found that model-12 with SNOWc SD (figures 6 and 7) as SEM moderator exhibited better performance in terms of lower AIC (4485.56 & 4458.43), SRMR (0.07 & 0.07), and RMSEA (0.1 & 0.09) values for NPP and ET, respectively (supplementary tables 1 and 2). We also found that socioeconomic drivers had a more substantial and direct influence on productivity (SPC of −0.51 and 0.81 for HINF-1 and HINF-2) when compared to climatic and water retention parameters.
We observed key differences in driving factor impacts on response variables (NDVI, NPP and ET) (figures 5-7). We found a strong direct positive influence of HINF-2 on NDVI, followed by VOD and P_MAM. We also found a strong indirect negative influence of HINF-1 on NDVI through SNOWc SD (figure 5). For NPP, we found that HINF-2 had a strong direct positive impact, followed by HINF-1 and VOD. In contrast, the influence of P_MAM on NPP was weaker (figure 6). Similar to NDVI, HINF-1 exhibited a strong indirect negative influence on NPP through SNOWc SD . While HINF-2 had a strong direct positive impact on ET, followed by all three moisture drivers (VOD, TWS and P_MAM), we Figure 5. Structural equation modeling examining percent snow cover variability (SNOWcSD) as the moderator between mean normalized difference vegetation index (NDVI) and precipitation (PRECP), water content (WATRc), water storage (WATRs), human influence-1 (HINF1) and human influence-2 (HINF1) as constructs (i.e. latent variables). Model fit-chi-square (χ 2 ; degrees of freedom = 19) = 58.31, comparative fit index = 0.96; Tucker-Lewis index = 0.90; standardized root mean square residual = 0.07. All parameter estimates are standardized (full forms in appendix). Statistical significance: p < 0.05 * ; p < 0.01 * * ; p < 0.001 * * * , n.s.-non-significant. Abbreviations same as in figure 2. The square elements in the first row represent the drivers, and the circles in the second row represent the latent constructs. In the last row, NDVI represents the response variable and the SNOWcSD is the SEM moderator.
found that the influence of LSK D on ET was weaker (figure 7). In summary, VOD and HINF-2 showed a strong direct positive influence and HINF-1 exhibited strong indirect negative influence on ESF variables.

Hypotheses testing
Model-12 (table 3, supplementary tables 1 and 2), which exhibited the best model fit, was used to test our framed hypotheses. Hypothesis 1: P_MAM and SNOWc SD showed a positive impact on NDVI (SPC of 0.38 and 0.08), NPP (SPC of 0.21 and 0) and ET (SPC of 0.36 and 0.09). These results suggest that the proposed hypothesis (H1) can be accepted as both P_MAM and SNOWc had a significant positive impact (p < 0.05) on NDVI. Hypothesis 2: Soil moisture was removed from the final model as it resulted in poor model fit statistics. However, both SM and VOD showed a strong influence on NDVI (SPC of 0.46-supplementary  figure S6(f)). These results suggest that the H2 can be accepted as SM and VOD had a stronger impact on ESF variables. Hypothesis 3: We can accept our third hypothesis as socioeconomic variables under HINF-2 showed a more substantial impact than climatic drivers when interacting with ESF variables (figures 5-7). In addition, HINF-1 showed a stronger impact than climatic drivers when interacting with NPP and a similar impact to climatic drivers when interacting with NDVI (figures 5-7).
Our research addressed the four research questions that were framed in the introduction. We found that (1) spring rainfall had a stronger impact on ESF variables than winter snow cover (figures 5-7). However, we also found that the variance in winter snow cover has a higher influence on NDVI than spring precipitation (Model 12 was better than Model 9 and 10), (2) socioeconomic variables (HINF-1 and HINF-2) had a more substantial impact on ESF variables than precipitation in Kazakhstan (figures 5-7), (3) summer ground water storage showed a stronger influence than SM on NDVI (SPC of 0.

Discussion
Our research employed SES modeling to better understand vegetation dynamics in semi-arid regions. First, we used latent constructs to develop indicators that represent meaningful and complementary characteristics of SES drivers. Second, we tested the role of SEM moderators in explaining joint interactions of drivers on response variables. While the estimated latent construct values are site-specific, our methods suggest an approach to reduce predictor variables to represent vegetation ESF. Our findings could help guide researchers in prioritizing data collection if only a limited number of variables can be obtained. Such dimensionality reduction can be implemented in other water-limited ecosystems to improve sampling and interpretation capabilities.

Joint interaction of seasonal rainfall and winter percent snow cover on ESF
Very few studies have focused on coupled natural and human systems to explore the complex ecosystem processes that include the climatic and socioeconomic drivers in Central Asia. Some efforts have been made at smaller scales in Mongolia  (Yang et al 2016). Most of these studies used either regression models or trend analysis to examine the contribution of climatic and socioeconomic drivers to vegetation dynamics , Guo et al 2021. In contrast, we identified the joint interactive influences between SES drivers and the seasonal effect of these drivers on ESF using SEM in Kazakhstan. Our results showed that rainfall has a strong and significant positive influence on ESF across provinces in Kazakhstan. We also found that spring rainfall (P_MAM) contributed to increased greenness when compared to summer rainfall (P_JJA). The MK results also showed that NDVI and ET had similar trends of P_MAM. This was due to the significant time it takes for rainfall to convert into sub-surface moisture and contribute to vegetation growth (Zhang et al 2012, Zhou et al 2015. These inter-seasonal relationships are in agreement with previous studies that have attempted to analyze the influence of seasonal rainfall on vegetation in semi-arid ecosystems (Wen et al 2019, Shi et al 2021, Venkatesh et al 2022. Our findings also highlight that winter SNOWc has a significant but weak positive influence on peak season ESF in Kazakhstan. This could be due to ephemeral snowpacks that persist for less than 60 d (Petersky et al 2019). These ephemeral snowpacks have minimal predictable timing, and the snowmelt before the end of winter lowers the soil water availability for the peak growing summer season vegetation. Solomon et al (2009) suggested that with increasing spring and winter temperatures, precipitation would occur as rain rather than snow, particularly at the start and end of the winter, thereby enhancing the impact of rainfall on vegetation rather than snow. A similar phenomenon was observed in parts of Central Asia, as described by Chen et al (2016) and Tomaszewska et al (2020), indicating that the transition of snow to rain results in reduced snow and glacier ice accumulation during the winter. The MK spatial trends support these results as we found a reduction of SNOWc in most provinces of Kazakhstan.
We found that SNOWc has significant negative interaction with LSK D . This negative influence can be attributed to livestock deaths due to increased competition for natural forage and decreased quality and quantity of winter food availability during harsh winters in Kazakhstan (Hauck et al 2016, Mirzabaev et al 2016. The increase in frequency and severity of anomalous snow cover events due to significant warming and drying trends over Central Asia (e.g. Kazakhstan, Mongolia and Tibet) have caused high LSK D mortality rates , Nandintsetseg et al 2018.

Joint interaction of SM, VOD and TWS on ESF
Our SEM modeling indicated that SM has a strong and significant positive influence on ESF and strong positive covariance with P_MAM, VOD and TWS drivers. The MK trends also showed that SM had a similar spatial trend as VOD, NDVI, NPP and ET variables. The positive interaction between SM and vegetation suggests that the plant canopy had the function of water retention and storage. A similar positive feedback mechanism between SM trend and vegetation greening was earlier suggested by Li et al (2018) and Deng et al (2020). This positive feedback mechanism prevents extended soil drying, improves productivity and further prevents desertification and droughts.
We found that VOD has a significant positive relationship with ESF and strong positive covariance with SM, TWS and precipitation. These findings from SEM are similar to the results of MK spatial trends. This suggests that the alterations in moisture and storage variables affect vegetation water content, thereby affecting plant growth and productivity (Konkathi and Karthikeyan 2022). These results agree with the studies that employed VOD to understand vegetation dynamics in different ecosystems (Andela et al 2013, Tian et al 2018. The reduction of grassland canopy cover can explain the strong negative covariance between VOD and LSK D through overgrazing which strongly reduces plant water content . Future studies could perhaps focus on measuring the VOD variability at the early plant growth stage, which might help estimate the grassland canopy productivity at the end of the peak growing summer season. Our findings show that TWS has a significant positive influence on ESF. These findings are in agreement with previous studies that found a strong positive relationship between TWS and NDVI compared to SM and precipitation during summer (A et al 2015, Ugbaje and Bishop 2020). This strong positive relationship might be due to summer precipitation failing to establish a stronger connection with TWS that contributes to vegetation growth. In addition, vegetation growth could either diminish the available SM or convert to ET due to heat stress. This forces vegetation to depend on root zone water storage for survival, thereby exhibiting a strong relationship with TWS during the peak growing summer season. Hence, we suggest that examining TWS dynamics were equally important compared to the changes in precipitation and SM. We also found strong negative covariance between POP D , GDPca and TWS. This suggests an increase in water withdrawal and intensive irrigation agriculture, driven by increased population growth (Sun et al 2020). Future research could focus on measuring groundwater withdrawal rates in Kazakhstan as this is a major knowledge gap.

Joint interaction of socioeconomic drivers on ESF
We found a general increase in LSK D anomalies across provinces in Kazakhstan. The MK trends also showed a significant increase in LSK D in most provinces of Kazakhstan. This LSK D increase suggests higher grazing intensities that might cause a substantial decline in grassland plant diversity, NPP and resilience (Fetzel et al 2017, Liang et al 2018. Our SEM results also indicated that LSK D had a strong negative influence on ESF variables compared to other driving variables. Our results agree with the previous findings that grazing was a significant factor contributing to the decline in primary production compared to precipitation (Dangal et al 2016, Liang et al 2018. However, more research is needed to determine whether certain grazing patterns affect vegetation growth, particularly those used during the peak growing season. We also found a strong positive covariance between GDPca, POP D and LSK D . The MK trends showed that these three socioeconomic drivers are significantly increasing in most provinces of Kazakhstan. This can be attributed to the rapid increase in livestock density ( figure 3(a)) to meet the demand for meat by a growing population (Flammini et al 2013, Sans andCombris 2015). Recent studies found that many countries have transitioned from cereal-dominated diets to a rise in meat consumption following an increase in economic growth (Qi et al 2017). Such an increase in livestock density and subsequent intensive grazing has led to degradation of grassland ecosystems in Kazakhstan.
Our findings highlight the positive anomalies in CROPh and GDPca that could be explained by increased crop yields in Kazakhstan. The country experienced a drastic increase in croplands from 8% to 54% during 1953-1990, especially following the Soviet Union's virgin lands campaign in the 1960s (Swinnen et al 2017). In the 1990s, Kazakhstan experienced a drastic decline in cropland use after transitioning from state-ownership to a market economy. This resulted in a loss of assured markets, disintegration of value chain supplies, and deterioration of pricing linkages between inputs and outputs (Frühauf et al 2020). In addition, a reduction in livestock density and the demand for fodder drove a reduction in crop production, thus contributing to cropland abandonment during the post-Soviet period (Kraemer et al 2015). Increasing government support and investment in agricultural management since 2000 has led to 81% of cropland recultivation in the northern steppe region (Meyfroidt et al 2016). Furthermore, in southern Kazakhstan, cropland recultivation around major cities and a high percentage of irrigation agriculture near the Syr Darya river regions contributed to a sharp increase in cropland area (Klein et al 2012, Xi andSokolik 2016). These findings can be further validated from the MK trends that showed a significant increase of CROPh in northern (Qostanay and Aqmola) and South Kazakhstan provinces. Though CROPh and GDPca are rising, there are serious concerns regarding long-term sustainability in Kazakhstan, as it has limited water and energy resources for food production (Wright et al 2012, Chen et al 2015a. Grassland conversion to cropland affects water resources and causes environmental degradation due to reduced soil carbon, enhanced soil salinity and reduced groundwater levels (Lal 2011, Kulmatov 2014.

Conclusion
This study examined the spatiotemporal changes in ESF variables (NDVI, NPP and ET) during 2000-2016 using standardized anomalies, Mann-Kendall trends, Theil-Sen's slope estimates and structural equation modeling (SEM). We quantified direct and indirect interactions, covariances and joint influences of climatic, anthropogenic, heat stress, water storage and moisture drivers on ESF variables at the provincial scale in Kazakhstan. We found that Kazakhstan experienced a significant decrease in NDVI, NPP, ET, SM, VOD, spring and summer precipitation and a significant increase of T air and LST in the western provinces of Kazakhstan. We also found a significant greening (increasing trends of NDVI, NPP and ET) in the eastern and northern provinces, along with decreasing temperatures. The socioeconomic trends showed that POP D , LSK D and GDPca significantly increased in most Kazakhstan provinces. However, crop yields showed a non-significant trend in most provinces, with only a few provinces having significantly increased CROPh trends. Our SEM results suggest that water content (vegetation and soil) and joint interaction of human influence factors drive vegetation changes in Kazakhstan.
(a) We hypothesized and tested the effects of spring rainfall and winter snow cover on peak season ESF. We found that spring rainfall has a dominating impact over summer rainfall on ESF changes. Snow cover did not show a strong direct positive influence on ESF variables; instead, it played a moderating role, altering the influence of socioeconomic drivers on ESF variables. (b) VOD and SM, the critical drivers for vegetation growth, exhibited a strong positive influence on ESF and strong positive covariance with P_MAM and TWS variables. We identified positive feedback between vegetation and SM, indicating that the plant canopy had the function of water retention and storage. Furthermore, we found a strong positive impact of TWS on vegetation changes in Kazakhstan, a dryland ecosystem with diminishing surface water resources. (c) We also found an increase in livestock density and its enhanced negative impact on land degradation and productivity as measured by ESF dynamics. We observed a strong positive relationship between livestock density, population density and economic growth. This is attributable to the increased demand for livestock numbers and a growing population that has shifted its dietary habits from cereals to meat in Kazakhstan. Furthermore, positive anomalies in CROPh were explained by cropland recultivation in northern and southern Kazakhstan regions. The increase in livestock density, grassland conversions, a transition from snow to rain and strong dependence of vegetation on subsurface moisture affect water availability and causes grassland degradation in Kazakhstan. Thus, there is a need to develop alternative approaches to limit overgrazing and grassland conversion and maximize forage production in Kazakhstan.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Ethical approval
Our study does not involve human subjects and/or animals. The manuscript in part or in full has not been submitted or published anywhere.

Funding statement
This study was supported by the LCLUC Program of NASA (80NSSC20K0410).

Variable
Full

Study area
Kazakhstan is the largest landlocked country (2.72 million km 2 ) in dryland Asia (40 −55 • N, 50 −85 • E). The country comprises 14 regions (i.e. provinces) with grasslands as the dominant land cover type, occupying over 1.45 million km 2 (figure 1). Kazakhstan has a typical continental arid climate with extreme winters and dry, hot summers (Yuan et al 2022). The annual rainfall ranges from 100 mm in arid grasslands to 200 mm in semi-arid grasslands, reaching 900 mm in montane grasslands in alpine regions (Yan et al 2020a). The average temperature over the country fluctuates widely, ranging from −20 • C in January in north and central regions, to 18 • C in the North and 29 • C in the South by July. It was estimated that ∼60% of Kazakhstan (i.e. ∼180 million hectares) is currently experiencing desertification (de Beurs and Henebry 2004). Prior to Soviet Union's disintegration in 1991, a few large stateowned organizations dominated agricultural production. After 1991, in spite of large cooperatives by private landowners, there were serious constraints on agricultural productivity in Kazakhstan. These were mainly due to a cessation of trading agreements; decline in regional demand for cereals and other food grains; improper planning for trade; transportation and storage; minimal government support for innovation and development in agricultural research; and higher costs for pesticides (Meng et al 2000, Shmelev et al 2021. However, by 1998 almost 98% of the croplands were managed by private landowners (Suleimenov andOram 2000, Schierhorn et al 2020). These changes have led to reductions in areas under cultivation, livestock units and production of food grains (Baydildina et al 2000). Kazakhstan also experiences frequent droughts resulting in lowered productivity and increased interannual variability of agricultural yields (Kim et al 2021).

Data sources
Ecosystem, climatic and socioeconomic data records: NDVI, NPP and ET datasets were obtained from the moderate resolution imaging spectroradiometer (MODIS) dataset (table 1). Precipitation and air temperature were obtained from ECMWF reanalysis (ERA-5) datasets. Cumulative spring and summer precipitation (P_MAM and P_JJA) and mean summer air temperature (T air ) were computed for each year over the study period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016). Normalized difference snow index (NDSI) data were used to derive a maximum percent snow cover (SNOWc) composite over winter for each province (Hall and Riggs 2016). A merged (active and passive sensor) SM dataset from the European space agency-climate change initiative (ESA CCI) project was used in the study (table 1). The vegetation optical depth (VOD) dataset developed from X-band Radar was used to examine the structural behavior of vegetation (Moesinger et al 2020). We obtained TWS data from GRACE and GRACE-FO datasets produced by NASA JPL (table 1). We also used the NASA global land data assimilation system (GLDAS) derived land surface temperature (LST) (table 1). To maintain consistency, z-scores were calculated for each dataset and processed to obtain each province's spatial mean. Socioeconomic variables, including livestock density (LSK D ), population density (POP D ), gross domestic product per capita (GDPca) and crop yield (CROPh), were obtained for each province from the Kazakhstan Bureau of National Statistics and normalized by the areal extent of the province. Our first hypothesis was to test whether the preceding winter snow cover and spring season precipitation impacts peak season (summer) greenness. Hence, only spring (March-April-May; MAM) and summer (June, July and August; JJA) precipitation and winter snow cover (December-January-February; DJF) were considered, as Kazakhstan is a semi-arid country with moisture being a major limiting factor (supplementary figure S8).
Furthermore, Kazakhstan like most countries, has annual records, while monthly or seasonal datasets are rare. Thus, we used annual socioeconomic data that was available. Grasses/crops in semi-arid regions use surface water available in the initial phase (April-May during sprouting) from precipitation/snow and later depend on soil moisture (SM during June-July-August from supplementary figure S8) and groundwater storage during the peak growing season. Hence, we have used summer SM and terrestrial water storage (TWS) datasets in the analysis.

Mann-Kendall trends and standardized anomalies
Mann-Kendall trend (MK) and Sen's slope estimator (SS) are non-parametric statistical tests used to identify monotonic upward or downward trends in the SES variables. The MK and SS are considered robust as they (a) avoid presumptions of data distribution and data skews, (b) correct for serial autocorrelation, (c) are less sensitive to outliers and (d) can handle abrupt changes due to non-uniform time series. The MK Tau (τ ) ranges from −1 to 1, with positive and negative values indicating an increasing and decreasing trend. The Sen's slope magnitude (change per unit time), either positive or negative, signifies the trend's strength. The significance of trends is determined based on the p-value. The null hypothesis indicating no significant trend will be rejected if the p-value is less than 0.05. The MK trend and Sen's slope was calculated using Kendall, trend and SpatialEco packages in R software, version 4.0.3. More details regarding the equations and usage can be found in Mann (1945), Sen (1968) and John et al (2016).
Anomalies are the difference between a variable's current value and the long-term mean. Standardized anomalies are obtained when these anomalies are divided by the long-term standard deviation of that variable (John et al 2013). We calculated standardized anomalies (Z) of all input datasets to maintain consistency as: where z is the standardized anomaly of a response/predictor variable (e.g. NDVI, SM, etc) for a specific year (2020) relative to the long-term mean (x) and standard deviation (σ) from 2000 to 2016. These standardized anomalies were used as input in SEM.

Structural equation model (SEM)
The multivariate causal relationships in our SEMs were modeled by employing two or more structural equations (Fan et al 2016). The hypothetical dependencies between variables based on path analysis were tested using the structural model, whereas the latent variables were measured using the measurement model in SEM. Latent variables in SEM are unquantified variables whose impact can be estimated using one or more predictor variables. Latent variables are important for capturing complex system properties that are difficult to estimate physically. These latent and measured variables in SEM are framed based on theoretical knowledge and developed to test competing hypotheses regarding processes accountable for dynamics in the data (Chen et al 2015a). A significant advantage of SEM is that it explains the covariances among different variables instead of correlations and can handle non-linearities, hierarchical paths, categorical variables and complicating factors Keeley 2006, Grace et al 2012). Furthermore, we have tested the role of SEM moderator in this study. An SEM moderator is a quantitative variable that impacts the strength of the association between a predictor and response variable (Giannico et al 2021).
We developed eight SEM models in the first phase of the analysis. Model-1 structure was described in section 2 and represented in supplementary figure  S4(a). Model-2 introduced P_JJA under PRECP2 LC to test its impact on summer greenness (supplementary figure S4(b)). Model-3 introduced VOD as a separate LC (termed water content or WATRc) to test its association with greenness (supplementary figure S4(c)). Model-4 presents SM parallel to VOD in WATRc LC as VOD and SM are proxies of plant and soil water content, respectively (supplementary figure S4(d)). Model-5 includes a new LC (WATRs) with TWS, a proxy of water storage (supplementary figure S4(e)). Model-6 tests the influence of T air (heat stress-HEATs LC), that might exert a negative impact on greenness (supplementary figure  S4(f)). Model-7 introduces LST parallel to T air in HEATs LC to test the combined effect of temperature on greenness (supplementary figure S4(g)). Model-8 introduces a proximity variable, i.e. distance to cities (D CITY ) in the HINF-2 LC (supplementary figure  S4(h)). We suggest that D CITY and greenness might have a positive association, as increasing proximity to cities decreases ecological disturbance. In total, 24 models were tested in the first phase of the analysis as we are testing the associations between driving variables and three response variables (NDVI, NPP and ET). The SEM was tested for 17 years across 14 provinces with a maximum of 11 driving variables in a model. Therefore, we maintained a ratio of 238 (17 × 14):11, i.e. ∼21:1 for a model with highest number of variables (table 3, supplementary tables 1 and 2). SEM goodness of fit was assessed with a chi-squared value ( χ 2 ) , comparative fit index (CFI), Tucker-Lewis index (TLI), standardized root mean square residual (SRMR), root mean square error of approximation (RMSEA) and akaike information criterion (AIC). A good model fit is defined when the CFI and TLI are >0.9, SRMR is <0.08, RMSEA is <0.06 and lower AIC among different models.

Standardized anomalies of response and forcing variables
NDVI showed significant positive anomalies across eastern and western Kazakhstan regions but negative anomalies in the southern (South Kazakhstan and Zhambyl) and over a few northern (Qostanay and Pavlodar) provinces ( figure 1 and supplementary  figure S1). NPP and ET showed similar trends with significant negative anomalies in the western (Atyrau and Mangghystau) and south-western (Qyzylorda) regions and positive anomalies across central and eastern parts of Kazakhstan (figure 1 and supplementary figure S1). P_MAM showed a significant increase in western and south-eastern provinces of Kazakhstan, whereas there was a significant decrease in northern regions (North Kazakhstan, Aqmola and Pavlodar). P_JJA manifested an opposite trend to P_MAM across the entire country. SNOWc revealed negative anomalies over the northern and northwestern (Aqtobe and Qostanay) part of Kazakhstan, with positive anomalies over central and northeastern provinces ( figure 1 and supplementary figure  S1). T air and LST exhibited a negative anomaly (decrease in temperature) over Kazakhstan's central portion, with a slight increase in temperatures over the other regions. SM and VOD exhibited similar trends, with a significant decrease in northern regions and a significant increase in other country regions (supplementary figure S1). TWS showed positive anomalies over northern and north-eastern parts of Kazakhstan, with significant negative anomalies over western and southwestern regions (Qyzylorda- figure 1 and supplementary figure S1).

Spatial trends of response and forcing variables
The MK test results for NDVI, NPP and ET revealed that the western provinces are experiencing decreasing trends in contrast to eastern provinces that have increasing trends. West Kazakhstan and Aqtobe provinces experienced a significant decrease (p < 0.05) in NDVI, NPP and ET, whereas East Kazakhstan and Almaty provinces experienced increased NPP (figure 4). NPP and ET had a magnitude ranging from −10 to 10 (g C m −2 /season and mm/season) across the country, whereas the magnitude of NDVI ranged from −0.03 to 0.03 (supplementary figure S2). P_MAM, P_JJA and VOD showed similar trends to NDVI, with decrease in West Kazakhstan, Mangghystau, Aqtobe, Almaty and Zhambyl provinces (figure 4). Spring precipitation had significant decreasing trends (magnitude <−5 mm/season) over the Pavlodar region and some portions of Aqtobe province. Summer precipitation had significant decreasing trends (magnitude <−5 mm/season) in West Kazakhstan, Aqtobe and some parts of South Kazakhstan, Zhambyl and Almaty provinces (figure 4 and supplementary figure S2). Though VOD had significant decreasing trends in West Kazakhstan, Aqmola, Pavlodar and some parts of Aqtobe, Qostanay, North and South Kazakhstan provinces, the magnitude of change was minimal across the country.
SNOWc exhibited significant increasing trends in West Kazakhstan and Atyrau regions. In contrast, there was decreasing trend in the central parts of Kazakhstan along with the mountainous regions of Almaty and East Kazakhstan provinces (figure 4). SM followed a similar trend to ET, with decrease in western provinces and in Zhambyl and Almaty provinces. Similar to VOD, SNOWc and SM had a minimal magnitude of change across the country (±2%/season and ±0.002, respectively; supplementary figure S2). T air and LST showed similar trends in all provinces of Kazakhstan. A significant increase in temperature was found in western provinces along with a few provinces in southern parts of Kazakhstan. LST showed decreasing trends in North Kazakhstan and over the mountainous regions of Almaty and East Kazakhstan provinces (figure 4). T air and LST had a magnitude of ±0.1 degree/season across the country.
The MK trend for socioeconomic drivers showed that POP D had increased in most provinces except Qostanay, East and North Kazakhstan (table 2  and supplementary figure S3). The increasing and decreasing trends are significant (p < 0.05) for all provinces except Pavlodar. South Kazakhstan was experiencing the highest rate of increase (0.45 persons km −2 yr −1 ), followed by Almaty and Aqmola (0.3 and 0.18 persons km −2 yr −1 , respectively) provinces (table 2). LSK D showed a significant increasing trend in all provinces of Kazakhstan, with a higher magnitude of increase in South Kazakhstan (1.48 heads/area of aimag/year), followed by Zhambyl and Almaty (0.51 and 0.45 heads/area of aimag/year) provinces. Mangghystau and Qyzylorda have a lower increase in magnitude and exhibited a non-significant trend (p > 0.05). GDPca showed a significant increasing trend in all provinces of Kazakhstan with a higher magnitude in Atyrau (1648.5 GDP per capita/year) followed by Mangghystau and Aqmola (791.42 and 685.95 GDP per capita/year) provinces (table 2 and  supplementary figure S3). GDPca had a minimal increasing trend in South Kazakhstan and Zhambyl (188.08 and 193.52 GDP per capita/year) provinces, where LSK D was higher. CROPh showed significant trends in eight out of 14 provinces with a higher magnitude in Aqmola (0.009 tons ha −1 yr −1 ), followed by South Kazakhstan and Qostanay (0.006 and 0.005 tons ha −1 yr −1 ) provinces. Aqtobe and West Kazakhstan showed decreasing CROPh trends, whereas Atyrau and Almaty had higher significant increasing trends.

Precipitation and snow cover as latent variables (model 2-model 7)
In model-2 (supplementary figure S4(b)), P_MAM (SPC of 0.64) exhibited a strong positive influence on NDVI compared to summer precipitation-P_JJA (SPC of 0.29). P_JJA was removed from further SEM computations as it negatively impacted fit statistics (table 3, AIC from 4143 to 4667). Model-3 introduced VOD (WATRc LC) and showed that VOD improved the model fit and had a strong positive influence (SPC of 0.58) compared to other driving variables (supplementary figure S4(c), model-3, table 3). SM, introduced in model-4, showed a significant positive influence on NDVI in combination with VOD (supplementary figure S4(d)) that explained approximately 75% of the variability in NDVI (SPC of 0.75). The joint interaction of SM and VOD reduced all other predictor variable's impacts on NDVI, including the strong positive influence of P_MAM (SPC from 0.69 to 0.22). However, SM was removed from successive SEM models as it resulted in poor model fit statistics (table 3, AIC from 4143 to 4589). TWS, introduced in model-5 (table 3, supplementary figure S4(e)), showed a significant positive influence on NDVI (SPC of 0.29). T air and LST (Heat stress LC) added in models-6 and 7 (table 3, supplementary figures S4(f) and (g)) revealed a significant negative interaction with NDVI. The joint interaction of T air and LST degraded the model fit (table 3, AIC from 4143 to 5326), so it was removed from further SEM analysis.