Introduction

Soil moisture (SM), as a key component of terrestrial water cycle, controls many processes and feedbacks in the climate system1,2,3. It not only plays a crucial role in the water and energy cycles between the land and the atmosphere, but is also associated with terrestrial ecosystems and thus influences carbon exchange4,5. Under the combined effects of global warming and human activities, soil moisture has shown a decreasing trend, which has triggered a water crisis1,6,7,8. Especially in the summer of the Northern Hemisphere, soil moisture deficit leads to weakening of the evaporative cooling effect, which triggers more extreme events like high temperature, drought, and both together, causing serious damage to agriculture, society, and ecosystems9,10,11,12. Therefore, it is urgent to explore the main drivers of summer soil moisture variability in the Northern Hemisphere.

Based on the water balance principle, the main factors affecting terrestrial SM are evapotranspiration (ET, including evaporation from bare soil and transpiration from vegetation), precipitation, surface runoff, and groundwater recharge (which is negligible compared to the first two). Transpiration of vegetation accounts for about 65% of ET, which is the main contributor to ET13,14. So, the dynamics of vegetation is closely related to the water cycle15. Evidence based on remote sensing and ground-based observations in recent years has demonstrated widespread spring greening in the Northern Hemisphere under global warming, which can be reflected by changes in leaf area index (LAI)9,16,17. Widespread greening alters bio-geophysical processes between the land and the atmosphere, and affects the local and larger-scale climate system. Lian et al.9 explored the role of spring greening based on climate modeling and observational data, and showed that the increase in LAI in the Northern Hemisphere between 1982 and 2011 triggered an additional SM deficit, which could lead to additional high temperature associated with heat waves due to land-atmosphere feedbacks18,19,20,21. This apparent land cover change-induced modulation of the climate system and ecosystems has attracted the full attention of ecologists and climatologists.

However, previous studies have neglected the role of the saturated water vapor pressure deficit (VPD) in summer SM variability. VPD is expressed as the difference between the saturated water vapor pressure and the actual water vapor pressure. As a composite indicator that is closely related to temperature and relative humidity, VPD is one of the controlling variables of atmospheric evaporative demand22. The increased VPD is a key driver of increased aridity and drought under global warming11,23,24,25. During extreme heat wave events, increases in VPD will increase atmospheric evaporative demand, which will allow further enhancement of evapotranspiration and accelerate soil moisture loss26,27. In addition to this, VPD and SM are frequently tightly coupled due to land-atmosphere interactions, i.e. lower soil moisture is often accompanied by drier air28. Therefore, VPD not only has an important impact on the water cycle, but also plays a key role in modulating photosynthesis and carbon uptake and plant growth. Due to the effect of global warming, the asymmetry of the change rate of actual and saturated vapor pressure has led to VPD showing an increasing trend, although its increasing trend has been alleviated recently29. Compared to the greening of vegetation in spring, the role of increased VPD for summer SM change has not won sufficient attention. What’s more, the relative impacts of increased VPD and spring greening on summer SM changes remain elusive, and identifying these impacts is important for better understanding changes in the terrestrial water cycle and helping ecosystems adapt to climate change.

Here, we assess the relative impacts of VPD and vegetation change on summer SM in the Northern Hemisphere over the past four decades based on state-of-the-art global reanalysis data and multiple statistical methods. We analyzed change trends in Northern Hemisphere Spring (LAI), summer VPD and summer SM. To determine the effects of spring LAI and summer VPD to summer SM, which is a challenging work28,30, we used ridge regression, an experimental simulation approach based on linear least squares regression methods and an interpretable machine-learning model coupled with a random forest model to conduct the sensitivity analysis. On this basis, the relative contributions of VPD and LAI to summer SM were calculated. We aimed to address the following questions: (a) What changes has happened to summer soil moisture in the Northern Hemisphere under the recent climate warming? (b) How can we accurately quantify the contribution of VPD to summer SM variability? (c) What are the relative effects of VPD and spring vegetation greening on summer SM? We hope this study will be helpful in the development of climate change adaptation strategies and the enhancement of surface water management.

Results

Relationships between LAIMAM, VPDJJA and SMJJA

We firstly analyzed the change trends of LAIMAM, VPDJJA and SMJJA in the vegetative areas of the Northern Hemisphere in recent decades. The results plotted in Fig. 1 show that SMJJA shows a decreasing trend, while both LAIMAM and VPDJJA show an increasing trend. Specifically, the standardized anomalies of regional mean SMJJA in the Northern Hemisphere show a significant decreasing trend with a slope of -0.059 ± 0.011 per year, while the standardized anomalies of LAIMAM and VPDJJA show significant increasing trends of 0.063 ± 0.010 per year and 0.074 ± 0.008 per year, respectively. These results reveal the vegetation greening in spring and the atmosphere and soil drying in summer in the Northern Hemisphere. In addition, the spatial characteristics of LAIMAM, VPDJJA and SMJJA trends were analyzed. As for SMJJA, there are large differences between regions. The decreasing regions are mainly distributed in the western and southern parts of North America, south-central Europe, western Asia and eastern Asia, accounting for 69.3 percent of the total, of which the decreasing rate is greater than 0.025 m3·m−3·decade−1 in the south-western part of North America, western Asia and northern Mongolia. In the high latitude regions and parts of southwestern China, SMJJA shows an increasing trend. The proportion of regions showing an increasing trend in LAIMAM is 69.8%, mainly in the south-east of North America and in most of the Eurasian continent. The proportion of regions showing an increasing trend in VPDJJA is 87.0%, with the exception of the high latitudes of North America, and parts of Asia and Europe. In particular, the trend of VPD rise exceed 1.75 hPa·decade−1 in southwestern North America, southern Europe and western Asia. We also analyzed the trend of evapotranspiration and precipitation (Fig. S1), and found that evapotranspiration and precipitation showed increasing trends in most of the regions, but there were large differences in spatial distribution.

Fig. 1: Trends in LAIMAM, VPDJJA and SMJJA during 1982-2020.
figure 1

a Interannual standardized anomalies of the area-weighted average of LAIMAM (green dotted lines), VPDJJA (orange dotted lines) and SMJJA (blue dotted lines) for all northern region. The dashed lines indicate the least-squares linear regression for the variables corresponding to the colors. (b–d) denote the spatial distribution of LAIMAM, VPDJJA and SMJJA trends. The dotted area indicates that the trend is significant (P-value < 0.05). The white areas represent areas without vegetation cover that are masked off.

Full size image

The relationship between LAIMAM, VPDJJA and SMJJA presented in Fig. 1 shows that both vegetation greening in spring and atmospheric drying in summer significantly contributed to soil drying in summer (corr(SMJJA, LAIMAM) = -0.62, corr(SMJJA,VPDJJA) = -0.81). We further explored their relationships using partial correlation analysis, and the results showed large spatial differences (Fig. 2). We selected SMMAM, LAIMAM, PreJJA, TaJJA and VPDJJA as the variables affecting SM changes in summer (Methods). When conducting the partial correlation analysis between summer SM and LAIMAM (VPDJJA), the other main influence factors mentioned above were set as control variables. The areas where LAIMAM and SMJJA were negatively correlated were mainly concentrated at high latitudes, whereas VPDJJA and SMJJA were negatively correlated in almost all regions except for southwestern North America and parts of southern Asia. To detect the potential effect of human activities on the relationship between LAIMAM, VPDJJA and SMJJA, we explored the results under the agricultural subsurface. In cropland, correlations between VPD(LAI) and SM can be disturbed by measures such as irrigation and field management, thus differing from what happens in the natural environment. Meanwhile, we assume that the greater the share of cropland, the more human activities will affect the correlation. The identification of farmland types was performed based on GFSAD data (Methods), the spatial distribution of which is shown in Fig. S2. The results showed that the corr(SMJJA, LAIMAM) of rainfed farmland was lower than that of irrigated farmland under the same farmland occupancy condition, indicating that the role of spring greening on summer soil drying was higher in rainfed farmland. The corr(SMJJA, LAIMAM) gradually shifted to a positive correlation under both types of farmland as the proportion of farmland within the grid point increased. This indicates that the effect of greening in the spring on soil drying in the summer is weakening. On the contrary, corr(SMJJA, VPDJJA) remains at numerically higher negative values as the proportion of farmland increases. These findings indicate that the role of vegetation greening is more affected by human activities, while the role of VPD will not be significantly changed.

Fig. 2: Partial correlations between SMJJA and either LAIMAM or VPDJJA.
figure 2

a The spatial distribution of the partial correlation coefficients between LAIMAM and SMJJA. b Spatial distribution of partial correlation coefficients between VPDJJA and SMJJA. c Partial correlations affected by agricultural extent. The blue and red colors represent the partial correlation coefficients of LAIMAM and VPDJJA with SMJJA respectively. The diamonds and dots represent irrigated and rainfed farmland, respectively. The percentage of grid point farmland was identified based on the GFSAD1KCD data. X-axis shows 11 categories in 10% intervals from 0% to 100% (including 0). The numbers in different colors above X axis indicate the area proportions occupied by the corresponding category.

Full size image

Relative impacts of LAIMAM and VPDJJA on SMJJA

In order to obtain reliable results on the sensitivity of SMJJA to LAIMAM and VPDJJA, we analyzed them separately using three statistical methods. The statistical models developed by the three methods were evaluated, and the results showed that the R2 exceeded 0.7 in 73.6%, 97.5% and 70.6% of all the vegetated grids respectively, indicating that the fitting results were robust (Fig. S3). The results show that the sensitivities of SMJJA to LAIMAM and VPDJJA in the northern region are negative. This indicates that summer soil become drier as spring vegetation greens up and the atmosphere becomes drier in the summer. Specifically, the results for the sensitivity of SMJJA to LAIMAM were -0.15 ± 0.72, -0.14 ± 0.65, and -0.06 ± 1.02 (0.1 m3 · m−3) for the three methods. At high latitudes, the sensitivity of SMJJA to LAIMAM is strong negative. In addition, the results of the three methods are consistently negative in southwestern North America, southern Europe, central Asia, and northeastern China. The proportions of negative values for all regions are 64%, 65% and 53%. In terms of the magnitude of the sensitivity values, the results of methods 1 and 2 are similar and larger than those of method 3. The sensitivity results of SMJJA to VPDJJA are negative in most regions (99%, 96% and 93% of the total regions have negative values). The results for the three methods were -1.49 ± 1.09, -0.38 ± 0.37 and -0.98 ± 1.15 (101 m3·m−3·hPa−1). In terms of the magnitude of sensitivity values, the results of method 1 were large, method 3 was the next largest, and method 2 was the smallest.

We further explored the joint distribution of the sensitivity of SMJJA to LAIMAM and VPDJJA and present the results in Fig. S4. We did not show the case where the sensitivity results of SMJJA to VPDJJA is greater than 0 since it owns little area proportion in the results of all three methods. As shown in Fig. S4a–d, the proportion of area where the sensitivity of SMJJA to LAIMAM bellows 0 and the sensitivity of SMJJA to VPDJJA bellows 0 in the three methods is 63.5%, 62.3% and 32.4%. The average sensitivities of SMJJA to LAIMAM are −0.38 ± 0.68, −0.35 ± 0.62, −0.31 ± 0.95(0.1 m3·m−3) and average sensitivities of SMJJA to VPDJJA are −1.55 ± 1.06, −0.40 ± 0.37, −1.03 ± 1.13(101 m3·m−3·hPa−1). In this case, the summer soil will become drier as the spring vegetation becomes greener and the summer atmosphere drier. Conversely, in those regions that the sensitivity of SMJJA to LAIMAM bellows 0 and the sensitivity of SMJJA to VPDJJA exceeds 0, the average sensitivities of SMJJA to LAIMAM are 0.26 ± 0.62, 0.24 ± 0.57, 0.25 ± 0.87(0.1 m3·m−3) and average sensitivities of SMJJA to VPDJJA are −1.40 ± 1.12, −0.39 ± 0.38, −1.06 ± 1.18(101 m3·m−3·hPa−1) (Fig. 3).

Fig. 3: Spatial distribution of the sensitivity of SMJJA to LAIMAM and VPDJJA based on the three assessment methods.
figure 3

a represents the mean values of the sensitivities obtained by the three methods and the corresponding standard deviations. Different points correspond to the mean values of the results of the different methods, and the line gives the magnitude of the corresponding one standard deviation. (b–g) represent the spatial distribution of the sensitivities, respectively. The insets in the figure show the proportion of area occupied by positive sensitivities (blue) and negative sensitivities (brown).

Full size image

Based on the sensitivity assessment, the contributions of LAIMAM and VPDJJA to SMJJA were further explored separately. In general, the contributions of LAIMAM and VPDJJA to SMJJA in the northern region were both negative, i.e., both spring vegetation and summer atmospheric changes made summer soil moisture dryer. Specifically, for the LAIMAM -induced changes in SMJJA, the results of the three methods were -0.06 ± 0.83, -0.07 ± 0.70, and -0.07 ± 1.12 (10−3m3·m−3·decade−1). From the results of the spatial distribution of the three methods, the regions that showed negative values for the contributions were the southwestern part of North America, the southern part of Europe, the central part of Asia and the northern part of China, in addition to the high latitudes. For the change of SMJJA caused by VPDJJA, the results of the three methods are -2.64 ± 2.77, -0.56 ± 0.71 and -2.00 ± 3.46 (10−3m3·m−3·decade−1), and the magnitude of the results varies greatly. From the results of the spatial distribution of the three methods, except for the northwestern region of North America and some parts of the Asian and European continents, which have negative contributions, and for the southwestern part of North America, where the results of the three methods are inconsistent, the contribution values of all other regions show negative contributions, and the magnitude of their values is much larger than that of the contribution of LAIMAM.

It is also worth noting that when both VPDJJA and LAIMAM contribute negatively to SMJJA, this can exacerbate the SM deficit, which is more likely to trigger extreme dry events. Due to this, we further explored the joint distribution of the contributions of LAIMAM and VPDJJA to SMJJA and present the results in Fig. S5. We did not show the results where contributions of VPDJJA were positive due to the relatively small percentage of area. As shown in Fig. S5a–d, the proportion of area where both the contribution of LAIMAM and VPDJJA to SMJJA bellows 0 in the three methods is 47.5%, 45.4% and 27%. These areas are concentrated in southern North America, parts of Europe, and in most of Asia. The average contributions of LAIMAM to SMJJA are − 0.54 ± 0.64, −0.48 ± 0.52, −0.54 ± 1.09 (10−3m3·m−3·decade−1) and average contributions of SMJJA to VPDJJA are −3.39 ± 2.34, −0.73 ± 0.61, −2.55 ± 3.38 (10−3m3·m−3·decade−1). These areas are more likely to have summer SM deficits, as both VPDJJA and LAIMAM contribute negatively, which is noteworthy. Conversely, in those regions that the contributions of LAIMAM to SMJJA exceeds 0 and the contributions of VPDJJA to SMJJA bellows 0, the average contributions of LAIMAM to SMJJA are 0.50 ± 0.63, 0.45 ± 0.54, 0.40 ± 0.81 (10−3m3·m−3·decade−1) and average contributions of VPDJJA to SMJJA are − 3.15 ± 2.34, −0.76 ± 0.65, −2.54 ± 3.41 (10−3m3·m−3·decade−1). In these regions, the negative contribution of VPDJJA can be somewhat mitigated by the positive contribution of LAIMAM, weakening the SMJJA deficit (Fig. 4).

Fig. 4: Spatial distribution of LAIMAM and VPDJJA contributions to SMJJA using the three assessment methods.
figure 4

(a) shows the distribution of the contributions obtained by the three methods. The different points correspond to the mean of the results of the different methods and the line gives the magnitude of the corresponding one standard deviation. (b–g) denote the spatial distribution of the three methods and the two contributions, respectively. Positive contributions indicate that a change in the variable increases the SMJJA during the study period, while negative contributions indicate that a change in the variable decreases the SMJJA. The insets in the figure show the proportion of area occupied by positive (blue) and negative (red) contributions.

Full size image

The above results show that overall, both LAIMAM and VPDJJA exert negative contribution to SMJJA in the northern region. Spatially, the contribution of LAIMAM can be positive or negative, but the contribution of VPDJJA exceeds that of LAIMAM in many regions. Based on the above contribution results, we further calculated the ratio of the contribution of LAIMAM and VPDJJA to SMJJA. Although the results of the three methods show a large spatial inconsistency, it can be concluded by the area share of each ratio that the areas where the absolute value of the ratio of contribution is less than 1 account for the majority of the ratio. The results of method 1 show that 36% of the area with \({Con}\left({{LAI}}_{{MAM}}\right)/{Con}\left({{VPD}}_{{JJA}}\right)\) in the range of 0–25% is mainly distributed in some parts of North America, central and northeastern Asia, and 30% of the area with the range of -25%-0% is mainly distributed in the southeastern part of North America, northern Europe, and a small part of the east-central part of Asia. The results of method 2 show that \({Con}\left({{LAI}}_{{MAM}}\right)/{Con}\left({{VPD}}_{{JJA}}\right)\) between 0% and 25% accounted for 15% of the area, and between -25% and 0% accounted for 14% of the area, which were mainly found in the western and northeastern parts of Asia. The results of method 3 show that \({Con}\left({{LAI}}_{{MAM}}\right)/{Con}\left({{VPD}}_{{JJA}}\right)\) has 56% of the area between 0% and 25%, and its distribution area is more extensive and dispersed (Fig. 5).

Fig. 5: Magnitude and spatial distribution of the contribution of LAIMAM and VPDJJA to SMJJA using the three assessment methods.
figure 5

(a–c) represent the result of method 1-3, respectively.

Full size image

Comprehensive comparison of LAIMAM and VPDJJA contributions to SMJJA

The spatial distribution of \({Con}\left({{LAI}}_{{MAM}}\right)/{Con}\left({{VPD}}_{{JJA}}\right)\) in the northern region is analysed by combining the assessment results of the three methods, and is divided into four categories, namely, [1,+∞), (0,1), (-1,0) and [-∞, -1). Where (0,1) and (-1,0) are the categories of Con(VPDJJA) over Con(LAIMAM). Meanwhile, the agreements of three methods (i.e., the number of methods with the same classification result) are given for the result of each grid point. Since Con(VPDJJA) is basically negative, the sign of contribution ratio is opposite to that of Con(LAIMAM). Among them, it is worth focusing on the regions with contribution ratios of (0,1) and (-1,0) with agreements of 3. Specifically, regions with a contribution ratio of (0,1) and agreement level of 3 account for 7.5% of the area and are mainly located in parts of south-western North America, central-eastern Europe, and central and eastern Asia. The proportion of areas with a contribution ratio of (-1,0) and agreement of 3 is 5.9%, and these areas have a positive Con(LAIMAM), which are mainly located in parts of northern North America, northern Europe, northern Asia, southern Russia, and parts of southeastern Asia. Besides, the contribution ratio in [1,+∞) and (-∞, -1] are the regions where Con(VPDJJA) is less than Con(LAIMAM). Among these regions, the proportion of regions with agreement level of 3 is small (not shown in Fig. b), which suggests that there are large variations in different methods to obtain more robust results (Fig. 6).

Fig. 6: Comprehensive analysis of the results of LAIMAM and VPDJJA contribution ratio to SMJJA.
figure 6

Represented in (a) is the spatial distribution of patterns with different levels of agreements for the three methods. The inset in the lower right corner shows the legend of (a). Where agreement is represented by the number of methods with the same results as the ensembled results, and a number of 3 indicates that the results of the three statistical methods are consistent and that the results have high reliability. (b) Indicates the proportion of area occupied by the 12 subgroups presented in a. The size of each circle corresponds to the color in the legend in the lower right corner.

Full size image

Discussion

In this study, we find spatial differences in the sensitivities of SMJJA to LAIMAM and VPDJJA in the Northern Hemisphere, with the overall contributions of VPDJJA to soil moisture exceeding that of LAIMAM. Although changes in spring vegetation (more greening) accelerate soil moisture loss by increasing plant transpiration and thus evapotranspiration from the surface31,32; moisture entering the atmosphere affects the distribution of precipitation, generating spatial differences in response to atmospheric circulation and large-scale water vapor convergence and divergence. Therefore, the variation of vegetation through influencing P-ET and thus soil moisture varies among regions. In this paper, the changes in summer SM due to spring vegetation changes separated by three methods show positive values in high latitudes and negative values in middle and low latitudes, which is basically consistent with the results simulated by the Earth system model of the previous study31. What’s more, by separating the contributions of different factors through multiple statistical methods, we find that the magnitude of the contribution value of VPDJJA has exceeded the contribution of LAIMAM in many regions. Specifically, in parts of northern North America, northern Europe, northern Asia, southern Russia, and parts of southeastern Asia, the contribution of VPDJJA exceeds the contribution of LAIMAM and the contribution of LAI is positive, and the positive contribution of LAI partially offsets the negative contribution of VPD, resulting in a relatively small downward trend in SM. In contrast, in parts of southwestern North America, central-eastern Europe, and central and eastern Asia, the contribution of VPD exceeded the contribution of spring LAI and the contribution of LAI was negative, and the negative contribution of both was superimposed, resulting in a rapid decline in soil moisture in these regions, which corroborated the previous studies that reported the occurrence of drought more frequently in these regions17,33,34,35,36,37,38. VPD and soil moisture are closely linked through land-atmosphere interactions, i.e., higher VPD increases atmospheric evaporative demand, allowing more water to enter the atmosphere from the land via evapotranspiration and leading dry soils. In addition, VPD is also the main factor influencing canopy conductance of vegetation, so its influence on regional moisture changes can exceed that of vegetation. In this study, the contribution of summer VPD changes to summer soil moisture analyzed by three methods is negative in roughly 80% (area share) of the area and positive in the rest of the area, which is also due to the difference in sensitivity and trend of changes. In addition to this, anthropogenic impacts such as irrigation of agricultural land can also interfere with the impact of LAI. In summary, in addition to the need to pay attention to the feedback of summer moisture from the apparent change in surface cover, such as the greening of spring, the increase in atmospheric water vapor deficit may also have an important effect on summer soil moisture, which has not received sufficient attention39.

The effects of vegetation and atmosphere conditions on SM are complex, and it is challenging to isolate the separate contributions of the two to summer moisture. The two mainly affect soil moisture by altering changes in precipitation and evapotranspiration and thus soil moisture. In this study, temperature, precipitation, LAI, VPD, and spring soil moisture were selected as the main influences on summer soil moisture from the consideration of changes in precipitation and evapotranspiration. Radiation level was not selected due to its high covariance with temperature. The inclusion of spring soil moisture takes into account the fact that the memory of soil moisture in the earlier period can have an impact on soil moisture as well as P-ET changes in the later period (the so-called SM-P-ET feedback), which should be separated from the analyses in this paper. In the analysis of sensitivity, this study takes a total of three main approaches: multiple linear regression, least squares-based imputation experiments, and the combination of a random forest model and an interpretable machine learning model. The method of multiple linear regression is easy to understand and operate to obtain the sensitivity of the target variable to the explanatory variables; the attribution experiment based on the least squares method can intuitively obtain the sensitivity through the change of the target variable caused by the change of the explanatory variables; and the random forest is also widely used due to the advantages of the good handling of the nonlinear relationship and the low requirement of the distribution of the data. In analyzing the results, the results of the three methods were considered and method recognition was used to evaluate the results of the contribution analysis. Although the method of this indicator is simpler, it ensures the reliability of our results to a certain extent. In addition to this, the ERA5-Land dataset we chose is the reanalysis data that fuses multiple source datasets. Although the model is driven in such a way that the relationship between VPD and LAI and summer SM is determined by a set functional relationship, this concern should be ignored because of the better agreement between the reanalyzed data and the measured site data.

In the context of future global warming, VPD is projected to increase further on a global scale as temperatures rise. In contrast, projected soil moisture changes are heterogeneous and uncertain. Because the effects of VPD on ET are likely to be amplified at higher temperatures and the effects on plants themselves will be further deepened, more effort needs to be devoted to changes in the contribution of VPD to LAI under future scenarios. Although recent studies have demonstrated the critical role of VPD on plant physiological processes and ecosystem productivity, the effect of VPD on soil moisture seems to be overshadowed by the obvious change of vegetation greening in spring. Our results, which are reliable through multiple statistical methods, suggest that for the northern region the impact of VPD in summer will outweigh the effect of vegetation greening in spring. This will help to deepen the understanding of the regional water cycle and provide theoretical support for regional water management.

Methods

Datasets

The new generation of ERA5-land reanalysis data provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) was utilized in this study. The reanalysis data was reliable from the comprehensive analysis of ground observation data and remote sensing data with a spatial resolution of 0.1°. The data provided continuous global surface variables and had a higher resolution and reliable accuracy and therefore was widely applied40,41. We selected monthly averaged data from 1982 to 2020. The variables included top layer (0–7 cm) of soil moisture (Volumetric soil water layer 1), air temperature (2 m temperature, Ta), dewpoint temperature, surface net solar radiation, surface net thermal radiation and total precipitation (Pre).

The vegetation index (LAI) is an important parameter for studying global changes in vegetation structure and function. Here, LAI is used as an observational proxy for vegetation greenness. The LAI dataset was derived from the PKU GIMMS Normalized Difference Vegetation Index product (PKU GIMMS NDVI, version 1.2)42, which was generated based on biome-specific BPNN models that employed GIMMS NDVI3g product and 3.6 million high-quality global Landsat NDVI samples. This global satellite product has a spatial resolution of 8 × 8 km2 at biweekly intervals and has been widely applied in environmental science. We resampled it to 0.1° and calculated the spring LAI average for the period 1982–2020.

The Moderate Resolution Imaging Spectroradiometer (MODIS) land cover type data set (MCD12Q1) with the classification scheme of the International Geosphere-Biosphere Program was used to determine land use type for 2001–2019 with spatial resolution resampled to 0.01°. We also used the Global Food Security-support Analysis Data (GFSAD) at nominal 1 km and resampled it to 0.01° resolution. The Crop Dominance (CD) product (GFSAD1KCD)43 included in this dataset provided the map of irrigated and rainfed cropland, which can be helpful to analyze the effect of VPD and LAI on summer SM. We aggregated the MCD12Q1 data and GFSAD1KCD data into 0.01°. The crop fractions of each 0.1° × 0.1°grid according to the spatial resolution of climate data were calculated from the 10 × 10 (0.01°) grids from datasets mentioned above. The non-vegetated areas were identified by the MCD12Q1 data and masked in following results. The northern region in this article referred to the spatial range of 30°N-90°N, so the MAM and JJA months represented the Spring and Summer season.

Trend Analysis

The trends in LAI, VPD and SM were estimated using the linear least square regression method with a two-tailed t-test at the pixel scale with the following equation:

$${slope}=\frac{n\times \mathop{\sum }\nolimits_{i=1}^{n}i\times {a}_{i}-\left(\mathop{\sum }\nolimits_{i=1}^{n}i\right)\mathop{\sum }\nolimits_{i=1}^{n}{a}_{i}}{n\times \mathop{\sum }\nolimits_{i=1}^{n}{i}^{2}-{\left(\mathop{\sum }\nolimits_{i=1}^{n}i\right)}^{2}}$$
(1)

where slope represents a linear trend; ai represents the variable in year i; n is the number of time series.

Besides, the normalized anomalies (so-called Z-score normalization) of area-weighted mean LAI, VPD and SM were calculated to further estimate their trends. We estimated the trends in their normalized anomalies using the linear least square regression method with a two-tailed t-test. The normalized anomalies were calculated using the following equation:

$${X}_{i}^{{\prime} }=\left({X}_{i}-\bar{X}\right)/{\sigma }_{X}$$
(2)

where \({X}_{i}^{{\prime} }\) represents the normalized anomaly of the variable X at time i; \({X}_{i}\) represents the value of the variable X at time i; \(\bar{X}\) and σX represent the mean and standard deviation of the variable X over the whole study period, respectively.

Partial Correlation Analysis

To determine the magnitude of the relationship between spring LAI, summer VPD and summer SM in the Northern Hemisphere, while controlling for the effect of other variables, we conducted a partial correlation analysis. Partial correlation analysis is a statistical technique used to analyze the correlation between two variables while controlling for the effects of one or more other variables44.

Prior to conducting the partial correlation, we identified the main factors affecting summer soil moisture. Based on water balance theory, precipitation and evapotranspiration are considered to be the main factors in soil moisture changes, while the latter can be considered as a function on energy, vegetation cover and atmospheric evaporative demand. Also, the effect of pre-existing soil conditions was taken into account (memory of soil moisture). Therefore, we have selected SMMAM, LAIMAM, PreJJA, TaJJA and VPDJJA as the variables affecting SM changes in summer.

After that, we conducted partial correlation analysis in each individual grid cell. Specifically, we set the other main influence factors mentioned above as control variables when analyzing the correlation between summer SM and LAIMAM(VPDJJA). The significance of the partial correlations was evaluated at a threshold of p-value < 0.05.

Sensitivity Analysis

To isolate the relative impacts of LAIMAM and VPDJJA on SMJJA, we carried out analysis using three statistical methods to gain robust results. Sensitivity is used to measure the extent to which the dependent variable is disturbed by the independent variable.

Method 1: To avoid the effect of high correlations among the independent variables, ridge regression analysis was performed at individual grid points. In ridge regression, the introduction of a regularity term in the standard least squares objective function, controlled by a tuning parameter λ,helps stabilize the estimated coefficients. The ridge regression objective function can be expressed as follows:

$${\hat{\beta }}^{\wedge }=\mathop{\sum }\limits_{i=1}^{n}{\left({y}_{i}-{\beta }_{0}-\sum {\beta }_{i}{x}_{i}\right)}^{2}+\lambda \sum {\beta }_{i}^{2}$$
(3)

where \({\hat{\beta }}^{\wedge }\) represents the estimated regression coefficients, yi is the dependent variable, β0 is the intercept term, and βi represents the regression coefficient for the independent variable xi. The tuning parameter λ determines the degree of shrinkage applied to the coefficients. When λ is set to zero, the regularization term exerts no effect, and ridge regression reduces to ordinary least squares regression. The larger λ becomes, the more pronounced is the shrinkage applied to the coefficients, which can effectively mitigate the impact of multi-collinearity. In each grid cell, and the tuning parameter λ was determined based on the Variance Inflation Factor (VIF) of the independent variables. The following formula was used to calculate the VIF:

$${VI}{F}_{i}=\frac{1}{1-{R}_{i}^{2}}$$
(4)

where \({R}_{i}^{2}\) denotes the coefficient of determination between the ith independent variable and all other independent variables. Here, a VIF value less than 3 suggests an acceptable level of multi-collinearity45.

Method 2: Based on method1, two types of simulation experiments were carried out to obtain the impacts of the change of LAIMAM and VPDJJA46,47. The first type of simulation experiment (Sall) was a normal model run, and all the variables (i.e., SMMAM, LAIMAM, PreJJA, TaJJA, and VPDJJA) were set to change with time. The second type of simulation experiments (\({S}_{{{LAI}}_{{MAM}}}\) and \({S}_{{{VPD}}_{{JJA}}}\)) were the simulation experiments conducted by holding LAIMAM and VPDJJA constant, respectively. Then the differences between the simulation results of the first type (Sall) and the second type (\({S}_{{{LAI}}_{{MAM}}}\) and \({S}_{{{VPD}}_{{JJA}}}\)) to estimate the sensitivities of SMJJA to VPDJJA (βVPD)and LAIMAM (βLAI). Specifically, βVPD and βLAI were calculated using the following equations:

$$\Delta {y}_{({S}_{{all}}-{S}_{{{LAI}}_{{MAM}}})}={\beta }_{0}+{\beta }_{{LAI}}\times \Delta {{LAI}}_{({S}_{{all}}-{S}_{{{LAI}}_{{MAM}}})}+\varepsilon$$
(5)
$$\Delta {y}_{\left({S}_{{all}}-{S}_{{{VPD}}_{{JJA}}}\right)}={\beta }_{0}+{\beta }_{{VPD}}\times \Delta {{VPD}}_{\left({S}_{{all}}-{S}_{{{VPD}}_{{JJA}}}\right)}+\varepsilon$$
(6)

where ∆y, \(\Delta {{LAI}}_{({S}_{{all}}-{S}_{{{LAI}}_{{MAM}}})}\) and \(\Delta {{VPD}}_{({S}_{{all}}-{S}_{{{VPD}}_{{JJA}}})}\) represent the differences in SMJJA, LAIMAM and VPDJJA between the two types of experimental simulations, respectively; β0 is the intercept; ε is the residual. This approach fixes all other variables to their center or baseline values, so that any changes observed in the output will be clearly attributed to changes in individual variables thus calculating the sensitivity of the variables.

Method 3: We use explainable machine learning (SHapley Additive exPlanations) to study the effects of LAIMAM and VPDJJA on SMJJA. Random forests are one of the data-driven machine learning algorithms based on a bootstrap aggregating strategy for improving results stability, and it requires no statistical assumptions on predictors and target variables using sufficient numbers of data48. For this purpose, we first train Random forests models at each grid and then apply SHapley Additive exPlanations (SHAP) to isolate the marginal contributions of each predictor on the target variable. Specifically, we treat the SMJJA as the target variable and corresponding SMMAM, LAIMAM, PreJJA, TaJJA, and VPDJJA as predictors by a common hyperparameter setting (tree numbers = 100, leaf size =5). For one trained model, we apply SHAP dependence method to isolate marginal contributions of one predictor on the target variable49. We define overall SMJJA sensitivity as the slope estimated from Theil-sen regression according to SHAP dependence of LAIMAM and VPDJJA and by assuming that grid cell-level interactions are nearly linear. This method combines the advantages of bootstrap aggregating and non-distribution-assumption by random forest modeling, as well as advantages of global interpretations being consistent with the local explanations in the SHAP algorithm49,50,51, hence strengthening the robustness of the results

Contribution Analysis

Based on the sensitivity results from three methods, the relative contributions of LAIMAM and VPDJJA were calculated using the following equation47:

$${Con}={\beta }_{{Trend}}\times {\beta }_{{Sen}}$$
(7)

where Con represents the contributions of LAIMAM and VPDJJA to SMJJA, βTrend and βSen represent the trends and sensitivities of LAIMAM and VPDJJA.

After that, the ratio of the changes of SMJJA induced by LAIMAM and VPDJJA were calculated respectively:

$${Con}{Ratio}={Con}\left({{LAI}}_{{MAM}}\right)/{Con}\left({{VPD}}_{{JJA}}\right)$$
(8)

Where, Con(LAIMAM) represents the contribution of spring vegetation change to SMJJA, Con(VPDJJA) represents the contribution of summer VPD to SMJJA. The magnitude of the two contributions was determined by Con Ratio.