Beyond “greening” and “browning”: Trends in grassland ground cover fractions across Eurasia that account for spatial and temporal autocorrelation

Grassland ecosystems cover up to 40% of the global land area and provide many ecosystem services directly supporting the livelihoods of over 1 billion people. Monitoring long‐term changes in grasslands is crucial for food security, biodiversity conservation, achieving Land Degradation Neutrality goals, and modeling the global carbon budget. Although long‐term grassland monitoring using remote sensing is extensive, it is typically based on a single vegetation index and does not account for temporal and spatial autocorrelation, which means that some trends are falsely identified while others are missed. Our goal was to analyze trends in grasslands in Eurasia, the largest continuous grassland ecosystems on Earth. To do so, we calculated Cumulative Endmember Fractions (annual sums of monthly ground cover fractions) derived from MODIS 2002–2020 time series, and applied a new statistical approach PARTS that explicitly accounts for temporal and spatial autocorrelation in trends. We examined trends in green vegetation, non‐photosynthetic vegetation, and soil ground cover fractions considering their independent change trajectories and relations among fractions over time. We derived temporally uncorrelated pixel‐based trend maps and statistically tested whether observed trends could be explained by elevation, land cover, SPEI3, climate, country, and their combinations, all while accounting for spatial autocorrelation. We found no statistical evidence for a decrease in vegetation cover in grasslands in Eurasia. Instead, there was a significant map‐level increase in non‐photosynthetic vegetation across the region and local increases in green vegetation with a concomitant decrease in soil fraction. Independent environmental variables affected trends significantly, but effects varied by region. Overall, our analyses show in a statistically robust manner that Eurasian grasslands have changed considerably over the past two decades. Our approach enhances remote sensing‐based monitoring of trends in grasslands so that underlying processes can be discerned.

Fractions (annual sums of monthly ground cover fractions) derived from MODIS 2002-2020 time series, and applied a new statistical approach PARTS that explicitly accounts for temporal and spatial autocorrelation in trends. We examined trends in green vegetation, non-photosynthetic vegetation, and soil ground cover fractions considering their independent change trajectories and relations among fractions over time. We derived temporally uncorrelated pixel-based trend maps and statistically tested whether observed trends could be explained by elevation, land cover, SPEI3, climate, country, and their combinations, all while accounting for spatial autocorrelation. We found no statistical evidence for a decrease in vegetation cover in grasslands in Eurasia. Instead, there was a significant map-level increase in non-photosynthetic vegetation across the region and local increases in green vegetation with a concomitant decrease in soil fraction. Independent environmental variables affected trends significantly, but effects varied by region. Overall, our analyses show in a statistically robust manner that Eurasian grasslands have changed considerably over the past two decades. Our approach enhances remote sensing-based monitoring of trends in grasslands so that underlying processes can be discerned.
Grasslands also provide many ecosystem services (Bengtsson et al., 2019), support wellbeing of more than 1 billion people, and regulate climate through carbon sequestration and albedo-driven air circulation (Bengtsson et al., 2019;Lioubimtseva & Henebry, 2009;O'Mara, 2012). Anthropogenic and environmental factors influence grasslands at different spatial and temporal scales and may lead to degradation and vegetation recovery (Bardgett et al., 2021;Cherlet et al., 2018;FAO, 2005). While identification of abrupt changes with moderate to high magnitudes and their causes is typically straightforward, detection of subtle changes taking place over longer time periods is more challenging, as is identification of significant drivers causing these trends. Appropriate and statistically valid assessment of changes in grasslands and identification of their drivers is pivotal to support biodiversity, food security, and climate change modeling and mitigation efforts (IPBES, 2018;IPCC, 2019;O'Mara, 2012), including meeting global environmental and development goals (Cowie et al., 2018), and UN ecosystem restoration efforts (UN General Assembly, 2021).
The grasslands of Eurasia are the largest continuous grassland ecosystem on Earth and are instrumental to economy, food security (O'Mara, 2012), global ecology Smelansky & Tishkov, 2012), air circulation (Lioubimtseva & Henebry, 2009), and the carbon cycle (Ni, 2002;Rolinski et al., 2021). However, grassland ecosystems in Eurasia vary considerably (De Keersmaecker et al., 2015;Smelansky & Tishkov, 2012), experience diverse climate change impacts (IPCC, 2019;Lioubimtseva & Henebry, 2009), and have different land use histories (Freitag et al., 2021;Kerven et al., 2021). Since the beginning of the 20th century, grassland ecosystems in Eurasia have sustained several extensive and intensive changes in land cover (Winkler et al., 2021) and land use intensity that have triggered long-term changes, often accelerated by climate change (Jiang et al., 2017). For example, many grasslands in Central Asia were ploughed during the Virgin Lands Campaign in the 1950s and 1960s (Karch et al., 1964;Prishchepov et al., 2020), others experience the desertification of the Circum-Aral region that began in the 1960s (Saiko & Zonn, 2000), and yet other widespread land abandonment after the collapse of the Soviet Union in 1991 (de Beurs & Henebry, 2004;Kraemer et al., 2015;Lesiv et al., 2018;Prishchepov et al., 2013). Those changes have resulted in long-term changes in grassland vegetation cover and species composition, altered fire regimes Dubinin et al., 2011;Freitag et al., 2021) and shrub encroachment (Li et al., 2015). In Mongolia grasslands are threatened by intensified and often sedentary grazing replacing mobile pastoralism, which results in degradation and desertification (Meng et al., 2021;Sanzheev et al., 2020). Similarly, in Inner Mongolia of China, current grassland species composition and ecology reflect effects of the agricultural intensification that occurred between the early 1950s and late 1990s, followed by rewilding and afforestation initiated by the Great Green Wall and Grain for Green programs in late 1970s and 1999s, respectively Song et al., 2014;Yin et al., 2018). Consequently, when analyzing changes in grassland cover in Eurasia it is necessary to recognize different ground cover trajectories of change to correctly attribute drivers, something that a change in "greenness" alone cannot capture.
Remote sensing data are well suited for grassland monitoring, including trend detection, due to their reliability and effectiveness for large-scale analysis, long-term image archives, and ability to capture even subtle changes (Ali et al., 2016;Dara et al., 2020;Reinermann et al., 2020). Vegetation indices derived based on active optical data are commonly used proxies for monitoring grassland vegetation cover, productivity and density (e.g., NDVI (Miao et al., 2021), LAI (Pasolli et al., 2015), or fPAR (Hobi et al., 2017)).
Several analyses based on time series of remotely sensed data have identified increasing and decreasing trends in grassland vegetation productivity and green vegetation cover in Eurasia often referred to as "greening" and "browning," respectively (Cortés et al., 2021;de Beurs et al., 2015;de Jong et al., 2012).
This "browning" has been attributed to climate variability (de Beurs et al., 2015;de Jong et al., 2013;Jiang et al., 2017) or the combination of climate and human activities (Zhang et al., 2021), though the effects of weather were not always statistically significant (de Beurs et al., 2009). Concurrently, grasslands in Mongolia and Eastern Mongolia are "greening," with the increase in productivity according to long-(e.g., 1982long-(e.g., -2008long-(e.g., : de Jong et al. (20121985-2015: Miao et al. (2021) and medium-term studies (e.g., : Chen et al. (2019-2016: Ding et al. (2020and : Zhang et al. (2021). The "greening" has also been attributed to climate (de Jong et al., 2013;Miao et al., 2021) and a combination of climate and anthropogenic activity (Zhang et al., 2021). Unfortunately, the different trend analyses often show different change patterns and sometimes even contradicting trends for a given region (e.g., "browning" in central and northern Mongolia in Miao et al. (2021)  The differences between identified trends may be due to differences in the timing and length of the observation records (Cortés et al., 2021) as well as due to different analytical approaches (Pan et al., 2018). Furthermore, although the long-term changes in grasslands in Eurasia have been extensively studied, most of the analyses are based on "one-dimensional" indices of photosynthetical activity (but see de Beurs et al., 2015;Hill & Guerschman, 2020), and do not account for temporal and spatial autocorrelation in the remotely sensed data. We sought to overcome both of these limitations here.
Spectral mixture analysis (SMA) quantifies all ground cover components of grassland cover (i.e., green vegetation, nonphotosynthetic vegetation, soil, and shade), which characterize grassland conditions more comprehensively compared with vegetation indices (Masiliunas et al., 2021). Furthermore, spectral unmixing is advantageous in sparse vegetation conditions where soil reflectance affects vegetation indices (Elmore et al., 2000;Huete et al., 1985;Smith et al., 2019). Cumulative Endmember Fractions (CEFs), calculated as sums of endmember fractions over a period of 1 year, or a growing season (Lewińska et al., 2020(Lewińska et al., , 2021, capture the full range of illumination conditions and phenology phases, and "normalize" those effects, which allows unbiased year-to-year comparisons. The "four-dimensional" information space of CEFs allows the identification of change trajectories (e.g., to separate a change from green vegetation to soil from a change from green vegetation to non-photosynthetic vegetation). To better capture trajectories, single fractions can be combined into ratios to illustrate transitions among fractions over time, which facilitates the mapping and understanding of change processes and their drivers.
Prior trend analyses of remotely sensed data have largely not accounted for temporal nor spatial autocorrelation leading to potentially incorrect conclusions (Ives et al., 2021). In per-pixel trend analyses, not accounting for temporal autocorrelation can lead to false positive identification of temporal trends (de Beurs & Henebry, 2004). Similarly, failing to account for spatial autocorrelation across the data (or a map), precludes recognition of true and statistically significant changes at the map scale (i.e., at the scale of the entire study area), which can result in both, false detection of clustered changes, and omission of subtle trends (Ives et al., 2021). Although the complications originating from autocorrelation in remotely sensed data are recognized (e.g., de Beurs et al., 2015;Tomaszewska et al., 2020;Zhou et al., 2001), few analyses have correctly accounted for it due to associated statistical and computational challenges (but see Bi et al., 2013;de Beurs et al., 2015;Xu et al., 2013). The Mann-Kendell test, despite frequent claim to the contrary, does not account at all for temporal autocorrelation, and the seasonal Mann-Kendall test accounts only for temporal autocorrelation among seasons, but not among years (Hirsch & Slack, 1984).
Our goal was to investigate systematic trends in ground cover composition of grassland ecosystems in Eurasia while accounting for temporal and spatial autocorrelation in remotely sensed time series using the Partitioned Autoregressive Time Series (PARTS) approach (Ives et al., 2021). Specifically, our objectives were to: (i) calculate CEFs for all grassland-specific ground cover fractions (i.e., green vegetation, non-photosynthetic vegetation, soil, and shade) based on MODIS 2002-2020 time series; (ii) identify per-pixel systematic temporal trends in grassland ground cover depicting diverse trajectories of change and evaluate the range of spatial autocorrelation of these changes; and (iii) analyze the detected trends at map-scale, testing hypotheses about the effects of different independent environmental variables and their interactions (i.e., climate zone, country, elevation, land cover class, and SPEI3 (Standardized Precipitation-Evapotranspiration aggregated over 3 months)) on ground cover composition while accounting for spatial autocorrelation.

| Study area
The Eurasian grasslands ( Figure 1) comprise 33 ecoregions from forest steppe to deserts (Olson et al., 2001) (Figure SA1), and represent 6 of 14 global biomes with a strong climate gradient from north to south ( Figure SA2). In the northern part and in the mountain regions of the Pamir and the Tian Shan, climate is continental, but that changes toward the south into cold semiarid and cold desert. The mean annual temperature varies from −23°C in the mountain regions to 20°C elsewhere. The maximum of the warmest month varies from 10° to 40°C, while the minimum of the coldest month varies from 0° to −35°C. The annual precipitation ranges from 60 to 600 mm/year, with extreme values from almost 0 mm/year in the Taklamakan and Gobi Deserts to 1000 mm/ year in the Pamir, the Tian Shan, and the Altai mountains (Beck et al., 2018;Karger et al., 2017).
The soils are mainly chernozems and dark kastanozems in the north, light-colored kastanozems in the center, and yermosols in the south ( Figure SA3). Sandy soils, along with salty soils, are common in both the Aral and the Caspian depressions, and in the Kyzylkum, the Karakum, the Taklamakan, and the Gobi deserts. Lithosols are found in the mountains. The vegetation varies from forest steppe and meadow steppe (forbs) to shrubby deserts, depending on humidity, with some enclaves of boreal forests in the mountains.
In terms of land use, the west, northwest, and far northeast, where soil and weather conditions are the most favorable, are largely ploughed. Some enclaves of croplands are also located along larger rivers in the desert parts of the Eurasian grasslands (Hankerson et al., 2019;Sulla-Menashe et al., 2019). Abandoned fields are more common toward the south of the former Soviet republics (Meyfroidt et al., 2016). Central Kazakhstan and most of Mongolia are covered by dry grasslands and used as pastures (Hankerson et al., 2019; Mongolian Statistics Office http://www.en.nso.mn/).
The deserts or mountains of the South are also used as pastures or are depopulated.

| Earth observation data
We analyzed the 2002-2020 time series of the 8-day 500 m MODIS surface spectral reflectance product (MOD09A1; Collection 6) available in Google Earth Engine (Gorelick et al., 2017;data accessed in March 2021) to calculate CEF (Lewińska et al., 2020(Lewińska et al., , 2021. The data are atmospherically corrected (Vermote et al., 2002) and adjusted for sensor degradation (Lyapustin et al., 2014). We selected only high-quality pixels, and excluded clouds and snow (based on the "StateQA" band). We did not include data acquired prior to 2002 due to their lower quality.
To identify endmembers of ground cover fractions, which are essential to derive CEFs, we analyzed 55 Landsat scenes acquired between 2003 and 2019 in 11 footprints ( Figure 1, Table SB1).
When selecting footprints we prioritized those representing the predominantly arid climate of our study area, while ensuring representation of the main ecosystems, soil types, and country-specific conditions ( Figures SA1-SA3). For each footprint, we selected cloudfree or nearly cloud-free scenes acquired during multi-year dry and wet spells ( Figure SB1, Table SB1) as well as different phenological phases. This approach ensured a broad selection of vegetation and soil conditions, hence maximizing the diversity in our endmembers' pool. Suitable winter-time scenes were not available sue to snow and cloud (Table SB1), but this did not hinder analyses due to the design of the CEF (Section 2.4). We downloaded all the Landsat data from the U.S. Geological Survey (USGS) service (data accessed in November 2020) as collection 1 tier 1 surface reflectance products that were atmospherically corrected with LEDAPS (ETM+: Masek et al., 2006) or LaSRC (OLI: Vermote et al., 2016). We removed all pixels flagged as cloud, shadow, or snow according to the pixel quality assessment bands (Zhu & Woodcock, 2012). Finally, we adjusted the OLI scenes to ETM+ surface reflectance following Roy et al. (2016).
All ancillary data were resampled to 500 m to match MODIS data using majority rule (for categorical data) and bilinear interpolation (for continuous data).

| Cumulative Endmember Fraction time series
For the SMA we identified four ground cover fractions characterizing grasslands, herein called endmembers: green vegetation, non-photosynthetic vegetation (i.e., dry leaves, shrub twigs), soil (or rock), and shade (i.e., vegetation micro-shadowing and topographic effect). Because 500-m resolution is too coarse to ensure spectrally pure pixels for endmember identification, we exploited in the similarity of spectral and orbital parameters of MODIS and Landsat Hilker et al., 2009;Zhu et al., 2019).
Following the methodology of Lewińska et al. (2020) we identified candidate "image endmembers" from Landsat images, making sure to cover a wide range of environmental and meteorological conditions across the study area (see Section 2.2, and Supplement B for details). By calculating the Minimum Noise Fraction and Pixel Purity Index for each selected Landsat scene, we identified and tested spectra of multiple soil types and green vegetation in different phenological stages (Table SB2). Additionally, we used the Ecological Spectral Information System (ECOsis) Spectral Library (https://ecosis.org; accessed in May 2019) to add the nonphotosynthetic vegetation endmember, which cannot be identified reliably from the Landsat images (Supplement B). We selected the set of final endmembers ensuring the lowest collinearity (Meer & Jia, 2012) and the lowest overall SMA RMSE, which quantifies how well endmembers explained the variance in the data (see Supplement B for details). Because our approach aims at quantifying the main ground cover characteristics of grasslands (i.e., green vegetation, non-photosynthetic vegetation, soil, and shade) and their change over time, rather than mapping the abundance of specific plant species, it allowed us to use more generically defined endmembers. The use of non-specie-specific endmembers definitions is common in SMA-based studies covering large grassland areas (Guerschman et al., 2015;Lewińska et al., 2020), and in other vegetation types (Bullock et al., 2018;Chen et al., 2021;Nill et al., 2022) because it allows to track changes in ground cover over time.
We applied the selected four endmembers in a constrained SMA model and ran it on each 8-day MODIS surface reflectance dataset available for 2002-2020 over Eurasia. Because MODIS band 5 has no equivalent in Landsat TM/ETM+, we excluded it from the SMA. We composited the resulting time series to monthly observations selecting from all corresponding 8-day datasets a set of endmembers with the highest green vegetation fraction. Finally, we calculated pixel-based endmember-specific CEFs, that is, the sums of monthly fractions aggregated for the snow-free period of each year. We identified the pixel-specific snow-free period based on the 2002-2020 MODIS data taking the median of the first and last snow-free months to define the aggregation period ( Figures SA4 and   SA5, respectively). For each year we normalized the four CEFs on a per-pixel basis to conjointly sum up to 100%, allowing for straightforward comparison among years.

| Change trajectories in grassland ground cover composition
The design of the CEFs reflects a four-dimensional feature space with each of the ground cover types being represented by one axis used to describe the ground cover conditions. Because the sum of all ground cover fractions is constrained to be 100%, this feature space can be represented by a tetrahedron where at respective apexes each ground cover fraction reaches 100% coverage ( Figure 2).
Hence, by definition a change in one ground cover fraction results in a change in at least one other ground cover fraction. This principle applies to both a change in ground cover composition between two time-steps and long-term changes. It is possible to analyze changes in individual ground cover fractions, yet such an approach misses important aspects of the changes. Instead, it is advantageous to analyze the relationship between two or more fractions captured through ratios between fractions, which greatly facilitates interpretation of change trajectories. Because statistical distributions of ratios are often highly skewed, and were in our data too, we applied a logarithmic transformation to the response variables: where e_i and e_j are CEFs for endmember i and j, respectively, in year n, and c is a small constant of 0.0001 (to prevent dividing by 0). Importantly, logarithmic transformation preserved the additive properties of our ratios allowing comparisons . We first assessed overall long-term changes in vegetation cover, which we estimated as the trend in the sum of green and non-photosynthetic vegetation over soil CEFs ( gv + npv soil ). Subsequently, we analyzed positive trends in green vegetation (gv), non-photosynthetic vegetation (npv), and soil CEFs. Since positive and negative changes are interconnected among all fractions, we focus in our results on positive changes. Finally, we investigated trends in green vegetation in relation to non-photosynthetic vegetation ( gv npv ) and soil ( gv soil ), which allowed us to distinguish areas where the balance between green and dry vegetation fractions changed, irrespective of overall trends of either revegetation or desertification. Overall, the use of ratios allows for better understanding of the changes and their drivers (Table 1). Although some of the results were visually similar (e.g., gv and gv soil ), the ratios provide a quantitative assessment of specific change pathways (e.g., change in green vegetation in relation to soil) that single ground fractions do not capture, especially when testing mapscale statistical hypotheses as we intended (see Section 2.6).
For clarity, we refer to the overall direction of changes in ground cover composition using full names of fractions, that is, "green vegetation," whereas when discussing a specific change trajectory we use the abbreviation, for example, gv. Finally, we present the shade ground cover fraction and all relevant results in the Supplement A, section SA.1.

F I G U R E 2
Conceptual 4D feature space of the Cumulative Endmember Fractions in grasslands showing change trajectories in ground cover composition captured with single fractions and their ratios. For process-based meanings of the trajectories, see Table 1.

TA B L E 1
Examples of process-based meanings of change trajectories in grassland ground cover composition.

Endmember Increase Decrease
gv Increase in vegetation productivity or vegetation cover due to natural (e.g., revegetation) and anthropogenic factors (e.g., recultivation, cultivation practices) Change in species composition Longer vegetation season Decrease in vegetation productivity or vegetation cover due to natural and anthropogenic factors (e.g., land abandonment, cultivation practices, grazing) Change in species composition Shorter vegetation season npv Desiccation due to natural (e.g., increasing aridity) and anthropogenic factors Overall increase in biomass leading to more dry matter in the senescence phase Change in species composition (e.g., shrub encroachment) Desiccation due to natural (e.g., decreasing aridity) and anthropogenic factors Overall decrease in biomass leading to less dry matter in the senescence phase Change in species composition soil Increase in soil visibility due to natural (e.g., desertification) and anthropogenic factors (e.g., introduction or change in cultivation practices) Decrease in soil visibility due to natural (e.g., revegetation) and anthropogenic factors (e.g., land abandonment, change in cultivation practices) shade Change in vegetation structure, for example, shrub encroachment Increase in soil crust Change in vegetation structure, for example, reversing of shrub encroachment Decrease in soil crust 2.6 | PARTS: Accounting for temporal and spatial autocorrelation in trends Developed by Ives et al. (2021Ives et al. ( , 2022 PARTS is a three-step statistical approach that allows for the analysis of trends in time series of raster data while accounting for temporal and spatial autocorrelation and in combination with additional covariates (e.g., meteorological data, or land cover information). First, by using an autoregressive model, PARTS derives temporally uncorrelated trends for each pixel 2.6.2 | Spatial autocorrelation in time trends: Generalized least square regression To investigate whether the spatial patterns of our trends could be explained by environmental variables, we tested map-scale statistical hypotheses with GLS regressions that accounted for spatial autocorrelation. We defined the magnitude of the spatial correlation in each trend dataset as the unexplained variation in trends among pixels where the magnitude of autocorrelation depends on the distance between pixels (i.e., pixels located closer are more similar): where cor{e i , e j } is the correlation between the random errors in pixels i and j, e i and e j , d ij is the distance between pixels, and nugget, range, and shape are parameters that determine the dependence of the correlation on distance. We calculated the range and shape of spatial correlation for each trend map, which give the distance beyond which two pixels are uncorrelated and the shape of the decay function, respectively (Ives et al., 2021(Ives et al., , 2022. We used this information in PARTS' GLS when we tested mapscale statistical hypotheses on the relation between the trends and the following explanatory variables: land cover (3 levels), climate zones (12 levels), countries (8 levels), elevation (continuous), and a temporally uncorrelated trend in SPEI3 (continuous). To examine changes in grasslands originating from regional land use legacies and policies, we further tested for statistical interactions between climate zone and country, land cover class and country, and SPEI3 trend and country factors. Due to the limited area of grasslands in some countries, we limited our country-specific analyses to China, Kazakhstan, Kyrgyzstan, Mongolia, Russia, Tajikistan, Turkmenistan, and Uzbekistan.
To overcome reduced statistical power arising from high correlation between adjacent pixels (Ives et al., 2021) we defined the magnitude and shape of spatial autocorrelation for each dataset and performed all subsequent GLS analyses using every third pixel in x and y direction of each map. Considering every third pixel aimed at preserving as much of the original data as possible, yet systematic sampling reduced data volume to about 11% of all pixels, which substantially accelerated the computations. Because remote PARTS divides the data points of the time trends in large remotely sensed data into random, non-overlapping partitions, then performs GLS regression on each partition and finally combines the statistical results from each partition (Ives et al., 2021(Ives et al., , 2022, preserving as much as much of the original data as possible is critical for sufficient represent all factors' levels within each partition. Our tests showed that using a systematic subsample of a dataset has negligible impact on differences in the estimated range and shape of spatial autocorrelation (Table SA1).

| Cumulative Endmember Fractions
We successfully derived annual 2002-2020 CEFs for grasslands in Eurasia. The four endmembers that we selected for the SMA represented well the grassland ground cover fractions for the entire study area. The pixel-specific CEF aggregation period accounted for local phenological activity and excluded spurious values from the off-season. The per-pixel mean annual cumulative RMSE aggregated for the corresponding snow-free period rarely exceeded 1% reflectance, with the standard deviation not greater than 0.5% ( Figure SA6), indicating how well our unmixing performed. The overall changes in vegetation cover were insignificant at the map scale (Table 3), but changes in vegetation cover were positively related to changes in SPEI3 (Table 4). At the map scale, the increase in vegetation cover was significant in herbaceous grasslands (Table SA3), regions under diverse variants of cold climates (Table SA4), and in Kyrgyzstan and Russia (Table SA5). Moreover, we confirmed a significant increase in vegetation cover in cold climates with dry winter and cold summers in China. We also found significant decreases in vegetation cover in cold climates without dry season in China, Kazakhstan, Kyrgyzstan and Russia, and in tundra areas in Kyrgyzstan and Russia (Table SA8). Individual factors had a significant overall effect on shaping trends in the gv + npv soil ratio, but the effects of their interactions were insignificant.

| Increase in green vegetation fraction
Positive trends in green vegetation fraction were widespread in Mongolia, China, northeastern Kazakhstan, and in the Tian Shan Mountains (Figures 4a and 5a,b). That pattern was apparent in gv , and in the gv soil ratio time series. The pixel-level positive trends in the gv npv ratio were the strongest in the Aral region (Figure 5a). The range of spatial autocorrelation in trends, which quantifies the maximum distance within which two observations are statistically dependent, was >110 km for gv, and around 50 and 33 km for gv soil and gv npv , respectively ( Table 2).
The overall map-scale positive systematic trend in green vegetation fraction was significant only for the gv soil ratio (Table 3). We found, however, a significant positive impact of elevation on gv trend (slope 0.096), and significant correlation between trends in gv and SPEI3 (0.089), and in the gv soil ratio and SPEI3 (Table 4). At the map scale, that is, across our entire study area, relations between environmental variables and spatial patterns in green vegetation trends were complex. Within the areas with herbaceous land cover we confirmed a significant increase in the gv soil ratio, and so did we in the "Transition" class in the gv fraction (0.107; Table SA3). Concomitantly, increases in green vegetation fraction were significant in regions with cold and dry winters as indicated by positive changes in gv and gv soil , as well as in cold climates without dry season according to the gv soil ratio (Table SA4). When comparing the increase in green vegetation fraction among countries, we F I G U R E 3 Systematic trends in the ratio of the sum of green vegetation and non-photosynthetic vegetation over soil ( gv + npv soil ) Cumulative Endmember Fractions. Corresponding map with the strength of the temporal autocorrelation in Figure SA11, and map showing pixel-level trends with p < .05 significance are presented in Figure SA7. detected significant but small increasing trends in gv in China and Mongolia (0.091 and 0.093, respectively; Table SA5), and in gv soil ratio in Kazakhstan, Kyrgyzstan, and in Russia. Importantly, increases in the gv soil ratio were significantly related to an increase in SPEI3 in China, Kazakhstan, Mongolia, and Russia (Table SA6), whereas the increases in the gv npv ratio were related to increase in SPEI3 in China and Mongolia. Although we found no significant increase in green vegetation fraction for the interactions of country and land cover (Table SA7), we did find significant increase in gv in cold and dry regions in China, Mongolia, and Tajikistan (Table SA8), and positive systematic changes in the gv npv ratio in arid deserts and steppes in China, Kazakhstan, Kyrgyzstan, Turkmenistan, and Uzbekistan. Importantly, we confirmed a significant effect of our environmental variables and their interactions on changes in green vegetation fraction (Chi-square results in Tables SA3-SA8). The only exceptions were the gv soil ratio for interactions between country and land cover, and the gv npv for interaction between country and climate zone.

| Increase in non-photosynthetic vegetation fraction
Increasing trends in the non-photosynthetic vegetation fraction were abundant in Eurasia at the pixel level, and especially common in Central Asia and Inner Mongolia (Figure 4b), with the latter having the strongest npv trends. An increase in the amount of nonphotosynthetic vegetation combined with decreases in photosynthetically active vegetation captured in the gv npv ratio was widespread in Central Asia, the Gobi Desert, and in the mountains (Figure 5a).
The range of spatial autocorrelation was of about 14 and 33 km for npv and gv npv , trends, respectively (Table 2). Across Eurasia, we found a strong and significant overall mapscale increase in the npv fraction (Table 3; trend slope: 0.246).
Moreover, both elevation and SPEI3 had positive relationships with changes in the npv fraction (Table 4; 0.247 and 0.248, respectively). The map-scale hypothesis testing confirmed the increase in non-photosynthetic fraction was significant across many factors' levels (Tables SA3-SA8). For example, the increase in npv was significant in the "transition" and herbaceous land cover, as was the decrease in the gv npv ratio in barren and "transition" classes (Table SA3). Increasing trends in npv were significant and strong in arid deserts and steppe, in polar climates, and in areas with cold climate with no dry season, or dry summer (Table SA4). Conversely, decreasing trends in the gv npv ratio were significant only in arid deserts (Table SA4). Our map-scale analyses confirmed significant positive trends in npv in all countries except Mongolia (Table SA5).
Interestingly, in Russia the increase in npv was significantly related to negative changes in SPEI3, whereas in Uzbekistan increasing npv was associated with increasing SPEI3 (Table SA6). Negative trends in the gv npv ratio were limited to China (Table SA5), where they were also significantly related with SPEI3 increases (Table SA6). We found significant increasing trends in npv for interactions between TA B L E 2 Range of spatial autocorrelation (r) in trends in change trajectories in grassland ground cover composition. Note: Abbreviation as in Table 3. country and land cover only in "transition" and herbaceous classes in Tajikistan and Turkmenistan (Table SA7). Map-scale npv trends for the interactions between country and climate were significant only in arid deserts in China, cold regions with dry and cold summers in Kyrgyzstan, and in temperate climates with dry and hot summer in Turkmenistan (Table SA8). A Chi-square test confirmed insignificant differences between npv trends among countries, and the combination of country and climate (for npv and gv npv ).

F I G U R E 4 Systematic trends in (a) green vegetation (gv), (b) non-photosynthetic vegetation (npv), and (c) soil (soil) Cumulative
Endmember Fractions. The slope is reported in percentage points. Corresponding maps with the strength of temporal autocorrelation in Figures SA11-SA13, respectively, and maps showing pixel-level trends with p < .05 significance in Figure SA7. Area statics for negative and positive trends are presented in Table SA2.

| Decrease in vegetation cover
Increases in pixel-level abundance of soil fraction were most widespread in Central Asia, Caspian Sea Depression, and the Gobi Desert (Figures 4c and 5b). While positive trends in the soil fraction were often clustered, changes in the gv soil ratio were more dispersed, which was reflected in the range of spatial autocorrelation identified at around 12 and 50 km, respectively (Table 2). We found no overall map-scale significant increase in the soil fraction abundance (Table 3), but higher elevation and SPEI3 were inversely related to the increase in soil fraction at the map scale (Table 4; trend slope of −0.206 and −0.276, respectively).
The increase in the soil fraction was insignificant at the map scale (Tables SA3-SA8). The only exception occurred in regions with a cold climate and dry and hot summers in Kazakhstan (for soil) and with a cold climate with no dry season and cold summer in Russia ( gv soil ; Table SA8). Notably, all environmental factors had a significant impact on trends in soil and gv soil when inspected separately. However, for factor's interaction, we found significant impacts only for interactions between country and land cover for soil and between country and climate for gv soil (Tables SA7 and SA8, respectively).

| DISCUSS ION
We analyzed systematic trends in grasslands ground cover composition in Eurasia between 2002 and 2020 using ground cover fractions and statistical methods that directly account for temporal and spatial autocorrelation in the remotely sensed data.
Across Eurasia, we found a statistically significant increase in  Figures SA14 and SA15, respectively, and maps showing pixel-level trends with p < .05 significance in Figure SA8. no statistical evidence for a map-level increase in the soil fraction across Eurasia. This shows that over the past 20 years there was an overall increase in vegetation cover in the grasslands of Eurasia, with grasslands now comprising greater non-photosynthetic fraction due to an overall increase in vegetation matter and potential shifts in vegetation type.

| Increase in vegetation cover
Our results demonstrated interesting differences in the change pathways of vegetation cover when analyzed as green vegetation, non-photosynthetic vegetation, and the sum of both. Through hypothesis testing we rejected the hypothesis of an overall maplevel increase in vegetation cover ( gv + npv soil ) and gv fraction, but confirmed significant increase in non-photosynthetic vegetation (npv) and increase in green vegetation paired with a decrease in soil ( gv soil ) (Table 3). This shows the importance of disentangling green and non-photosynthetic vegetation fractions and analyzing their specific change pathways separately. Although both fractions are closely related, and even though both together represent vegetation cover, the relation between gv and npv can alter due to successional changes in vegetation type or management practices (e.g., due to restoration programs like in China; Chen et al., 2015;Song et al., 2014;Yin et al., 2018). Consequently, we do not interpret the increase in non-photosynthetic vegetation, often detected as  (Table 4), which suggests that both fractions were similarly affected by meteorological conditions. The observed positive trends in vegetation-related fractions and their relation to SPEI3 may reflect a return to normal and wet weather conditions after the prolonged dry period in 1990s and 2000s (Guo et al., 2018;Sheffield et al., 2009).

Pixel-level increases in vegetation were common across
Eurasia. Our trend maps showed strong gv + npv soil increases in northern Kazakhstan, northern Mongolia, northeast China, and Inner Mongolia (Figure 3). When compared to other trend maps, the increase in vegetation was similar pixel-level positive trends in gv soil and the increase in npv (Figures 4b and 5). Positive trends in gv, concentrated in the northern parts of our study area and in the mountains, also added to the overall pattern, though their contribution was less pronounced.
The overall footprint of our vegetation increase-related pathways corresponded well with the pixel-based "greening" patterns identified by de Jong et al. (2011de Jong et al. ( , 2013, Jiang et al. (2017), Munier et al. (2018, Zhang et al. (2021), and Chen et al. (2019), but that greening was less pronounced in our maps. The difference is most likely due to our method accounting for temporal autocorrelations ( Figures SA10, SA11, SA12, and SA15), but could also be due to differences in metrics and study periods. Explicitly, the pixel-level increase in the gv soil change trajectory in northeastern Kazakhstan corresponds with "greening" attributed to land abandonment (de Beurs & Henebry, 2004;Robinson, 2016), whereas "greening" in Mongolia, Inner Mongolia, and northeast China has been associated with beneficial combination of grazing and precipitation (Miao et al., 2021;Zhang et al., 2021) and land use changes arising from the Grain for Green reforms Song et al., 2014;Yin et al., 2018).
At the same time, the increase in the non-photosynthetic vegetation fraction, which was the most widespread change we detected at the pixel level (Figure 3), agreed with many areas previously identified as "browning" hotspots (de Beurs et al., 2009(de Beurs et al., , 2015de Jong et al., 2012de Jong et al., , 2013Ives et al., 2021;Jiang et al., 2017). Some differences among the patterns we attributed to a different "browning" definition, monitoring approach, length of the time series, and treatment of temporal autocorrelation (Figures SA12 and SA14).
However, not all of the pixel-level trends were significant at the map level. For example, the increase in gv + npv soil ratio was significant only in Kyrgyzstan and Russia, gv was in China and Mongolia, whereas increasing trends in gv soil were significant in Kazakhstan, Kyrgyzstan, and Russia (Table SA4). Concomitantly, the increase in npv was significant, strong, and uniform in all countries expect Mongolia (Table SA5). Moreover, in China, Kazakhstan Mongolia, and Russia increasing trends in gv soil matched positive changes in SPEI3 (Table SA6), which complements results from Jiang et al. (2017), Zhang et al. (2021) and Miao et al. (2021), and is likely due to vegetation recovery after severe droughts in late 1990s and early 2000s (Guo et al., 2018;Sheffield et al., 2009). Significant map-scale increases in photosynthetic activity in regions under cold climates align with easing of temperature-related vegetation growth limitations and a shift toward warmer and wetter conditions (Zhang et al., 2021) and with shorter snow seasons in the mountains (Tomaszewska & Henebry, 2018).
The distinction between the increase in green vegetation and non-photosynthetic vegetation was possible through the gv npv ratio (Figure 5a). Although we found neither an overall map-level "greening" nor "browning" in Eurasia (Tables 3 and 4), map-level "browning" was statistically significant in barren areas, transitional grasslands (Table SA2), and arid cold deserts (Table SA4). Moreover, in China we found a surplus of non-photosynthetic vegetation over green vegetation (Table SA5), which was related to meteorological conditions (Table SA6), and may reflect a change in vegetation composition, such as shrub encroachment or afforestation, which have occurred in Central Asia, Northwest China (Li et al., 2015), the Tian Shan Mountains (Zhumanova et al., 2021), and Inner Mongolia (Bardgett et al., 2021;Hu & Nacun, 2018). At the pixel level, we found particularly good alignment of negative trends in the gv npv ratio with "browning" detected during the last two decades, which could be due to suggested turning points in vegetation response from "greening" to "browning" (Horion et al., 2016;Li et al., 2015;Pan et al., 2018), or vegetation recovery after persistent dry spells in the late 1990s and early 2000s (Guo et al., 2018;Sheffield et al., 2009). Furthermore, our results are similar to the increase in non-photosynthetic vegetation mapped by Hill and Guerschman (2020). Strong interannual autocorrelation in gv npv time series indicates a gradual process, typically associated with accumulation of dry vegetation mater and change in meteorological conditions.

| Decrease in vegetation cover
Prior studies suggested widespread grassland degradation and desertification in Eurasia (Cai et al., 2022;Hu et al., 2020;Zhang et al., 2018). However, our hypothesis tests found no significant map-scale decrease in vegetation combined with an increase in the soil fraction (Table 3) even though our trend maps showed pixellevel decrease in vegetation cover in favor of soil in Central Asia ( gv + npv soil , soil and gv soil ), and the Gobi Desert ( gv soil ) (Figures 3, 4c, and 5b). Statistical tests showed that these trends were significant only in areas with cold climates in Kazakhstan, and in regions under cold climate with hot summers and without dry season in China, Russia, and Kyrgyzstan (Table SA8), which are typical for high mountains, such as the Altai Mountains and Caucasus.
Our pixel-level trend maps resembled to some extent the "browning" hotspots reported by Jiang et al. (2017), the decrease in vegetation indices and increase in Tasseled Cap brightness demonstrated by de Beurs et al. (2015), negative NDVI trends after 1998 observed by Li et al. (2015), and increases in bare soil mapped by Hill and Guerschman (2020). Overall, we expected to find only limited coincidence between our soil-specific results and generic "browning" maps, because most of the previous research considered "browning" as an overall decrease in vegetation indices, not an increase in soil presence, which is relatively sparse in Eurasia (Hu et al., 2020). As such, the detected pixel-based changes in soil, gv soil and gv + npv soil fractions corresponded best with desertification in Central Asia and Kazakhstan arising from agricultural exploitation and climate change (Cai et al., 2022;Hu et al., 2020), and salinization (Robinson, 2016;UNEP, 2011). The high inter-annual autocorrelation in the detected increase in soil fraction indicates the steadiness of desertification processes, presumably propelled by rising temperatures and aridity (Li et al., 2015). Furthermore, the small range of spatial autocorrelation and well-defined footprints of areas showing an increase in soil fraction may suggest land use drivers, such as localized overgrazing  and chemical degradation (Robinson, 2016

| Significance and limitations of the approach
We identified systematic changes in grasslands in Eurasia by applying a new and statistically rigorous approach. First, our approach is based on spectral unmixing, which provides physically based information on four grassland-specific ground cover fractions derived from multiple spectral bands. This is an advancement compared with vegetation indices and allowed us to better analyze main characteristics of grassland ecosystems and their changes (Masiliunas et al., 2021). However, we acknowledge the potential limitations of SMA in separating soil and non-photosynthetic vegetation fractions (Kowalski et al., 2022) apparent in higher RMSE in regions with sparse vegetation cover ( Figure SA6).
Second, by combining multiple ground cover fractions into ratios, we were able to identify distinct change trajectories better reflecting change processes. For example, ratio-based analyses permitted us to identify desertification-related changes constituting a decrease in vegetation accompanied by an increase in soil fraction.
However, since various processes can result in the same change trajectory final process attribution still requires knowledge on local land use changes.
Third, by accounting for temporal autocorrelation at the pixel level our results give reliable pixel-based trends in ground cover fractions and their ratios. We found that temporal autocorrelation was very common in our time series (Figures SA10-SA15). When comparing the locations of autocorrelation hotspots in our data with the locations of long-term changes reported previously we noted a troublingly high degree of correspondence suggesting that some of the previously reported trends and hotspots may not be significant if they did not account for temporal autocorrelation.
Fourth, the range of spatial autocorrelation in the data varied between the datasets but was estimated consistently regardless of the spatial density of sample points ( Table 2, Table SA1). This confirms the persistent character of spatial autocorrelation, the need to account for it, and highlights the spatial variability of different processes (e.g., wide-scale increase in gv vs. small-scale increase in soil fractions). Although the uniform range of spatial autocorrelation is a simplification, accounting for it in the PARTS' GLS regression allowed us to examine different map-scale statistical hypotheses.
This is an important advancement because accounting for autocorrelation in the data both lends confidence to our results and demonstrates how the significance of the changes depends on the level of stratification and the hypothesis being tested (e.g., country vs. combination of country and climate zone).
Fifth, the PARTS method was designed specifically for regression-type analyses, identifying time trends in the response variable associated with explanatory variables. Therefore, PARTS makes it possible to predict changes in the response variable if there are changes in the future values of the environmental variables.
PARTS, however, does not assess the unexplained variations in the data. Specifically, it does not assess whether areas of, for example, "browning" trends seen on a map will continue if they are not explained by environmental variables used in the regression analysis.
Finally, we focused in our presentation on the most commonly studied changes and change trajectories, and discussed only a subset of the results stemming from our analyses to maintain clarity of the manuscript. However, our approach provides opportunity to explore other change trajectories, such as changes in the shade fraction, which may indicate shrub encroachment and development of soil crust (Chen et al., 2005) (Supplement A).

| CON CLUS IONS
Our trend analyses showed that there was no significant decrease in vegetation cover in Eurasia over the past 20 years. On the contrary, the increase in non-photosynthetic vegetation fraction was by far the most widespread and statistically significant trend across Eurasia, and in all countries expect Mongolia. Increases in green vegetation were more localized and moderate, but mostly significantly related to changes in meteorological conditions. Our results corroborate previous studies, especially those analyzing changes in Eurasian grasslands in the past two decades, but expand further upon the direction and character of changes, as well as the underlying process. Importantly, our state-of-the-art statistical approach accounting for temporal and spatial autocorrelation in the data ensures statistically robust assessment of the trends in the region. This is of great importance for understanding the changes and correct attribution of their drivers, especially considering recently arising ambiguity around trends analyses (Cortés et al., 2021;Ives et al., 2021).
Furthermore, we tested a wide range of statistical hypotheses about the effect of different variables on changes in grasslands and found that all variables affected grasslands in Eurasia, but the map-scale significance of responses and relations varied greatly. This complexity highlights the wide range of change drivers and processes in grassland in Eurasia and thus the importance of an analysis design that allows for hypothesis testing. Lastly, our approach is scalable and transferable to other time series of satellite data and regions, and can be implemented in any computational environment, assuring accessibility and reproducibility.

ACK N OWLED G M ENTS
We gratefully acknowledge support from the Advanced Information Systems Technology (AIST) program, grants NASAAIST-80NSSC20K0282, and the Land Cover and Land Use Change (LCLUC) Program of the National Aeronautic Space Administration (NASA), grants 80NSSC18K0316 and 80NSSC18K0343. We thank two anonymous reviewers for their constructive comments, which allowed us to greatly improve the manuscript.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare that they have no known competing financial interests or personal relationships that could have influence the work reported in this paper.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available in Dryads as https://doi.org/10.5061/dryad.hmgqn k9nt and at https:// silvis.forest.wisc.edu/data/euras ia-trends. These data were derived from the MOD09A1 (Collection 6) data available in the public domain (http://doi.org/10.5067/MODIS/ MOD09 A1.006).