Ecological Response to Climate Change Across China From Combined Soil Temperature and Moisture Changes

A coupled soil temperature (ST) and moisture (SM) balance reflects a synthetic climate regime, having huge ecological impacts. This paper used ST and SM data from the European Center for Medium‐Range Weather Forecasts climate reanalysis‐Land and the Coupled Model Intercomparison Project Phase 6 and leaf area index (LAI) data from the Global Land Surface Satellite Product Suite. The focus was on understanding joint ST‐SM changes and the resulting ecological response across China. The results show that during 2000–2020, 24.5% of the land area in China experienced a warming‐drying trend resulting in a 9.7% LAI decrease, while 6.4% of the area experienced a warming‐wetting trend leading to an 8.6% LAI increase. During 2015–2100, 30.6% of the land area in China will be warmer and drier, while 55.2% of the area will be warmer but wetter across three shared socioeconomic pathways (SSP126, 245, and 585). Superimposed on the long‐term trends, there are also significant spatiotemporal variabilities in ST and SM on annual to decadal timescales. The LAI also showed substantial short‐term fluctuations in both typical regions and ecosystems despite consistent long‐term increases. Our findings suggest that ecosystems could be impaired on annual to decadal scales by adverse soil conditions in the twenty‐first century, but in terms of long‐term trends, ecosystems may be resilient partly because of the compensating effects of global warming and regional hydrological changes. Impact studies should thus focus more on annual to decadal soil‐ecosystem anomalous events.

land-atmosphere coupling and thus affect regional climates and ecosystems. Global warming is pumping increasingly more heat into Earth's surface and is modifying the terrestrial water cycle (Wu et al., 2013), which in turn is changing soil properties. Importantly, such changes are spatially inhomogeneous such as the so-called "wet gets wetter and dry gets drier paradigm" (Chou et al., 2009;Feng & Zhang, 2015). Potentially significant changes are expected across transition zones such as China (IPCC, 2021). As global warming continues to intensify, how the soil-plant-atmosphere continuum (SPAC) responds is of paramount importance to terrestrial habitats and human societies. It is therefore important to quantify the changes in soil temperature-moisture regimes, for example, warming-drying, warming-wetting, cooling-drying, and cooling-wetting, for the assessment of regional climate change and the preparation for adverse impacts on ecosystems and other terrestrial elements.
Because of the importance of SM to climate, ecosystems, and so on, there have long been studies of SM changes in respective disciplines (Denmead & Shaw, 1962;Philip, 1957). However, the scarcity of in situ measurements prevents the demonstration of large-scale SM-climate feedbacks (Dorigo et al., 2011;Robock et al., 2000). It was not until the 1970s that the development of numerical models created a new opportunity for the study of SM-climate interactions (Charney et al., 1977;Yeh et al., 1984). The early numerical sensitivity experiment showed a linkage between SM and precipitation on global and regional scales (Yeh et al., 1984). Importantly, Koster et al. (2004) conducted a multimodel experiment and proved that the hot spots with a strong coupling between SM and precipitation were mainly located in the transition zones between wet and dry climates. North China is among these hot spots, although the coupling is less intense than that in North America, Sahel, equatorial Africa, and India. From the perspective of physical processes, the long memory of SM links the spring variability to summer anomalies of precipitation in eastern China (Liu et al., 2017). In North China, the prolonged memory of SM favors the propagation of precipitation deficits into the soil and amplifies the impacts of droughts on the terrestrial environment .
Another role played by soil is associated with the heat cycle, which is interactive with the water cycle through evapotranspiration, albedo, and geochemical cycles in the climate system. Especially at high latitudes and altitudes, heat conditions in the form of ST, snow variations, and thawing and freezing cycles are closely linked to SM dynamics (Zhang et al., 2005). For the interdependence on land processes, on the one hand, the mean annual ST shows up to a 10 K difference from the corresponding air temperature, with marked variations in space and time (Lembrechts et al., 2022). On the other hand, ST can largely increase surface air temperature variability and persistence (Zhang et al., 2005). Meanwhile, a modeling study indicated that increased ST is favorable to convective development and precipitation (Fan, 2009;Xue et al., 2012).
More importantly, the combined effects of ST and SM changes may lead to intensified impacts of climate change (Grillakis, 2019). For example, warming soil can amplify drought at midlatitudes; in turn, drought can largely alter the effects of soil warming on terrestrial water cycles (Samaniego et al., 2018). With warming, soil droughts in China have intensified from the perspectives of drought severity, duration, and frequency during the past several decades (Wang et al., 2010). In addition, such joint changes in ST and SM can aggravate the impacts of regional climate change on microbial activity, soil respiration, and other environmental elements (Schindlbacher et al., 2012). Furthermore, the joint changes in ST and SM appear markedly inconsistent in space and time due to vegetation processes (Berg & Sheffield, 2018).
For the interdependency of water transfer in the land-atmosphere system, Philip (1966) proposed the concept of the SPAC, which highlights the influence of vegetation on ST and SM changes. On the global scale, extensive studies have found that more than half of moisture is lost from the soil into the atmosphere by vegetation (Lian et al., 2018;Zhang et al., 2016). Such a sizable proportion of SM loss is supplied by deeper soils where moisture is transferred by plant roots and thus occurs more slowly and profoundly affects regional SM-climate feedbacks. Across China, the role of vegetation in controlling water transfer in the SPAC increases with climate warming as a result of the extended growing seasons and enhanced water use efficiency . Furthermore, during past decades in China, satellite monitoring has detected significant vegetation changes following climate change along with other direct and indirect factors (e.g., land-use, carbon dioxide fertilization, and nitrogen deposition) (Chen et al., 2019;Piao et al., 2020). For example, on the Tibetan Plateau (TP), significant greening trends caused by climate factors have affected over 7.63% of the total territory (Zhong et al., 2019). In the Yellow River Basin, vegetation greening and browning trends have expanded more than 94% of the areas of the region, accompanied by a warming climate, afforestation, and increasing anthropogenic activities (Tian et al., 2021). The region of Inner Mongolia has experienced increasing vegetation coverage but degradation of meadow vegetation with 10.1029/2022EA002640 3 of 24 intensified temperature and precipitation variations . However, how vegetation changes across China are associated with the coupled changes in ST and SM has yet to be adequately investigated.
From a holistic perspective of soil temperature-moisture-vegetation and not considering only the singular changes in these elements, these coupled changes in the soil-vegetation-atmosphere system tremendously influence the terrestrial response to climate change. To our knowledge, however, inadequate attempts have been made to investigate the interactive changes in soil water-heat processes and the effects on vegetation across China in the context of intensified climate change. Therefore, we hypothesized that climate change would affect the spatiotemporal patterns of SM and ST regimes and result in ecological responses from the past to future decades featured by global warming. The objective was to contribute to understanding climatic changes in terrestrial conditions and the impacts on vegetation growth status. Here, we tested the hypothesis using SM and ST data sets taken from a replay of the land component of the European Center for Medium-Range Weather Forecasts climate reanalysis (ERA5-Land) and six models of the Coupled Model Intercomparison Project Phase 6 (CMIP6), along with remote sensing-based leaf area index (LAI) and CMIP6 simulations, for the future. We quantified the warming-drying and warming-wetting changes in the SM and ST regimes and the LAI responses. Section 2 describes the data sets and methods, and Section 3 shows the data evaluations. The results on the spatiotemporal changes in SM and STs are shown in Section 4 along with LAI responses from the present-day climate to future scenarios. Section 5 discusses uncertainties in the data sets, plausible mechanisms, and implications for eco-environments. Section 6 summarizes the paper with conclusions.

Reanalysis Data Sets of Soil Temperature and Moisture From ERA5-Land
The reanalysis data sets of ST and SM used to assess changes in soil water-heat regimes in the present-day climate were taken from the ERA5-Land data set (Muñoz-Sabater et al., 2021). To produce ERA5-Land, the land model, namely, the Hydrology-Tiled ECMWF Scheme for Surface Exchanges over Land (H-TESSEL) (Balsamo et al., 2011), was run at a finer spatial resolution of 0.1 latitude-longitude degree (approximately 9 km; in ERA5, it is 0.25° and 31 km) grid spacing, forced by meteorological fields from ERA5. Observations act only indirectly on the simulation through the forcings from ERA5. The soil was represented with four vertical layers with depths of 0.07, 0.28, 1.0, and 2.89 m below the surface. The monthly mean data were used for the 1980-2020 period in this study. To test soil water-heat processes from the perspective of a continuum, the variables of precipitation, runoff, surface air temperature, and sensible and latent heat fluxes were also evaluated using observations in this study.

Observational Data Sets of Soil Temperature and Moisture, Sensible and Latent Heat Fluxes, Runoff, and Leaf Area Index
In situ observations of SM and ST data were taken from the International Soil Moisture Network (https://ismn. geo.tuwien.ac.at) (Dorigo et al., 2011) and China Meteorological Administration (CMA, http://www.cma.gov. cn/2011qxfw/2011qsjcx), measured with gravimetric or sensor methods. The gravimetric measurements were performed on 3 days (8th, 18th, and 28th days) each month, and frozen periods and rainy days were excluded from sampling. The observational network known as CMA_Agr records measurements in relative SM based on field capacity to be easily used in agricultural practices; the CMA network, however, uses the form of volumetric SM. The SM data measured with sensors in the other four networks were also recorded as volumetric fractions. These in situ measurements were subjected to rigorous control of data quality. When the data sets were compiled the effects of sensor properties, accidental interference, and outliers were taken into account according to network-specific technical reports and the references therein (https://ismn.geo.tuwien.ac.at). For the present study, the available observations involve 956 sites during 1981-2016, and Table 1 shows detailed information on these data sets.
For the network CMA_Agr, the volumetric SM data were derived from observation-based relative SM according to the following relationship (Maidment, 1993): where SM v denotes the volumetric SM (m 3 m −3 ) and SM w denotes gravimetric SM (kg kg −1 ); β w and β d denote water and soil densities (kg m −3 ), respectively, and β w is normally taken as 1,000 kg m −3 ; SM r is the relative SM; and FC represents the field capacity with a value taken from a data set by Dai et al. (2013).
The 0-1.0 m observations (averaged throughout 10 layers, Table 1) from the CMA network were compared with the reanalysis means of the uppermost three layers. The 5 and 0.1 cm observations from the MAQU and SW-WHU networks were directly compared to the 0-0.07 m layer reanalysis data. The 0.1, 0.2, and 0.4/0.5 m observations from the CTP_SMTMN, HiWATER_EHWSN, and CMA_Agr networks were used to obtain the 0-0.28 m values, corresponding to the first two layers in ERA5-Land, by interpolation the same as that used by Li et al. (2005). Measurements were not conducted during frozen periods, except those from the CTP_SMTMN and MAQU networks on the TP affected by soil freeze-thaw processes.
The in situ observations of ST at depths of 0.15, 0.2, and 0.4 m below the surface were measured by CMA using soil thermometers, according to the observational requirements of the World Meteorological Organization. The observations are conducted at 2, 8, 14, and 20 o'clock (UTC+8) every day. Except for the sites with over 5 days of consecutive missing records and 30 days of total missing records, the data from 81 sites were used herein during the period 1980-2000. The interpolations for the 0-0.28 m layer were compared to those of ERA5-Land on a monthly basis.
The sensible heat and latent heat flux observations were from the Chinese flux observation and research network (Yu et al., 2006) and the integrated land-atmosphere interaction observation network on the TP . There were 9-site records from the first network during 2003-2010 available for this study. These measurements were subjected to data quality control and postprocessing and thus were commensurate with the other records compiled by the international FLUXNET sites (Yu et al., 2006). There were an additional 6-site observations from the second network on the TP measured during 2005-2016. The eddy covariance technique was used in both networks, and data quality was controlled following the steps by Foken et al. (2004) or using the TK3 software package in the TP network (Mauder & Foken, 2011). The ERA5-Land fluxes were evaluated using monthly mean data in this study.
The monthly mean streamflow data were measured at the Lijin and Datong sites for the Yellow River and Yangtze River basins, with data available for the periods 1980-2016 and 1980-2001, respectively. These records are referred to as natural streamflow and were obtained from measured streamflow data conducted by the respective river conservancy commission by estimating the human-caused upstream water loss (Dai et al., 2009;Lv et al., 2019). The ERA5-Land streamflow at the sites was estimated by adding up runoff across the upstream watershed at the monthly scale, with negligible influence induced by the lagged time of the flow confluence process.
The LAI data set was taken from the Global Land Surface Satellite (GLASS) Product Suite (Liang et al., 2021). The GLASS LAI product is synthesized from Moderate Resolution Imaging Spectroradiometer and Advanced Very High-Resolution Radiometer data, with long temporal coverage (1981-2018 available at the time), high spatial resolution (0.05°), and spatial continuities without missing pixels. This product (available from  www.glass.umd.edu) has been validated extensively using in situ observations and has been compared with other satellite remote-sensing data sets and used widely in ecological and climate change studies (Liang et al., 2021 and the citing references).

CMIP6 Projections Based on Three Socioeconomic Pathways
To assess future changes in soil temperature-moisture regimes with climate warming, we examined the two variables from all CMIP6 outputs for three shared socioeconomic pathway (SSP) scenarios (SSP126, SSP245, and SSP585). Since ERA5-Land can regenerate the observed terrestrial water and heat dynamics on various spatiotemporal scales (see Section 3.1), we evaluated and selected the CMIP6 ST and SM outputs by comparison with ERA5-Land on a regional scale. According to the pattern correlation coefficients (PCCs) of linear trends during 2000-2020, the CMIP6 historical simulations of ST and SM outputs were compared to those of the ERA5-Land reanalysis. Given both the temporal and the spatial signals involved, the PCC of linear trends acted suitably as a metric to evaluate the terrestrial water and energy dynamics of the CMIP6 outputs for the present-day climate.
Here, all CMIP6 outputs were selected with positive PPCs of SM. Table 2 shows the selected models in this study, and detailed model information can be seen online (at https://wcrp-cmip.github.io/CMIP6_CVs). At the time, to address the unavailability of layered ST and SM from the CAS-ESM2 and FGOALS-g3 simulations, the interpolation values were generated using those of column and 10-cm soil layers by the same practice as that described in Section 2.1.2.
The LAI projection data used in this study were from the outputs of three CMIP6 models, namely, MPI_ESM1-2-HR, TaiESM1, and UKESM1-0-LL, because the other models used a monthly climatology, that is, with only monthly variability but not interannual variability.

Ensemble Mean of CMIP6 Outputs
For the ST and SM outputs of the six selected CMIP6 models (Table 2), as well as LAI outputs, first, their data with identical resolution (0.25°) were generated by bilinear interpolation. The ensemble data sets were then yielded by equally weighted multimodel means, as commonly described in the IPCC sixth assessment report (Chapter 1). In this study, given that the interest lies in the long-term variations with respect to the historical period (2000-2014), we obtained ensemble data of the six CMIP6 outputs by two steps, separately calculated means and anomalies, and then integrated them into the ensemble data, holding equal individual long-term trends and variabilities.

Seasonal Decomposition
Given the scarce in situ observations of ST and SM, herein, the reanalysis counterparts from ERA5-Land were used as reference data, and thus, their fidelity was first assessed by comparison with available observations across China.  Of high importance was the long-term variation in this study. Thus, we integrated all available SM observations into a time series by averaging observations across sites along their time axes, constituting a subset of the statistical population for the real climate over the entire observational period, to assess ERA5-Land SM from monthly to decadal scales. Accordingly, the time series of the ERA5-Land SM was derived according to the corresponding timings and locations. The synthesized time series may not be equal to a real climate, but they were more robust and competent in terms of assessing the long-term variations in model simulations . The resulting monthly time series were then decomposed into long-term trends, seasonal cycles, and residual signals using a seasonal-trend decomposition method based on locally estimated scatterplot smoothing (Cleveland et al., 1990) as follows: where Y t indicates the synthesized observations or reanalysis time series. Y t is composed of the long-term trend (T t ), seasonal cycle (S t ), and the remaining (R t ) signal in an additive model at time t. Using linear least squares regression, the linear and nonlinear signals can be further derived from T t .

Statistical Analysis
The Pearson correlation coefficient (r) was employed to measure the consistency between observation and model data sets, and r is defined by the following equation: where COV(X, Y) is the sample covariance between the two time series of X and Y with sample size n, and they are normally distributed with means and and standard deviations s x and s y respectively. The p value method was used to measure the significance of r according to a t-distribution with the degrees of freedom equal to n − 2, and the test statistic t was calculated as follows: The PCC was obtained from the Pearson product-moment coefficient of linear correlation (Equation 2) between two arrays with the values of the same variables at corresponding locations. The PCCs were calculated directly (uncentered) according to the discussion in the IPCC report (IPCC, 2021).
Regarding the trend analysis, this study used the nonparametric Theil-Sen estimator (Theil, 1949), which is a robust linear regression method, and calculates the slope as the median of all slopes between pairs of values. In contrast to the ordinary least squares estimator, the Theil-Sen estimator has strong robustness to outliers. The significance of the trend was tested by the Mann-Kendall test (Mann, 1945).

Evaluation of Reanalysis Data
The evaluation of ERA5-Land is herein focused on terrestrial water and energy interactions, involving critical variables of SM, ST, sensible and latent heat fluxes, and runoff, and the goal was to assess the physical adequacy of water-energy processes in the model system against in situ observations.
At the point scale, the SM comparisons showed that ERA5-Land captured variations in the observations in terms of correlations ( Figure 1). Across 80% of the sites from networks 1-5 (Figure 1a), the coefficients were significant (p < 0.05) on the daily scale. Moreover, 86% of sites from network 6 were significant on the monthly scale ( Figure 1b). These sites were distributed evenly across the networks, especially in the regions to the east of 100°E, showing reasonable spatial representativeness of the adequacy of ERA5-Land SM. Notably, on the TP (networks 2 and 4), strong correlations appeared, indicating the reasonability of the major dynamics of SM even with soil freeze-thaw processes. The low correlations on the point scale were largely attributed to the lower spatial representativeness of in situ observations with respect to the large model grid cells (9 × 9 km), considering the complex and heterogeneous local soil properties, topography, and climates.
To depress the adverse influence of such a spatial mismatch on evaluations, the observations from two networks (CMA and CMA_Agr), with relatively long and complete records, were integrated into two time series, and then the ERA5-Land SM was evaluated using a method of seasonal decomposition ( Figure 2). The comparisons showed correlation coefficients up to 0.66 and 0.77 (p < 0.01) for the two networks on the monthly scale (Figures 2a and 2e). The seasonal cycles represented the most skillful signal in ERA5-Land SM, with correlation coefficients of 0.73 and 0.80 and comparable averaged seasonal cycles (Figures 2c and 2g). The long-term signals were also in good agreement with each other, although they had lower correlation coefficients (0.47 and 0.72, both at the p < 0.01 significance level), with the same sign in their linear trends over the two periods (Figures 2b  and 2f). The high-frequency variations also showed high correlations (0.72, 0.74, Figures 2d and 2h).
Regarding ST, ERA5-Land showed high skill levels compared with observations ( Figure 3), with correlations greater than 0.5 at 75 of 81 sites even after averaged annual cycles were removed. The lower correlations were distributed on the TP or in western mountainous areas (Figure 3a). The integrated time series across all sites agreed well with the observed behaviors on the monthly scale, and the correlation coefficient was as high as 0.97 for Network 6, the data are at a monthly time interval. The colored bars indicate the correlation coefficient distribution, and orange denotes significance (p < 0.05) as the minimum sample size n is equal to 30, accounting for 33% and 50% for (a) and (b); red indicates a significant correlation with values larger than 0.5, accounting for 47% and 36% for (a and b), respectively.

of 24
(p < 0.01, Figure 3b). Although slightly lower, especially during warm seasons, the mean annual cycle markedly captured the temporal structure of the observations (Figure 3c).
Apart from SM and ST, the comparisons of the other components of terrestrial water and energy processes further showed reasonable levels of realism in ERA5-Land data. For instance, the sensible and latent heat fluxes were in good agreement with the observations with respect to annual cycles averaged across the networks, with correct peak timing across eastern sites, indicating reasonable water and heat partitioning over the region. On the TP, however, 1-month lagged peaks appeared in both fluxes (Figure 4). With regard to runoff variations, ERA5-Land can regenerate the measured variations in the Yellow River and Yangtze River basins on a monthly scale, with correlations of 0.35 and 0.82 that were significant at the p < 0.05 and <0.01 levels, respectively ( Figure 5).
The preceding evaluations were justified using ERA5-Land as a reference data set in terms of the physical adequacy of water and energy processes. Nevertheless, the adverse influence of wet biases in ERA5 precipitation in China, especially before 2000, on the long-term trends of SM and ST was considered (Jiao et al., 2021; Z. Q. Ma et al., 2022). Hereafter, we focused on the period 2000-2020 to examine soil temperature-moisture change regimes. CMA_Agr networks and the corresponding reanalysis of the ERA5-Land for the 0-1.0 and 0-0.28 m soil layers, respectively. Here, r denotes the correlation coefficient, and the superscript **indicates significance at p < 0.01 with respect to sample sizes n = 228 and 253 in the two networks, respectively.

Evaluation of CMIP6 Historical Simulations
Given the high adequacy of water and energy balances, the ERA5-Land data set was used to evaluate the six CMIP6 model ensembles ( Figure 6). In terms of the spatial pattern of linear trends during 2000-2014, the ensemble SM was in good agreement with that of ERA5-Land, with a significant PCC at the p < 0.05 level (up to 0.2, n = 3,816). Notably, the strong drying trends in North China were well captured by the ensemble. Furthermore, from northwest to southeast, the spatial structure of decreasing-increasing-decreasin g-increasing trends was regenerated by the ensemble SM, despite the lower intensity overall. For the ST ensemble, the warming trends dominated the same period across China, consistent with those in ERA5-Land (PCC = 0.21, p < 0.05). The local cooling trends not generated in the ensemble would cause fewer large-scale changes in soil temperature-moisture changes. Taken together, the ensemble ST and SM data sets regenerated the main variations in their reanalysis counterparts in the present-day climate.

Changes in Soil Temperature and Moisture During 2000-2020
In the present-day climate, changes in ST and SM in the 0-1.0 m layer were examined using the ERA5-Land data on the monthly scale with seasonal cycles removed from 2000 to 2020 (Figure 7). Compared with the significant (p < 0.05) soil warming across China, the SM changed unevenly in space, decreased in the arid northwest, and increased in the semiarid regions from the northeast to the north of the TP. From North China to Southwest China, the SM significantly decreased, and then increased across the Southeast.
These interactive changes between the univariates fostered the regional soil temperature-moisture change regimes. In North China, the soil energy-water changes fall into the warming-drying regime, with the strongest trends over −0.03 m 3 m −3 K −1 in the period 2000-2020. This regime extends southwest and dominates the regions where major droughts have increasingly occurred during recent decades. Similarly, in the western TP and the arid northwest, warming-drying regimes dominate energy-water changes, with mean trends of 0.018 and 0.008 m 3 m −3 K −1 , respectively. In contrast, the warming-wetting regime appeared in the mountainous regions of the northwest, the eastern northeast, and parts of the TP and the south. To the west of the northeast, the cooling-drying and cooling-wetting regimes appeared locally. On average, the areas of the warming-drying, warming-wetting, cooling-drying, and cooling-wetting regimes were approximately 24.5%, 6.4%, 0.08%, and 0.03% of the land area of China, respectively. Because of the requirement of significance in both changes, large-area regions, for instance, Northeast and South China, fall into a nonsignificant change regime.
From the perspective of seasonality (Figure 8), the strongest regime, namely, the warming-drying changes, stood out during summer (June-August) and autumn (September-November). The maximum trend was −0.5 m 3 m −3 K −1 in North China. In the southwestern regions, the strong warming-drying regime endured from spring (March-May) to summer, with the steepest trend of −0.6 m 3 m −3 K −1 in spring. In terms of coverage, the four change regimes extended to their maximums in spring, covering approximately 20.5% of the land area of China. During winter (December-February), they contracted to their minimums, with approximately 4.0% of the land area, and summer and autumn were in between, with 12.8% and 10.6% of the land area, respectively. The areal percentages for the warming-drying regime were 2.1%, 17.0%, 10.8%, and 7.6% in winter, spring, summer, and autumn, respectively, and those for the warming-wetting regime were 0.3%, 3.5%, 1.4%, and 2.5%, respectively. It is worth noting that in the northwest, the cooling-wetting regime occurred in winter, and   transformed into a warming-wetting regime in the deeper layers. That is, the warming-drying regime contracts or weakens with increasing soil depth. Nevertheless, the warming-drying regime characterized the main changes in the ST and SM in response to present-day climate change and variability.

Ecological Effects of Soil Temperature-Moisture Changes During 2000-2018
Responding to the changes in ST and SM, ecological changes indicated profound effects of climate variability and change. The joint distributions of ST and SM with the LAI (Figure 10) indicated that the warming-drying regime adversely, but the warming-wetting regime favorably, affected the LAI variations during 2000-2018. Based on the regression analysis, across the regions dominated by warming-drying changes, the LAI decreased on average by 16.5%, 6.2%, 12.2%, and 3.7% with respect to that in 2000 in northwest (NW), North China (NC), the TP, and southwest (SW), respectively. The NE, SW, and SC regions were mainly covered by forestlands with high absolute LAI values. From the perspective of absolute variations (Table 3), the LAI in SW was more sensitive to such a soil temperature-moisture change, followed in turn by NC, NW, and TP. For the regions controlled by the warming-wetting regime, the LAI increased by 9.3%, 2.9%, 19.2%, and 3.0% in NW, NE, TP, and SC, respectively. Based on absolute variations, the LAI was more susceptible to warming-wetting changes in SC, followed by those in NE, TP, and NW (Table 3). The correlations of the LAI with SM and ST under the warming-drying regime had positive and negative values, respectively (Table 3). Their absolute values indicated a greater control on LAI variation by SM than ST variations. In contrast, under the warming-wetting regime, the effects from ST were positive, but both controls decreased, implying that water-heat stresses decreased.
The annual cycle changes in the anomalies of SM, ST, and LAI throughout 2000-2018 ( Figure 11) showed distinct interactions among the three variables across regions. In NW, during the growing season, the LAI decreased from 2000 to 2018 and was even faster in the early time, followed by a faster rise in spring ST and a consistent decrease in SM (WD regime, Figures 11a-11c). As soil experienced warming-wetting changes (WW regime), the LAI increased quickly in the early months and decreased slowly in the later months in this region (Figures 11d-11f). In NC, the warming-drying regime tended to inhibit leaf area developments around April and August but favored June, with the dynamics of ST and SM (Figures 11g-11i). The warming-wetting regime in NE led to LAI increases during the later months of the years because of the positive anomalies in SM and the lagged ST drop in autumn (Figures 11j-11l). The response curves of the LAI on TP were analogous to those in NW under both warming-drying and warming-wetting regimes (Figures 11m-11o and 11p-11r). In contrast, across southern China, following the quick drops or accumulated variations in SM, the LAI decreased in spring and autumn (April, September, and October) in SW (Figures 11s-11u) and increased in SC (Figure 11v-11x). Taken together, in terms of the averaged ranges of the annual LAI anomalies (0.24 and 0.21 for WD and WW regimes, respectively), the warming-drying regime exacted stronger impacts on ecosystems than did the warming-wetting regime under present-day climate change regardless of the ecosystems.

Projected Changes in Soil Temperature and Moisture During 2015-2100
Following three future socioeconomic pathways, interactive changes in ST and SM in the 0-1.0 m layer showed high sensitivity to the emission scenarios and phases ( Figure 12). The spatial structure of the present-day climate  persisted to an extent during 2015-2039 in that of the SSP585 scenario, particularly the warming-drying regime in North China (Figure 12a). Entering the middle of the twenty-first century (2045-2069), the warming-wetting regime expanded from the northwest to the southeast, with a 9.9% increase in its areal percentage (from 32.1% to 42.0% land area of China) at the expense of a 9.5% decrease (from 30.7% to 21.2%) in the area of the warming-drying regime. During the last 25 years in the twenty-first century (2075-2099), the consistent southeast expansion of the warming-wetting regime (55.5%) featured intensification in North China and a remarkable areal contract of the warming-drying regime (8.8%). Overall, throughout the twenty-first century (2015-2100), ST and SM jointly trended toward warmer and wetter states across China under the SSP585 scenario, although it was slightly warmer but drier in the southeast and on parts of the TP.
Under the medium pathway (SSP245), multidecadal shifts characterized the change regimes, with 26.6%, 32.2%, and 25.3% land areas for the warming-drying regime and 38.8%, 33.0%, and 27.5% for the warming-wetting regime during the early, middle, and latter 25-year periods, respectively. However, the eastern semiarid region was persistently influenced by the warming-drying regime, and this regime also tended to have a leading influence on the Yangtze River basin and its southern regions from the perspective of long-term trends.
Under the low pathway (SSP126), the two predominant regimes (warming-drying and warming-wetting) attenuated from the first to the last 25-year phases in the twenty-first century (from 34.4% to 10.6% and 37.3% to 1.6%, respectively). After the emission peak in the 2040s, they were partly replaced by the cooling-drying and cooling-wetting regimes across the semiarid regions (4.2% and 9.3%), although the warming-drying regime intensified in the southeast. However, in terms of long-term changes (2015-2100), the warming-drying and warming-wetting regimes featured the soil response to climate change under the low emission scenario.
Despite the distinct spatial shifts over the three phases, the long-term soil temperature-moisture changes were dominated by a northwest-southeast transition from the warming-wetting to warming-drying trends under three socioeconomic pathways (with 71.4%, 39.9%, and 54.3% of land areas for warming-wetting and 19.3%, 47.0%, and 25.6% of land areas for warming-drying trends during 2015-2100).

Ecological Effects of Soil Temperature-Moisture Changes During 2015-2100
Because the ecosystems in North China and the southwest are experiencing strong adverse impacts of present climate change, we thus focused on these two regions to assess the ecological impacts of soil temperature-moisture change under future scenarios ( Figure 13). The dynamics of ST, SM, and LAI showed distinct responses of LAI to soil temperature-moisture changes under the three scenarios. In North China, following increases in both ST and SM, the LAI showed a remarkable increase (up to 1.1) over 2015-2100 under SSP585, with decadal fluctuations (standard deviation, s = 0.3); for example, the decrease during 2025-2035 was consistent with that during 2000-2020 ( Figure 10). Given its positive correlation with SM (r = 0.79, p < 0.01), these fluctuations were more forced by soil water stress. Under SSP245, with the extended soil warming-wetting changes (Figure 12), the SM anomaly fluctuated more (s = 0.003), and the LAI showed a slower uptrend (slope is 0.004), especially after the 2040s. With negative moisture anomalies and slower rising temperatures under SSP126, the uptrend in the LAI shifted to a downtrend after the 2050s in North China.
The southwest was influenced by the warming-drying regime under the three scenarios, especially from the 2050s on. The LAI increased consistently except for the decadal fluctuations during the entire period under SSP585 and SSP245. The limited negative SM anomalies had yet to force a decrease in the LAI. However, entering the last 25 years in the twenty-first century, the LAI was increasingly depressed by water stress along with intermittent negative anomalies in ST. Overall, future soil climate changes favor ecological developments in the two typical regions in terms of long-term trends despite the adverse fluctuations on decadal scales from the medium to high emission scenarios. In contrast, under the low emission scenario, the developments are depressed during the second half of the twenty-first century.
In light of the annual cycle anomalies (Figure 14), the soil presented warming-wetting trends in the SSP585 scenario but warming-drying trends in the SSP245 and SSP126 scenarios from the 2020s to the 2090s in North China (Figures 14a, 14g and 14m), with peak temperatures in spring and autumn (Figures 14b, 14h and 14n). With overall positive anomalies, the LAI increased far more strongly in the first than in the second half of the years in response to the lower water stress and advancing soil warming in early spring (Figures 14c, 14i and 14o).
In the southwest, with soil warming-wetting changes, especially in spring and autumn months, the LAI showed increasing peaks around May and after November.   11. Temporal evolutions of annual cycles of soil moisture, soil temperature, and leaf area index anomalies (with mean annual cycles removed) in four typical regions (same as Figure 10) under the warming-drying and warming-wetting change regimes during 2000-2018.
With future changes in soil water and heat conditions, the enhancement of the LAI in spring and autumn has intensified remarkably over time in the twenty-first century. Moreover, for the two most important ecosystems, grasslands and forestlands, such LAI responses to soil water and heat changes were similarly true across China under the three scenarios (figures not shown).  -2039, 2045-2069, 2075-2099, and 2015-2100 under three shared socioeconomic pathways (SSP126, SSP245, and SSP585).

Associations Between Soil-Atmospheric Moisture and Temperature Changes
Based on ERA5-Land ST and SM during 2000-2020, approximately 24.5% of the land area experienced warming-drying changes in the soil, especially affecting North China and southwestern China. In the context of an ever-increasing warming climate with decreasing precipitation during the past 21 years, meteorological droughts were frequent, as were concurrent heatwave events (CMA, 2021;Kong et al., 2020;Li & Chen, 2020;Xu et al., 2022). Given that the two regions are in hotspot regions of land-atmosphere coupling (Li et al., 2017), drying trends in the atmosphere can largely propagate into the soil and thus lead to soil warming-drying changes.  In Northwest China, the western TP, Northeast China, and South China, increasing precipitation provides moisture and results in soil warming-wetting changes (Li & Chen, 2020). However, in the context of warming-wetting changes in climate, especially in Northwest China (Zhang et al., 2021), the significant increases in both potential and actual evapotranspiration lead to soil drying trends in parts of those regions (Su et al., 2022). As a result, soil temperature-moisture changes are not only modulated by atmospheric anomalies but also associated with terrestrial hydrological and ecological processes (Zhang et al., 2021).

Intensified Soil Temperature-Moisture Changes With Soil Depths
The areas that experienced significant changes in ST and SM expanded with soil depths from 21.9% to 24.7%, 31.0%, and 41.2% of land areas in the 0-0.07, 0-0.28, 0-1.0, and 0-2.89 m layers, respectively.  revealed substantially increased persistence of SM anomalies with soil depth in transition climate zones. Accordingly, the expanded changes in deep soil were mainly contributed by more persistent SM anomalies than those in the upper layers. Such changes in deep soil imply that soil may amplify atmospheric anomalies (e.g., droughts in drylands) and thus prolong and intensify adverse influences on overlying ecosystems under dry climates.

Contrasting Changes of Warming-Wetting Northwest and Warming-Drying Southwest
Under three climate scenarios from 2015 to 2100, regions of 30.6% and 55.2% land area on average experience significant soil warming-drying and warming-wetting changes, respectively; that is, in China, soil gradually undergoes warming-wetting changes from northwest-southeast but warming-drying changes in the southeast. Not unexpectedly, Wang and Zhai (2022) reported warming-wetting changes in Northwest China, especially under SSP585, based on 16 CMIP6 model projections. The patterns of a wetting northwest and a drying southeast also agreed well between CMIP5 and CMIP6 precipitation projections (Du et al., 2022); the more intensified the climate warming is, the more contrasting the precipitation patterns are (e.g., Figure 4.32 in IPCC (2021)). Accordingly, the nationwide characteristics of soil wetting-drying changes are, in effect, in agreement with precipitation trends.

Possible Mechanisms for the Ecological Effects of Soil Temperature-Moisture Changes
Concerning ecological effects, during the past 21-year period, soil warming-drying and warming-wetting changes jointly exerted a significant influence on the overlying ecosystems. For example, the decreasing LAI trends in warming-drying North China and the southwest and the increasing trends in the warming-wetting northeast and northwest mountainous regions are consistent with previous reports (Chen et al., 2019;Piao et al., 2020). According to theories on plant water use and temperature stress, soil water and heat conditions directly regulate plant physiological and biochemical activities and the ability to cope with environmental stresses (Chaves et al., 2002;Hasanuzzaman et al., 2013;Lian et al., 2018;Nelson et al., 2020). The soil warming-drying changes may stress LAI development in the main growing seasons (e.g., May-September) once SM decreases toward wilting moisture. Importantly, in spring, the LAI shows resilience to drying changes and increases earlier ( Figure 11) in autumn, the LAI decreases later, that is, there is lagged senescence, for warming soil and climate. These changes in vegetation phenology partly contribute to the vegetation greening observed in China because climate warming eases spring and autumn constraints on temperature-limited ecosystems (Chen et al., 2019;Piao et al., 2020). Vegetation process changes, as a subsystem of land-atmosphere interactions, reflect both the response and the feedbacks to soil water and heat changes with climate change.

Ecological Implications of Future Soil Temperature-Moisture Changes
Under the three future pathways, annual-decadal fluctuations in the LAI, indicating vegetation impairments, are substantial in both typical regions and ecosystems (Figure 13), resulting largely from the stresses from soil drying and warming. In contrast, in terms of long-term trends, the LAI consistently increased as the ST increased. The spring-month increases were much larger than those in the other growing months in North China, highlighting the trends of advanced phenology with climate warming. In the southwest, the autumn increases were comparable to those in spring. Although vegetation would reduce photosynthesis as ambient temperature exceeds 32°C, the extended growing seasons, among other favorable factors (e.g., carbon dioxide fertilization effect globally explains the largest contribution to the satellite-observed LAI trend from 1982 to 2009 (Zhu et al., 2016)), compensate for the reduction and further enhance the LAI in the future. These projected changes in the LAI in the present study agree with the estimates for 2021-2050 using geographically weighted regression based on a regional simulation driven by 5 CMIP5 projections (Gao et al., 2017;Shi et al., 2018). The results from multiple linear regression with 11 CMIP6 projections similarly reported that the northeast is likely to become warmer and wetter, and hence have an increase in regional LAI (Yuan et al., 2021). Likewise, throughout the twenty-first century, vegetation in China has a greening trend according to another study with 12 CMIP5 models . Taken together, the vegetation-soil interactions under climate warming are projected to favor vegetation greening in the twenty-first century. Nevertheless, multiyear drought events in the future are projected to be intensified with longer durations (Rakovec et al., 2022), and thus, we should be prepared for their impacts on soil-ecosystem systems and ecosystem services.

Limitations and Uncertainty
Note that the in situ observations of ST and SM, sensible and latent heat fluxes, runoff, and remote sensing retrieval of LAI data sets have a variety of spatiotemporal representativeness. As a result, the corresponding statistics may indicate unequal skill levels of the evaluated reanalysis data. Furthermore, the observed objects are always influenced by human-induced processes, for instance, land use changes and agricultural irrigation. These processes have yet to be adequately represented in Earth system models and reanalysis systems. In addition, the spatial scale mismatch between point-scale observations and model grid data may weaken the effectiveness of model evaluations based on in situ observations.
The present study summarized the skill levels of CMIP6 models according to the spatial patterns of linear trends in ST and SM in comparison with their reanalysis counterparts. It is noteworthy that Earth system models are physically complex systems, and thus, the assessment of CMIP6 projections simply in terms of SM and ST is inadequate due to both the fidelity of reanalysis data and the availability of observations. For this reason, different model groups likely contribute to discrepancies in soil water and heat changes (Srivastava et al., 2020). Nevertheless, nationwide soil warming-drying and warming-wetting changes and thus ecological impacts with climate change deserve attention in climate impact research, as well as model development, especially from the perspective of the systematic adequacy of land-atmosphere interactions.
The representations of water and heat balance processes, as well as vegetation activities, remain a challenging task in the CMIP6 models. For example, the absence of groundwater, runoff routing, and irrigation processes and the inadequacy of plant function types and their sensitivity to climate warming would add uncertainty to the understanding of soil-vegetation interactions under future climate scenarios.

Summary and Conclusion
Global warming modifies the hydrological cycle and consequently the regional soil properties through combined soil temperature-moisture changes. Joint soil temperature-moisture changes directly affect ecology and ecosystem services, which is particularly important for plant stress-sensitive regions such as China. This paper systematically investigated the spatiotemporal characteristics of joint soil temperature-moisture trends and regional ecological responses across China during the recent past and projected future. The main findings are as follows: In the past 20 years, soils in over 24.5% of land areas in China have experienced a warming-drying trend, and soils in over 6.4% of land areas have experienced a warming-wetting trend. Prominently, a warming-drying regime consistently appeared in all seasons and expanded with soil depth. These opposing trends in soil properties have led to widespread ecological changes, corresponding to a 9.7% decrease in LAI in drying areas and an 8.6% increase in LAI in wetting regions. These are typical regions covered dominantly by grass, crop, savanna, and forest ecosystems. Under future scenario pathways from 2015 to 2100, the area of soil warming-drying expanded to 30.6% and that of warming-wetting expanded to 55.2%, on average. Warming-wetting regions expanded from northwest to southeast at the cost of contracting warming-drying regions, despite significant spatiotemporal variabilities on annual to decadal timescales. Following combined soil temperature-moisture changes, the LAI also showed substantial annual to decadal scale fluctuations. The long-term increasing trends of the LAI remained the leading signal, benefiting largely from the increase in growing-season length, among others.
Anthropogenic global warming alters the terrestrial water cycle; temperature increases generally extend the length of the growing season and lead to the greening of the Arctic and the northern extratropical land surface. In response to more substantial uncertainties regarding the subtropics, this study examined the systematic impacts of climate change on soil water-heat changes and thus ecosystems, providing insights into the responses of terrestrial ecosystems to climate change. Although focused on China, which stretches from the northern subtropics to extratropical latitudes with contrasting land environments, the findings from this study have wide implications for other regions around similar latitudes and climate characteristics. The results also help understand potential changes in terrestrial carbon sinks and environments under a warming climate. We have emphasized the importance of joint warming and hydrological changes in ecological responses through combined soil temperature-moisture trends. Methodologically, it is helpful to evaluate Earth system modeling from the perspective of physical consistency. Future work should entail improving the physical, physiological, and biochemical representations of terrestrial water cycles and vegetation water use processes in Earth system models.

Data Availability Statement
The reanalysis data sets of ERA5-Land are available online at https://www.ecmwf.int/en/era5-land. The CMIP6 projection data sets are available online at https://esgf-data.dkrz.de/search/cmip6-dkrz. The monthly observations of ST (Li, 2022a) and SM (Li, 2022c) used in this study are available from Zenodo at https://doi.org/10.5281/ zenodo.7024791 and https://doi.org/10.5281/zenodo.7030324, respectively. The monthly observations of sensible and latent heat fluxes (Li, 2022b) from the nine sites used in this study are available from Zenodo at https://doi.org/10.5281/zenodo.7024171 and the other observations from the sites on the TP are available from the Science Data Bank (Ma et al., 2020a) at https://doi.org/10.11922/sciencedb.00103. The monthly runoff data sets (Dai, 2017) are available at https://rda.ucar.edu/datasets/ds551.0, and the LAI data set from the Global Land Surface Satellite product suite is available at http://www.glass.umd.edu/LAI/AVHRR. The figures were made with Matplotlib version 3.3.2 (Hunter, 2007) on the Python 3.8 (https://www.python.org) platform, for which we would like to express our gratitude.