Prominent vegetation greening in spring and autumn across China during the 1981–2018 period

Vegetation greening in China has been extensively examined, but little is known about its seasonal characteristics and its association with soil moisture (SM) and temperature changes. Using high-resolution (0.1°, 8 d) datasets of leaf area index (LAI), together with SM, soil temperature (ST) datasets, and the dominance analysis method, this study is designed to detect seasonal vegetation changes across China during 1981–2018 and its links to climate change. The results show that 56.8% of land area across China experienced a greening trend while 6.6% browning trend through 1981–2018. LAI increasing area expanded to a maximum of 59.3% in spring and the decreasing area reached a maximum of 10.6% in autumn. Spring increases in LAI in main vegetation regions were significantly correlated with positive ST anomalies, while autumn decreases in LAI except sparsely vegetated regions were correlated with negative SM anomalies. Combined SM and temperature anomalies explain 10.9% of the observed LAI changes, which is 4 times larger than that directly explained by precipitation and surface air temperature (2.7%). The warming of soil under climate change was driving the LAI increases, while drying was largely responsible for LAI decreases. These findings provide further evidence of climate change impacts on regional ecosystems and highlight the importance of soil heat and water conditions in translating global warming signals.


Introduction
Vegetation response is a direct indicator of climate change impacts on terrestrial environments and ecosystem services [1]. Increasing evidence is pointing to the reality of a rapidly changing living environment globally and regionally [2]. As global warming continues to intensify, it is vital to assess its emerging impacts on ecosystems and human habitats for supporting adaptation and policy-making [3,4].
Globally, the observed vegetation greening trends are believed to be driven mainly by CO 2 fertilization effects [5] and the browning trends are a result of increased drought stress [8]. Global warming is also expected to extend the growing season length [9] and intensify the hydrological cycle [10,11]. The risks of drought are also increasing as climate warms [12,13]. Detailed vegetation response to climate change is complex and dependent on both global and regional factors. Based on offline model simulations, a previous study suggests that warming is responsible for 8%-20% of the increases in leaf area index (LAI) data detected over 28.4% of terrestrial vegetated areas and CO 2 fertilization effects can explain 70% of the greening trends [5]. However, coupled Earth system model simulations show that rising surface air temperatures play a dominant role in vegetation greening [14]. Considerable uncertainties remain in current climate model simulations because of deficiencies in representing vegetation and land surface processes and the lack of observational data (e.g. LAI) at a high spatiotemporal resolution to constrain model simulations [15,16]. Data-driven statistical models still play an important role in understanding climate-vegetation interactions. Based on atmospheric variables such as precipitation (Pr) or potential evaporation, statistical analyses have indicated that droughts regulate vegetation browning on both regional and global scales [17,18]. According to the soil-plant-atmosphere continuum theory by Philip [19], soil plays a bridging role between the atmosphere and vegetation in terms of water cycles. Statistical analysis based on atmospheric-only variables and soil conditions can sometimes come to drastically contrasting conclusions. Taking areal changes in global drylands as an example, based on atmospheric indicators, drylands are expected to expand by 23% at the end of the 21st century relative to the area during the past 30 years [20]. In terms of soil moisture (SM), the predicted expansion is only 11% [21]. When LAI and transpiration are considered, there are virtually no expected changes (−0.7%) in dryland areas [22]. Arguably, soil plays a pivotal role in translating climate change into vegetation responses.
Previous studies have mostly focused on annual mean vegetation trends and large-scale global driving factors. This study is designed to look into seasonal characteristics of vegetation response to climate change in China by considering the soil-plant-atmosphere continuum. We consider two questions: (a) what are the spatiotemporal characteristics of seasonal vegetation changes across China during the recent decades? (b) How and to what extent are vegetation changes associated with anomalies (with 38 year means removed on the 8 d scale) in SM and soil temperature (ST)? Section 2 describes data and method; section 3 presents the main results; with a discussion on potential caveats placed in section 4 and conclusions in section 5.

Data collection and fidelity
The LAI data were obtained from the Global Land Surface Satellite (GLASS) product suite, version V50 [23]. This LAI product was synthesized from Advanced Very High-Resolution Radiometer, Moderate-resolution Imaging Spectroradiometer (MODIS) data, using a method of general regression neural networks [24]. The LAI product covers the period 1981-2018 with an 8 d interval and 0.05 • resolution. Extensive validation against in situ measurements and intercomparisons with other satellite products have shown high quality and accuracy of the GLASS LAI product, despite discrepancies due to varying algorithms and original remote-sensing data used [5,[25][26][27]. To date, this product has been widely used to examine vegetation changes and their interactions with climate change [2,5,28].
SM and ST from the European Centre for Medium-range Weather Forecasts (ECMWF) reanalysis version 5 (ERA5)-Land dataset [29] are used to examine LAI relationship with soil conditions. ERA5-Land is a replay of the land component of the ERA5 reanalysis system. The Hydrology-Tiled ECMWF Scheme for Surface Exchanges over Land (H-TESSEL) [30] was rerun at a 0.1 • spatial resolution. The soil scheme in H-TESSEL includes four vertical layers at depths of 0.07, 0.28, 1.0, and 2.89 m below the surface. The rerun simulations were not directly influenced by SM or ST observations, despite the indirect signals resulting from the assimilation of remote sensing-retrieved SM data in ERA5 [31].
The combined MODIS land-use product was used to classify vegetation types. Based on the consistent types during 2001-2010, the respective land cover types were determined according to the International Geosphere-Biosphere Programme classification scheme. The LAI, SM, and ST data are re-grided to a 0.1 • resolution using bilinear interpolation. Monthly LAI data are aggregated by average. If an 8 d data point overlaps two months, it is weighted by the number of days falling into the given month.

Linear-trend estimations of LAI
Linear LAI trends are estimated using the nonparametric Theil-Sen estimator, a method developed for robust linear regression, which computes the final slope as the median of all slopes between paired values [32]. Compared to the ordinary least-squares estimator, the Theil-Sen estimator is robust against outliers. Here, the significance of linear LAI trends was tested using the Mann-Kendall test [33].

Relative contributions of SM and ST anomalies to LAI changes
The relative contributions of SM and ST anomalies to the changes in LAI were quantified by the method of dominance analysis (DA) [34,35]. DA measures the contributions of independent variables in the regression model via the increase in the variance ratio produced by adding the variable to the regression model. By establishing a hierarchy of variable contribution levels in a pairwise way, DA measures the direct effect of an independent variable, its total effect together with other variables, and its partial effect considered together with various combinations of variables, thus allowing direct comparison of the contributions of variables in a model. When the input variables are correlated, DA uses hierarchical linear modeling instead of multiple linear regression to meet the assumption of variable independence [36]. Therefore, DA provides a robust and straightforward way to interpret the results and has been widely used to measure the relative contribution of a specific variable to the variability of air quality, flood, and groundwater storage in various regions of the world [37][38][39][40][41]. Figure 1 shows the mean states of hydrological and ecological backgrounds in China depicted in terms of annual mean Pr, land cover, SM, and LAI distributions. Hydrological regimes in China range from the arid northwest to the humid southeast, with annual Pr from less than 50 mm to over 2000 mm. SM largely follows the distribution of Pr from northwest to southeast, with distinct regional characteristics under the influence of terrestrial character and ecosystems. For instance, the spatial structure of LAI is largely regulated by land-cover types, which are closely associated with hydrological regimes [42].

Spatial structures of monthly LAI variations across China
From 1981 to 2018, LAI increases dominated the year-round variations, covering an average of approximately 56.8% of the total land area of China (figure 2). Remarkable LAI increases (with maximum LAI increases topping 1 m 2 m −2 over 38 years) were exhibited during the main growing season (April-September) spanning grassland, cropland, and savanna regions from northwest to southeast China. The overall LAI increases reached the maximum affected area (59.3%) in May and the minimum area (10.6%) in November (supplementary figure S1). In cold months (October-March), LAI increases were mainly supported by cropland and savanna regions. Note that LAI increases reached the maximum affected area in spring.
In contrast, in approximately 6.6% of the land area, the annual mean LAI decreased during the same period, prevailing in the transition zones between dry and wet climates from northeast to southwest China, along with parts of humid southeast China. In North China and southwestern China, these decreasing trends were quite strong (up to −1 m 2 m −2 over 38 years). The area with significant decreasing trends reached its minimum in July (3.8% of the land area of China) and its maximum in October (11.1%, supplementary figure S1). That is, large-area LAI reductions occurred in autumn.

Monthly LAI changes across six types of vegetation
The interannual shifts in LAI from 1981 to 2018 (figure 3, supplementary figures S2 and S3) indicate distinct changes across the six types of vegetation on the monthly scale. Overall, in the case of increasing LAI trends (figures 3(a)-(f)), the increases occurred in the main growing season in sparsely vegetated regions (Spars. Veg.), grasslands, and croplands, while seasonal increases were relatively even year-round in the savanna, deciduous forest (Decid. forest), and evergreen forest (Everg. forest) regions. Specifically, in sparsely vegetated regions, the maximum LAI increase (0.2 m 2 m −2 over 38 years) appeared in summer, as it also did in grasslands and croplands, with absolute LAI increases up to 0.5 and 0.9 m 2 m −2 , respectively (supplementary figure S2). However, the maximum LAI increases occurred in spring and summer in savanna regions and in spring in both forest type regions, with increases up to 0.5, 0.6, and 0.4 m 2 m −2 , respectively (supplementary figure S2). The percentages of the maximum increases relative to their 1981-2018 means (figures 3(m)-(r)) occurred in summer up to 56.9% in sparsely vegetated regions and 48.3% in grasslands. However, the maximum increases occurred in spring in croplands, savannas, and the two types of forest regions up to 45.5%, 25.6%, 34.8%, and 11.8%, respectively (supplementary figure S2). The relative increases arguably tended to quicken in spring.
Moreover, in the case of decreasing LAI trends, both the absolute (figures 3(g)-(i)) and relative variations (figures 3(s)-(x)) show similarly intensified decreases in the early and late periods of the growing seasons, namely, in spring and autumn. The decreases are −9.0% and −13.4% on average in spring and autumn, respectively, but the value is −4.6% in summer (supplementary figure S3).
Given that the LAI changes tended to intensify in spring and autumn from 1981 to 2018, we compared the LAI increases in these two seasons between the 2010s (2011-2018) and the 1980s (1981-1988). The probability density functions for the percentages of the two season LAI anomalies (figure 4), relative to the 1981-2018 means, show that LAI increased markedly in the 2010s compared to the 1980s, with an increment of 16.7% points on average. The probabilities of strong LAI increases (e.g. >25%) occurring were higher in the 2010s than in the 1980s, suggesting that    both the trends of and variances in monthly LAI were intensified over these 38 years.
To compare spring and autumn LAI increments between the first and last 10 year periods, we computed the incremental differences in LAI using linear regressions on a higher (8 d) temporal scale (supplementary figure S4). The maximum relative increments of the five types of vegetation appeared in spring or autumn (or even in winter for the two forest types), while only croplands saw the maximum increment in summer. The LAI increments reach 14.1, 18.9, 11.7, 12.1, 17.3, and 12.3 percentage points for sparsely vegetated areas, grasslands, croplands, savannas, deciduous forests, and evergreen forests, respectively. These relative increments further indicate that vegetation changes were more intensified in spring and autumn than in summer.

Associations of LAI changes with SM and ST anomalies
For the six analyzed vegetation types, the linear trends of underlying SM anomalies generally tended to be drying throughout 1981-2018 (figures 5(a)-(i)) following both LAI increases and decreases, except wetting in sparsely vegetated regions (northwestern China). However, the annual mean anomalies were extremely small (-0.005 and 0.002 m 3 m -3 over 38 years on average for drying and wetting regions, respectively). On the seasonal scale, however, the averaged SM anomaly reached 0.02 m 3 m -3 , thus enhancing the effects on LAI changes.
Moreover, the ST linear trends rose by 0.25 • C and 0.21 • C in the cases of LAI increases and decreases, respectively, from 1981 to 2018 on average across the regions comprising the six vegetation types. However, these trends, in turn, increased by 1.53 • C and 1.29 • C in spring, higher than the summertime (1.30 • C and 1.26 • C) and autumn (0.87 • C and 0.81 • C) trends.
To disentangle the associations between anomalies in LAI, SM, and ST, their correlation coefficients were computed on the monthly scale ( figure 6). As LAI increased, the correlations of LAI were negative with SM but positive with ST and were almost entirely significant (p < 0.05), especially for grasslands, croplands, savannas, and deciduous forests. On the seasonal scale, the summer correlations were higher than those in spring and autumn. As LAI decreased, the signs of the correlation coefficients were reversed, i.e. there were negative correlations of LAI with ST but positive correlations of LAI with SM. The autumn correlations were higher than those in spring and summer. The significant inverse correlations between SM and ST point to tight water and heat coupling processes in soils. Accordingly, SM and heat exert important but varying effects on LAI greening and browning trends.
To quantify the contributions of SM and ST anomalies to LAI changes, dominance analyses were conducted based on regression models (figure 7). The estimates show that SM and ST jointly contributed 10.9% of the total LAI trends during 1981-2018. For increasing LAI trends, SM and ST anomalies accounted for 7.8%, of which ST alone contributed 5.1%. For decreasing LAI trends, SM and ST jointly contributed 13.9%, of which SM alone contributed 10.2%. On the seasonal scale, the joint contribution of the soil variables to increasing/decreasing LAI trends reached 13.5%/11.4%, of which ST/SM alone contributed 9.3%/9.6%, on average. Therefore, ST dominates the joint contribution to increasing LAI trends, while SM dominates the joint contribution to decreasing LAI trends. Furthermore, the ratios of contribution to LAI increases by ST to those by SM reached 2.2, 3.1, and 0.9 in spring, summer, and autumn, respectively. In contrast, the ratios of contribution to LAI decreases by SM to those by ST reached 4.6, 2.9, and 10.0 in spring, summer, and autumn, respectively. In general, the contribution of soil conditions to LAI  increases was considerable in spring, although less than in summer. Soil condition contributions to LAI decreases were substantially higher in autumn than in summer.
The same dominance analyses were also carried out using Pr and surface air temperature (T2m) anomalies (figure S5). The joint contribution (2.7%) by Pr and T2m to LAI trends was only 24.8% of that (10.9%) by SM and ST. In the case of LAI increases, the joint contribution by Pr and T2m accounted for 50% of that by SM and ST. In the case of LAI decreases, the joint contribution by Pr and T2m accounted for 10.1% of that by SM and ST. The direct contributions of Pr and T2m anomalies to the LAI trends are much lower than those of SM and ST.

Discussion
The quickened vegetation changes in spring and autumn under global warming are associated with the sensitivity of photosynthetic enzymes to rising temperatures [43]. Enzyme activity increases with temperature until reaching a thermal optimum [44]. Over the past four decades, the surface air temperatures over northern hemisphere continents have risen faster in spring and autumn than in summer [45]. Warming in early spring and late autumn has greatly eased the constraint of temperature on vegetation germination and growth [2,7]. However, photosynthesis decreases as ambient temperatures approach or exceed the thermal optimum (e.g. 32 • C), slowing relative increases in summer LAI. The seasonal feedbacks of vegetation to climate warming contribute to the dumbbell-shaped relative annual variations in LAI during 1981-2018 (figures 3 (p)-(r) and (v)-(x)). Simultaneously, Pr increases alleviate water stress and enhance thermal effects on vegetation development in spring, especially in dry western China [46]. Our findings of seasonal changes in LAI are consistent with records indicating advancing spring phenophases and delaying autumn phenophases in China [47].
According to observations, the spring and autumn browning trends that occurred in northwest China were mainly the result of water stress [48]. This finding is qualitatively aligned with our results (figure 5); that is, vegetation greening changes are more closely related to ST anomalies, while browning changes are more closely related to SM anomalies. Furthermore, North China and southwestern China, which have experienced vegetation browning changes, have had major droughts during past decades [49][50][51].
In the case of LAI increases (figure 6 upper panel), the positive correlations of LAI with ST anomalies but the negative correlations with SM anomalies were generally caused by the processes in which rising ST leads to LAI increases and then enhances evapotranspiration, thus resulting in SM decreases. Conversely, in the case of LAI decreases (figure 6 lower panel), SM decreases drive concurrent decreases in LAI, which results in opposite changes in ST for the reduced latent heat loss from the soil. Although the relationships between the processes in soil-plant systems are not simply linear, close connections between soil conditions and vegetation have also been found in middle latitudes [52][53][54].
The soil-vegetation system is composed of living organisms and thus complicates disentangling the respective contributions of soil condition anomalies to LAI changes. Nevertheless, the statistical modelbased DA quantified that approximately 7.8% of vegetation greening trends could be attributed to SM and ST anomalies, 50% higher than that by Pr and T2m. Soil condition anomalies are much more directly associated with vegetation changes than with atmospheric variable anomalies. Based on a vegetation model, Zhu et al [5] estimated that 8% of global vegetation greening during 1982-2009 can be explained by climate change. Considering the regionally varying impacts of climate change [2,5], our result is generally comparable to the estimates by process-based vegetation models. In addition, the other 70% of the greening trend can be attributed to CO 2 effects, 9% to nitrogen deposition, and 4% to land-use change, according to Zhu et al [5].
Notably, to examine vegetation changes in sparsely vegetated regions, we excluded nonvegetation areas using an annual LAI < 0.05 threshold, which is lower than the 0.08 used by Liu and Xue [55] and the barren lands used by Chen et al [7]. This threshold may expand the area affected by climate change. In addition, the ERA5 model system has yet to incorporate the representations of groundwater, runoff routing, and irrigation processes, as well as dynamic vegetation activities [29]. The model inadequacies may add uncertainty to the quantification of regional soil-vegetation interactions.

Conclusions
In this paper, we report observed seasonal characteristics of vegetation changes across China and the associations of vegetation changes with anomalous soil conditions during the past 38 years. The main conclusions are: (a) the recent relative changes in vegetation tended to intensify in spring and autumn, reflecting a lengthening of the growing season and land hydrology changes under the influence of global warming. Spatially, the maximum LAI increases in spring occurred in four of six main vegetation types, except sparsely vegetated areas and grasslands reached their maximum in summer. In contrast, the maximum LAI decreases in autumn prevail in the transition zones between dry and wet climates from northeast to southwest China. (b) Observed vegetation changes were significantly correlated with soil condition anomalies, highlighting the importance of soil and land hydrology in translating climate change signals into vegetation and the wide terrestrial ecosystems.
This study has added further evidence that anthropogenic climate change is rapidly changing the terrestrial environment. The reported close connection between seasonal vegetation trends and anomalous soil conditions highlights the significance of considering the climate-soil-vegetation continuum as a whole, instead of only atmospheric variable-based indices, for assessing climate change impacts on the terrestrial ecosystems.

Data availability statement
All data used in this study are free and publicly available. LAI data are available through the Global Land Surface Satellite (GLASS) product suite (www. glass.umd.edu/). Soil moisture and soil temperature data are available through the European Centre for Medium-Range Weather Forecasts (www.ecmwf. int/en/era5-land). The Terra and Aqua combined MODIS land cover product (version 6, MCD12C1) was distributed by the Land Processes Distributed Active Archive Center (LP DAAC) of the National Aeronautics and Space Administration (NASA), U.S.A. (https://lpdaac.usgs.gov/products/ mcd12c1v006).
No new data were created or analysed in this study.