Introduction

As the main renewable freshwater resource, streamflow is important for sustaining human society, the ecological environment, and agricultural production, especially in the dry season when hydrological droughts can easily occur1,2,3. Although there is a strong agreement that global mean streamflow will increase in the future under global warming in response to increased precipitation, the changing patterns are likely to be highly heterogeneous4. Therefore, continued efforts to understand future streamflow changes at the regional scale are essential for developing local adaptive and sustainable water resource management strategies5.

The responses of regional water cycles to global warming extend beyond simple atmospheric adjustments and also include land–atmosphere interactions, such as the soil moisture–atmosphere feedback (SAF)6,7,8. For example, decreasing soil moisture weakens the cooling strength of evapotranspiration (ET) and increases the near-surface temperature, which subsequently enhances ET demand, dries the soil, and reduces runoff and streamflow9,10,11,12. Such a soil moisture-temperature feedback is generally positive in transition zones between wet and dry climates8. Additionally, soil moisture anomalies can modify precipitation by influencing the local ET supply and atmospheric circulation, which in turn affect streamflow. However, the interactions between soil moisture and precipitation are highly complex and both positive and negative relationships have been observed in dry soils13,14.

Since the establishment of the Global Land–Atmosphere Coupling Experiment-Coupled Model Intercomparison Project Phase 5, several efforts have been made to clarify the effects of SAFs on future water cycle changes15,16,17. In particular, SAF was projected to intensify the warming speed and extreme temperatures in extratropical regions (e.g., Europe, North America, and Australia), because decreased soil moisture over these regions in a warming future would decrease evaporative cooling, increase sensible heat, and warm the atmosphere18,19. The influence of SAF on precipitation changes is complicated owing to the presence of both direct and indirect effects, and the discrepancy among the climate models is much larger than that for temperature and ET14,15,20,21,22. A previous work suggested that the effect of SAF on precipitation varies by season and tends to be more pronounced during winter14. Recent studies have focused on precipitation minus ET (P-E) and found that SAF consistently increased P-E in drylands (e.g., Central North America, the Middle East, and North and South Africa) and diminished the seasonality of P-E in subtropical regions (e.g., Australia, the Mediterranean, and North America) and the Amazon, which would potentially mitigate declining water availability (e.g., streamflow) in the future16,17. However, because streamflow projections are not directly provided by climate models as they additionally require the use of hydrological models4, the influence of SAF on future streamflow changes has yet to be quantified. Moreover, unlike runoff and P-E, the streamflow is an integrated variable for a specific catchment. Previous studies have shown that although gridded runoff changes are not robust in some regions (e.g., East Asia and Central Africa), streamflow shows marked changes4,23. This indicates that climate models may have higher consistency at catchment scales, and it is necessary to investigate the influence of SAF on future streamflow changes that might be different from runoff responses at the grid scale.

Using meteorological forcing data from the Land Feedback Model Intercomparison Project with prescribed Land Conditions (LFMIP-pdLC) experiment in the Coupled Model Intercomparison Project Phase 6 (CMIP6) to drive an advanced land surface model named Conjunctive Surface–Subsurface Process version 2 (CSSPv2)24 (Supplementary Fig. 1), here we performed a series of hydrological simulations over China with the aims of: (1) investigating the connections between future streamflow changes and SAF, and (2) revealing the corresponding mechanisms. Our results show that the SAF significantly contributes to about 50% of the future streamflow changes in spring season with a spatially varying pattern (e.g., increasing streamflow in mid-latitude China but decreasing it in other regions). The underlying mechanisms are mainly attributed to soil moisture-precipitation interactions through thermodynamic and thermal atmospheric adjustments to SAF-induced warming.

Results

Reliability of CMIP6 and CMIP6/CSSPv2 simulations

Ensemble mean surface air temperature, precipitation, and evapotranspiration from the five CMIP6 models show a similar climatology to the state-of-the-art reanalysis products (ERA5&MERRA2), and that root mean square errors are within 3.5 K, 60 mm, and 20 mm at monthly time scales in most regions (Supplementary Fig. 2). Individual CMIP6 members also capture the monthly variations of temperature, precipitation, and evapotranspiration well, with correlation coefficients larger than 0.5 (p < 0.1) over most grids (Supplementary Figs. 3–5). The CMCC-ESM2 and EC-Earth show better performance than the other three models, partly due to higher model spatial resolutions25. In addition, these models also have good performance in modeling soil moisture dynamics26. The empirical statistical method16,17 (Supplementary Method) was used to quantify the strength of soil moisture feedbacks on precipitation and evapotranspiration. The soil moisture feedback on precipitation is highly heterogeneous at regional scales, and reanalysis datasets show differences over northern and western China (Supplementary Fig. 6). In general, the five CMIP6 models show similar patterns to ERA5, with positive and negative patterns over southern and northern China, respectively. And the spatial correlation is 0.3–0.5 (p < 0.1) between different CMIP6 models and the ERA5 reanalysis. The soil moisture feedback on evapotranspiration is generally positive and shows consistent spatial patterns across CMIP6 models and reanalysis datasets.

The land surface model simulation was evaluated based on the Kling–Gupta efficiency (KGE), an integrated metric that considers the correlation, variation, and bias of the simulation27. Figure 1b shows the KGE between modeled and observed monthly streamflow during 1961–1990 when human intervention was less prevalent. The KGE values are larger than 0.5 for both the OBS/CSSPv2 and CMIP6/CSSPv2 simulations over 25 catchments, with the median KGE values being 0.82 and 0.71, respectively. The relative simulation error of CMIP6/CSSPv2 is smaller than 10% for 85% of the basins. Although CMIP6/CSSPv2 shows smaller KGE values than the OBS/CSSPv2 simulation which may be related to the modeling uncertainties of the annual dynamics of meteorological forcings25, it is considered satisfactory as the KGE value exceeds 0.5 and the relative error is within 25%28. Moreover, CMIP6/CSSPv2 captures the seasonal cycle of monthly streamflow during 1961–1990, especially during the dry season (Fig. 1c). Compared with that of the individual members of the CMIP6/CSSPv2 simulation, the ensemble mean increases the KGE values by 2–8%. Therefore, the CMIP6 models reasonably represent the hydrological cycles and soil moisture-atmosphere feedbacks over China, and CMIP6/CSSPv2 shows high modeling accuracy in monthly streamflow.

Fig. 1: Station locations and evaluations of simulated monthly streamflow across China.
figure 1

a Streamflow stations and their record lengths during 1961–1990 in different catchments. b Boxplot of the Kling–Gupta efficiency (KGE) values between simulated and observed monthly streamflow during 1961–1990. OBS/CSSPv2 indicates that the meteorological forcings are from China gridded meteorological forcings combined with European Centre for Medium-Range Weather Forecasts Reanalysis v5 (CN05 + ERA5) datasets, while CMIP6/CSSPv2 indicates that the ensemble mean of five sets of streamflow simulation was driven by meteorological forcings from different Coupled Model Intercomparison Project Phase 6 (CMIP6) models. c Observed and simulated annual cycle of monthly streamflow (m3/s) averaged during 1961–1990. Blue shadings represent the range of all individual CMIP6 members.

Full size image

Future streamflow changes and contributions from SAF

CMIP6/CSSPv2 projection indicates a widely increasing trends of seasonal streamflow in 14 large river basins that contribute more than 90% of the total streamflow in China (Supplementary Fig. 7a). The SAF is found to significantly influence future streamflow changes in more than half of the basins in the spring season, including those of the Tarim, Haihe, Yellow, Pearl, Brahmaputra, Nu, and Lancang rivers and the upper and middle reaches of the Yangtze River (Supplementary Fig. 7b). In other seasons, however, a high discrepancy occurs among different simulations, and a significant effect of SAF occurs only in 2–3 basins (Supplementary Fig. 7b). Compared to the streamflow, SAF’s impacts on grid-scale seasonal runoff are insignificant throughout the seasons in most regions according to both CMIP6/CSSPv2 simulations and original CMIP6 outputs (Supplementary Fig. 8). This suggests the necessity to move from runoff to streamflow so as to filter out the grid-scale noise (due to different model resolutions and parameterizations) and enhance the detectability of SAF’s impacts on terrestrial water cycles. In addition, although precipitation minus evapotranspiration (P-ET) can represent water availability in land, it is not the same as runoff/streamflow because the terrestrial water storage (including water stored in snow, soil layers and groundwater) is changing under climate change29,30. Therefore, seasonal changes of P-ET show larger variations than the runoff (Supplementary Fig. 9), and the SAF’s influence on P-ET differs from that on runoff both in magnitude and spatial patterns (Supplementary Fig. 8). Such a difference further emphasizes the direct focus on streamflow to understand the impacts of SAF on freshwater availability instead of relying solely on P-ET.

Figure 2 shows the mean relative spring streamflow changes during 2080–2099 compared with those during 1980–1999, with the spatial patterns shown in Supplementary Fig. 10. In general, spring streamflow is projected to increase by approximately 20–130% in most basins under future global warming, except in that of the Pearl River, which shows a slight (−4%) decreasing trend (blue bars in Fig. 2). These results are consistent with those of previous studies that have projected future streamflow in various basins in China31,32,33,34,35. Basins located in northern China (e.g., Songhua, Liaohe, Haihe, and Yellow rivers) and the Tibetan Plateau (e.g., Brahmaputra and Nu rivers) show a higher increasing trend than the other basins, which may be related to the climatology of their smaller streamflows (Supplementary Fig. 10a) and colder climates that are more sensitive to global warming36. The insignificant streamflow changes in the Yangtze river basin are associated with decreasing spring streamflow over the southern part of its middle and lower reaches, whereas the decreasing trend in the Pearl river basin is generally found throughout the basin (Supplementary Fig. 10b). Compared with the results of CMIP6/CSSPv2, the increasing trend of streamflow is smaller in the data from the CMIP6-pdLC/CSSPv2 experiment, suggesting favorable effects of SAF on the increased spring streamflow (Supplementary Fig. 10c). As shown in Fig. 2, approximately 28–58% of the increased spring streamflow in mid-latitude China can be attributed to SAF. However, SAF tends to have a negative effect in the northwestern (e.g., Tarim Basin), northeastern (e.g., Liaohe Basin), and southern (Pearl Basin) parts of China, causing a 12–48% decrease in spring streamflow (Fig. 2 and Supplementary Fig. 10d). For the Yangtze and Pearl river basins, SAF and the other factors exert opposite influence on streamflow, making the total changes (blue bars) smaller than the SAF-forced changes. For example, SAF causes strong precipitation decline over the Pearl river basin (Supplementary Fig. 11c) while the other factors induce an increasing precipitation (Supplementary Fig. 11b). In general, most CMIP6 models (≥80%) are consistent in showing a positive effect of SAF on future spring streamflow changes over the mid-latitude China, but a negative effect over northern and southern China.

Fig. 2: Multi-model mean results of relative changes (%) of spring streamflow over different basins during 2080–2099 under different scenarios.
figure 2

The reference is the monthly streamflow climatology during 1980–1999. Different colors on the map designate different river basins, with red stars denoting river outlets. The blue bars show the results of the CMIP6/CSSPv2 experiment, indicating the streamflow changes induced by all factors. Red bars show the effects of soil moisture–atmosphere feedback (SAF), assessed using CMIP6/CSSPv2 minus CMIP6-pdLC/CSSPv2 experiments. Yellow bars show only the effect of SAF on precipitation. Hatched bars denote significant changes (p < 0.05). Markers show the results for the individual members.

Full size image

Although SAF can alter numerous near-surface meteorological forcings (e.g., temperature, shortwave radiation, and wind speed), thereby modifying streamflow, the results show slight differences when only the influence of SAF on precipitation was considered (Fig. 2). In addition, the spatial patterns of precipitation changes induced by SAF are highly consistent with those of the streamflow (Supplementary Fig. 11). For the alpine mountainous regions (e.g., the Tibetan Plateau), while changes of snowpacks (e.g., the time of snow melting) and the snow/rainfall partition are critical for spring streamflow37, the SAF causes slight changes of snow/rainfall and snow depth in winter (Supplementary Fig. 12), thus exerting minor impacts on the delay of winter snowpacks melting. In spring, SAF causes an increase of snow/rainfall, and the increased magnitude of rainfall is more than twice of snowfall. The warming effects of SAF, however, favor the snow melting and lead to little changes in spring snow depth. Therefore, the increased snowfall induced by SAF will melt rather than accumulate and, together with the increased rainfall (Supplementary Fig. 12), will contribute to the significant increase in surface runoff generation, which ultimately increases streamflow. This indicates that the influence of SAF on spring streamflow changes is mainly attributed to its effects on precipitation rather than other meteorological forcings and snow-related processes.

Mechanisms for the impact of SAF on future streamflow changes

The low-level moisture advection indicated by the water vapor flux at 925 hPa (WVF925), and mid-tropospheric vertical wind indicated by the vertical velocity at 500 hPa (ω500), are used to represent the moisture transport and precipitation processes, respectively38. The climatology of WVF925 and ω500 during 1980–1999 show that water vapor is mainly transported by the southerly winds from the South China Sea, and separately, an upward motion occurred over the regions south of 40° N (Fig. 3a). However, SAF is projected to accelerate land surface warming in the future15,18, which would reduce the surface pressure over land through heating, enhance the surface pressure gradient between the land and ocean, and increase the strength of southerly winds (Fig. 3b). As a result, more water vapor is transported to mid-latitude China (region I in Fig. 3c) and the cyclonic circulation wind anomaly over mid-latitude China forces an upward-moving anomaly (Fig. 3c). At the same time, a downward-moving anomaly occurs over southern and northern China (region II in Fig. 3), which depresses the formation of precipitation.

Fig. 3: Effects of soil moisture–atmosphere feedback (SAF) on circulation.
figure 3

a Mean horizontal water vapor flux at 925 hPa (WVF925) and vertical speed at 500 hPa (ω500) during 1980–1999 in the CMIP6 historical experiment. b Influence of SAF on changes in surface pressure (Ps) and horizontal wind vector at 925 hPa (V925) during 2080–2099 according to the CMIP6 high-end forcing scenario (SSP585) experiment. c Influence of SAF on the changes in WVF925 and ω500 during 2080–2099. Black boxes show the region II that includes south China (17–27°N, 108–122°E), northeastern China (45–55°N, 115–137°E), and northwestern China (36–50°N, 71–92°E). The other region in China is masked as region I.

Full size image

Figure 4 shows the precipitation changes induced by SAF from the perspective of the atmospheric water vapor budget. The precipitation changes are generally balanced by the changes in ET and meridional, zonal, and vertical water budgets. Furthermore, increased precipitation in mid-latitude China is mainly caused by changes in the vertical water vapor budget that can be further attributed to the significantly increased dynamic and nonlinear vertical terms (WB_D and WB_N in Fig. 4a). The small magnitude of the horizontal water budget is due to the contrary effects of horizontal dynamic and thermodynamic terms. However, the meridional and vertical water vapor budget terms are both responsible for the decreased precipitation in southern and northern China (Fig. 4b). In contrast to the vertical dynamic processes shown in Fig. 3, meridional thermodynamics are related to changes in the specific humidity gradient. Specifically, SAF has a stronger depression effect on the increase of specific humidity over the South China Sea than that over southern China, and it decreases the meridional specific humidity gradient (\(-\partial q/\partial y\, < \,0\)) between sea and land (Supplementary Fig. 13), which reduces the meridional water vapor budget over southern China (\(-\bar{\nu }\cdot \partial q/\partial y\, < \,0\)).

Fig. 4: Effects of soil moisture–atmosphere feedback on water cycles.
figure 4

Changes in precipitation (P), evapotranspiration (ET), and column integrated moisture budgets during 2080–2099 and 1980–1999 in a China, b mid-latitude China, and c southern and northern China. Moisture budget at the zonal direction is represented by UB and is further decomposed into the thermal (UB_T), dynamic (UB_D), and nonlinear terms (UB_N). Moisture budgets at meridional and vertical directions are similar, but are represented by VB and WB, respectively. Colored bars are ensemble means while error bars represent one standard deviation among different model members.

Full size image

Discussion

Compared to the previous works that focused on P-ET and found significant effect SAF over subtropical regions16,17, this study goes a step further by moving from grid-scale P-ET to spatially integrated variables (e.g., streamflow). The results reveal that future soil moisture changes can exert significant impacts on streamflow in specific seasons in regions that are not located in areas commonly considered SAF hotspots (e.g., southern and mid-latitude China). Moreover, the atmospheric thermal adjustment is also revealed to be important, except for the widely recognized dynamical adjustment. This suggests that model developers and scientists should pay more attentions to the soil moisture-atmosphere coupling processes when projecting and understanding the future streamflow changes.

The variations among CMIP6 models (markers in Fig. 2 and error bars in Fig. 4) are related to their different climate sensitivities (e.g., equilibrium climate sensitivity; ECS) that varying from 2.6 K (MIROC6) to 4.56 K (IPSL-CM6A-LR)39. However, current best estimate of the ECS still has large uncertainty (5–95% range is 2.3–4.7 K)40, underscoring the complexity of the climate system and the danger of trusting one CMIP6 in isolation4. This study uses multi-model ensemble means to reduce uncertainties among different models following previous works. Meanwhile, different CMIP6 models show similar signs of SAF’s impacts on future streamflow and results pass the significance test, suggesting the conclusion is robust and reliable. However, there is still a higher discrepancy among CMIP6 models during other seasons, which might cause an insignificant influence of SAF on streamflow (Supplementary Fig. 14). Future studies should take steps to reduce the uncertainties in representing the soil moisture–atmosphere coupling processes among the CMIP6 models41,42.

In addition, this study focuses on the SAF’s effects on natural streamflow without considering human impoundment systems (e.g., reservoirs), as the natural streamflow is a critically meaningful variable to forecast since it directly relates to upstream inflow into a reservoir system. Recent work has shown that the human impoundment systems, especially over the eastern China, will cause complex influences on future streamflow changes related to their different operation rules and can also modulate local climates43,44,45,46. However, although there are numerous reservoirs parameterizations, they are still not included in either CMIP6 models or CSSPv2 land models for national to continental-scale simulations because the highly heterogeneous reservoir operation rules need extensive parameter calibrations with considerable quantity of data that are usually unavailable47,48. Continued efforts are needed to better understand the interactions of human water regulation and SAF in determining spatial patterns of future fresh water resource changes.

Finally, considering that the LFMIP-pdLC experiment was performed at global scale, the current influence of SAF includes both local and non-local effects18. Therefore, future modeling studies using regional climate models are required to quantify the independent and joint contributions of local and non-local SAF to streamflow changes.

Methods

Data

Monthly streamflow observations during 1961–1990 at 25 stations throughout China were collected from local water resource departments and the global runoff data center49,50. Streamflow stations cover a wide range of climates and have long-term (>10 years) records (Fig. 1a). Daily temperature and precipitation from the 0.25° China gridded meteorological forcings (CN05) during 1961–201451, together with the 0.25° hourly surface air temperature, humidity, pressure, shortwave and longwave radiation, and 10 m wind speed from the European Centre for Medium-Range Weather Forecasts Reanalysis v5 (ERA5)52, were used as meteorological observations. Monthly precipitation, evapotranspiration and top 1 m soil moisture from ERA5 and the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA2)53 were used to evaluate CMIP6 models.

Daily meteorological forcing data from the CMIP6 historical experiment that considers all the external forcings (ALL) during 1961–2014 and the high-end forcing scenario (SSP585) during 2015–2099 were collected. The LFMIP-pdLC experiment shared the same settings as ALL and SSP585 during 1980–2014 and 2015–2099, respectively, except that the soil moisture was prescribed to a climatological state derived from the ALL experiment during 1980–2014 (Supplementary Fig. 1). Five CMIP6 models (CMCC-ESM2, EC-Earth3, IPSL-CM6A-LA, MIROC6, and MPI-ESM1-2-LR) were selected because they provided daily outputs in all three experiments. The Cumulative Distribution Function (CDF) quantile-based mapping method was used to correct the CMIP6 ALL simulation, which ensured that the simulated meteorological forcing data had the same CDF as the observations (CN05 and ERA5). The Equidistant CDF matching method, which directly considers distribution changes between future and historical periods54, was used to correct the SSP585 and LFMIP-pdLC experiments. Detailed information on the model is provided in Supplementary Table 1.

Experimental design

The CSSPv2 land surface model was used to simulate the streamflow. This model was well calibrated for the conditions in China and was shown to perform well in simulating terrestrial water cycles (e.g., snow, ET, runoff, and soil moisture)50,55.

Observed meteorological forcings were used to force CSSPv2 model (OBS/CSSPv2) to investigate the performance of the CSSPv2 model in modeling streamflow in China. Bias-corrected meteorological forcing data from the CMIP6 ALL and SSP585 experiments were used to drive the CSSPv2 model for 1961–2014 and 2015–2099, respectively, and the simulation was denoted as CMIP6/CSSPv2. Similarly, meteorological data from LFMIP-pdLC experiment was used to drive the CSSPv2 model for 1980–2099, and the simulation was named CMIP6-pdLC/CSSPv2, with the initial conditions derived from the last day of 1979 in the CMIP6/CSSPv2 simulation (Fig. S1). Each CMIP6/CSSPv2 and CMIP6-pdLC/CSSPv2 simulation had five members due to the use of different CMIP6 models, and future streamflow changes indicated differences between the ensemble mean streamflow climatology for 2080–2099 and 1980–1999. The changes were considered significant (p < 0.05) when at least four members showed the same positive (or negative) sign as the ensemble mean according to the binomial distribution17. The differences between the CMIP6/CSSPv2 and CMIP6-pdLC/CSSPv2 experiments represent the influence of SAF18. In addition, we performed another simulation that was the same as CMIP6/CSSPv2 except that it used the precipitation data from the CMIP6-pdLC/CSSPv2 simulation to investigate the relative importance of precipitation and other forcings induced by SAF. All the simulations were performed at a spatial resolution of 0.25°.

Atmospheric moisture budget

Based on the atmospheric moisture budget (Eq. (1)), changes in precipitation (\({P}^{{\prime} }\)) were approximately equal to changes in ET (\({E}^{{\prime} }\)) and horizontal (\(- < \,{V}_{h}\cdot \nabla q\,{ > }^{{\prime} }\)) and vertical (\(- < \,\omega \cdot {\partial }_{p}q\,{ > }^{{\prime} }\)) water vapor convergence:

$${P}^{{\prime} }\approx \;{E}^{{\prime} }+\frac{1}{{\rho }_{w}g}\left[- < \,{V}_{h}\cdot \nabla q\,{ > }^{{\prime} }- < \,\omega \cdot {\partial }_{p}q\,{ > }^{{\prime} }\right]\\ \approx \;{E}^{{\prime} }+\frac{1}{{\rho }_{w}g}\left[\underbrace{- { < \,\overline{{V}_{h}}\cdot \nabla {q}^{{\prime} }\, > }}_{{{\mbox{horizontal}}}\atop {\mbox{thermodynamic term}}}\;\;\underbrace{- { < \,{{V}^{{\prime} }}_{h}\cdot \nabla \bar{q} \, > }}_{{\mbox{horizontal}}\atop {\mbox{dynamic term}}} \;\;\;\;\underbrace{- { < \, {{V}^{{\prime} }}_{h}\cdot \nabla {q}^{{\prime} } \, > }}_{{\mbox{horizontal}}\atop {\mbox{nonlinear term}}}\right. \\ \qquad\qquad \left. \underbrace{- { < \,\bar{\omega }\cdot {\partial }_{p}{q}^{{\prime} } \, > }}_{{\mbox{vertical}}\atop {\mbox{thermodynamic term}}}\;\;\underbrace{- { < \,{\omega }^{{\prime} }\cdot {\partial }_{p}\bar{q} \, > }}_{{\mbox{vertical}}\atop {\mbox{dynamic term}}}\;\;\; \underbrace{- { < \,{\omega }^{{\prime} }\cdot {\partial }_{p}{q}^{{\prime} } \, > }}_{{\mbox{vertical}}\atop {\mbox{nonlinear term}}}\right]$$
(1)

where \({\rho }_{w}\) is the air density, \(g\) is gravity acceleration, “\( < \, > \)” means the vertical interpolation from land surface to the top of atmosphere, and \({V}_{h}\), \(\omega\), and \(q\) are horizontal vector wind, vertical wind, and specific humidity, respectively. The overbar means climatology. The water vapor convergence can be further separated into thermodynamic, thermal, and nonlinear terms56,57. Here we used monthly outputs at pressure levels to diagnose the mechanism of precipitation changes.