Global influence of soil texture on ecosystem water limitation

Linking stomatal regulation and soil–plant hydraulics

As the soil dries, its water potential, ψ (which is the negative work to extract a unit of volume of water), and its hydraulic conductivity, k, decrease by several orders of magnitude. Therefore, there is a critical soil water threshold (moisture or potential) at which the soil can no longer supply water to the roots at the rate required to sustain the transpiration demand^29,64. The relationship between water supply and demand provides the mechanistic link between soil–plant hydraulics and stomatal regulation in water and carbon exchange between vegetation and atmosphere^20,25. The nonlinear decrease in soil–plant hydraulic conductance with soil drying or increasing VPD has the key role in this framework. Note that different definitions and meanings of hydraulic ‘conductance’ and hydraulic ‘conductivity’ may create confusion. Here we refer to hydraulic conductance, K, as the ratio between the water flow, such as T (millimetres per day), and the water potential difference, Δψ (megapascal), yielding millimetres per day per megapascal, while we refer to soil hydraulic conductivity, k (metres per second), to the ratio of the water flux density, q (metres per second), to the gradient in soil water potential along the flow path, dψ/dx (metres per metre). It is believed that stomata downregulate transpiration and photosynthesis at the point where transpiration comes at a disproportional ‘cost’ to water transport^20,25,30,35. Despite the popularity and mechanistic strength of the framework, and models derived from this, it remains unclear which is the limiting hydraulic element of the soil–plant continuum. There is also discrepancy in the rules applied to derive stomatal regulation from the disproportionate cost of water transport (for example, centred on stomatal optimality, hydraulic conductance or physiological mechanisms). Here we used the supply–demand framework, implemented according to ref. ²⁵.

General model description

Critical soil water thresholds were calculated based on the hydraulic framework formulated in ref. ²⁵. The premise is that stomata downregulate transpiration when the relationship between transpiration and leaf water potential becomes nonlinear (that is, loss in hydraulic conductance). The onset of nonlinearity is the uppermost limit of transpiration that can be supplied by soil–plant water flow before stomatal closure restricts the transpiration rate and photosynthesis. Hence, θ_crit and ψ_crit are defined as the minimum soil moisture and water potential at which the soil–plant hydraulic system can supply water at the rate of the potential transpiration (T_pot). In other words, θ_crit and ψ_crit are defined where T_pot (here, 4 mm d^−1) is at the edge of the linear zone of the T (ψ_soil, ψ_leaf) surface (green zone in Extended Data Fig. 1a,b) intersecting the stress onset limit (SOL)²⁵.

The surface T(ψ_soil, ψ_leaf) is the physical space of plant water use, and the SOL delineates between the linear and nonlinear zones (yellow in Extended Data Fig. 1a,b), thereby defining the onset of water limitation²⁵. In wet soils, the surface is planar for a large range of transpiration rates (green zone). As the soil dries, the surface bends (brown zone) as the relationship between transpiration rate and leaf water potential for a given soil water potential becomes nonlinear (dotted black lines in Extended Data Fig. 1a). This nonlinearity corresponds to a substantial decline in the hydraulic conductance of the soil–plant continuum. The transition from the linear to the nonlinear zone—the SOL—is defined as the point where dT/dψ_leaf reaches 80% of its maximum (for each ψ_soil). This indicates a substantial loss of hydraulic conductance in the soil–plant continuum and is hypothesized to trigger stomatal closure²⁵. The SOL presumes a plant physiological optimization to minimize the trade-offs between gas-exchange benefits and hydraulic losses. The SOL is the uppermost limit of transpiration. It sets the maximum transpiration and corresponding stomatal conductance that plants could sustain under given soil water and VPD conditions. Obviously, stomatal conductance can be lower than this value, for instance, limited by light and elevated CO₂. In other words, our model does not aim to reproduce stomatal closure driven by factors other than hydraulic limitation, which are crucial in predicting stomal functioning below the SOL.

The model assumes that stomata progressively close when the hydraulic supply does not match the water demand—a process that has a time scale of minutes to hours. We used this model to predict the onset of water limitation during soil drying, a process that is comparatively slower and has a time scale of weeks. During this time, the plant hydraulics may change (particularly for grasses and annual crops). A key assumption of our analysis was that the relevant hydraulic variables changed proportionally, that is, T_pot, K_plant and L_root were assumed to change proportionally.

Modelling steps and parameter estimation

To test our hypothesis, that critical soil moisture thresholds show soil texture specificity on an ecosystem scale, we simulated soil water thresholds for each soil textural class (12 US Department of Agriculture (USDA) classes) and compared them to both flux-tower and sap flow-derived observations. We simulated θ_crit and ψ_crit by solely varying the soil hydraulic properties (Supplementary Table 2), while keeping plant and climate parameters constant (Extended Data Table 1). We justified the constant set of plant and atmospheric model parameters using the insignificant relationships of differences between observed and simulated θ_crit to site-specific latent heat fluxes (that is, to the absolute evapotranspiration rates determined by the climate of each site, T_pot) across climates and biomes (Supplementary Figs. 5 and 6). The few required model parameters were set to average values from the literature, except for the effective L_root. The potential transpiration per land surface, T_pot, was set to 4 mm d^−1 (during daytime)⁶⁵. The maximum plant hydraulic conductance was set as K_plant-max = T_pot/−ψ_leaf-max (mm d^{−1 }MPa^−1), which gave a leaf water potential of −1 MPa (ψ_leaf-max) when the soil was wet (ψ_soil ≈ 0 MPa) and transpiration was at a maximum, that is, T_pot. The plant water potential threshold (ψ_x*) at which the stem hydraulic conductance is reduced to 50% (approximately ψ_x50) was set to −2.8 MPa, based on reported values of xylem embolism in woody species⁶⁶, while the effective root and rhizosphere radii, and the slope of the decrease in K_plant with increasing water tension, were set to default values²⁵. The effective L_root (m m^−2), defined per land surface area, was the only fitting parameter, and was inversely estimated by fitting θ_crit over all soil textural classes. In other words, the difference between observed (θ_obs) and simulated (θ_sim) critical soil moisture was minimized over both datasets across all soils (least absolute deviations) by varying only L_root (where n = 149 is the number of θ_crit observations).

Main data collection

To test our hypotheses, we compared our model simulations to both eddy covariance and sap flow data across climates and biomes. To estimate critical soil moisture thresholds (θ_crit) from the eddy covariance data, we acquired daily data comprising the soil volumetric water content and the latent and sensible heat fluxes from the eddy covariance sites provided by Integrated Carbon Observation System (ICOS) (https://www.icos-cp.eu/), AmeriFlux (https://ameriflux.lbl.gov/) and FLUXNET (https://fluxnet.org/), all of which have undergone standardized quality control and gap filling⁶⁷ (we used both measured, quality flag = 0, and good-quality gap-filled, quality flag = 1, data). Only sites for which in situ estimates of soil texture were available were selected (n = 44). These sites either reported the soil textural class or provided the fractions of sand, silt and clay from which we could classify the soil texture based on the USDA soil texture classification system⁶⁸. Given the high correlation of critical soil moisture thresholds in the surface soil layer with the θ_crit observed in deeper layers³⁶, only the surface layer soil moisture was considered in our analysis. To estimate critical soil moisture thresholds (θ_crit) using sap flow data, we acquired the sapwood-area-based sap flux density (cm³â€‰cm^−2_Asw h^−1) and soil volumetric water content time series from SAPFLUXNET⁶⁹ (https://sapfluxnet.creaf.cat/) using the provided sapfluxnetr package v.0.1.4 (ref. ⁷⁰). Only sites that reported local estimates of θ (shallow soil layer), soil texture (USDA classification, and/or sand, silt and clay fractions) and soil depth along the sap flux density values were kept for further analysis.

Main data analysis

Eddy covariance data

The evaporative fraction was calculated as the ratio of the latent heat flux to the sum of the latent and sensible heat fluxes for each day^71,72. θ_crit was determined from the relationship of the evaporative fraction to θ by applying a regression between evaporative fraction and θ using a linear-plus-plateau model (lin_plateau.R function from ref. ⁷³). The good correlation between the onset of a decline in evaporative fraction and a decline in gross primary production⁴ justified the approach to interpret these data as a decline in transpiration, although bare evaporation can substantially contribute to the evaporative fraction. θ_crit is defined as the soil moisture at the breakpoint between the linear increase phase and the plateau of the model. As in ref. ³⁶, θ_crit was estimated from periods of soil moisture dry-downs during the respective summer seasons (that is, June–July–August for the Northern Hemisphere and December–January–February for the southern hemisphere). Soil moisture dry-downs were considered to be periods in which the soil moisture decreased consecutively for at least 10 days after a rain event^33,74. We could determine θ_crit for 36 out of the 44 sites (see ‘Main data collection’) using the dry-down definition from ref. ³⁶. For five sites, we determined θ_crit for periods beyond the summer season, but still applying the 10-day drying criterion. For two sites, this dry-down criterion also did not result in a θ_crit estimate, and thus we neglected the 10-day dry-down criterion, rather applying the summer criterion. Finally, we discarded one site for which no θ_crit could be determined at all. The remaining 43 sites formed the basis for all further analysis. We justified the different dry-down criteria by the high coefficients of determination (‘Summer’ versus ‘Full-criterion’: R²_adj = 0.97, n = 11; ‘10 days’ versus ‘Full-criterion’: R²_adj = 0.98, n = 10). The eddy covariance sites span around the globe, encompassing all continents (excluding Antarctica) (Supplementary Fig. 7).

Sap flow data

Similarly to ref. ⁷⁵, (sub)hourly sap flux density time series were aggregated to daylight averages (06:00 to 20:00) using daylight_metrics from sapfluxnetr⁷⁰. As for the eddy covariance data, we estimated dry-down periods as periods where the daily (24 h) averages of θ decreased for at least ten consecutive days (site level). After intersecting summer and dry-down periods, θ_crit was determined as for the eddy covariance data, but on a tree level (multiple trees in the site sharing the same environmental data) using the linear-plus-plateau model (now part of the soiltestcorr package v.2.2.0 (ref. ⁷⁶)), given that a positive linear slope was determined and that the breakpoint determination met the standard significance criterion (P < 0.05). Finally, 14 sites (Supplementary Fig. 8) and 106 trees (multiple tree individuals per site, either from the same or a different tree species) resulted in 106 sap flow-derived estimates of θ_crit, spanning six soil textural classes (Supplementary Table 1).

Soil hydraulic properties

The parameters of the soil water characteristics as a function of soil textural classes were taken from ref. ⁷⁷, which reported the mean and standard deviation. Because the saturated conductivity data in ref. ⁷⁷ were from another data source and without information on its variability, we took the values from ref. ⁷⁸ that provided a recent global data collection.

The relative importance of soil and plant hydraulics

Analysing the relative importance of soil and plant hydraulics is key to identifying the dominant controls on ecosystem water limitation. We approached this in two ways: (1) by comparing the simulated soil and plant hydraulic conductance as a function of soil texture (Fig. 1b); and (2) by comparing the differences in θ_crit variability between observations and simulations (Fig. 2c).

Simulating the physical space and transpiration downregulation (SOL) in each soil textural class allowed us to disentangle, by means of the soil–plant hydraulic model, whether the soil or plant hydraulics would have, in relative terms, a stronger impact on soil water thresholds. We calculated the soil and plant hydraulic conductance as K_soil = T/(ψ_soil − ψ_soil-root) and K_plant = T/(ψ_soil-root − ψ_leaf), respectively, where ψ_soil-root and ψ_leaf are water potentials at the soil–root interface and in the leaves, respectively. To analyse the θ_crit variability, we quantified the variation in θ_crit in response to variations in soil hydraulic properties for each soil textural class. First, we determined, for each soil textural class, the minimum and maximum values of a hydraulic property by subtracting or adding the standard deviation from the mean value (the geometric mean, in the case of saturated hydraulic conductivity and the shape parameters of the soil water characteristics curve). Next, we defined the hydraulic properties of the ‘coarse end’ of a soil textural class by combining the minimum air entry value, maximum slope parameter of soil water characteristics curve and maximum hydraulic conductivity for each soil textural class. For the ‘fine end’ of a soil textural class, the maximum air entry value, minimum slope parameter and minimum soil hydraulic conductivity values were chosen. Note that τ, which is the slope of k_soil over ψ_soil when both are expressed in logarithmic scale, is positively correlated with K_sat and inversely correlated with the air entry value, h_b (the correlations between log(h_b) and τ and log(K_sat) and τ in the data presented in ref. ⁷⁷ are 0.78 and 0.88, respectively). To test the effects of variable plant traits (L_root density and plant vulnerability) and atmospheric conditions (increasing VPD) on soil water thresholds, we modelled θ_crit and ψ_crit by varying L_root from 1/30 (minimum) to 30 (maximum) times the reference, and ψ_x* from −1.5 (minimum) to −5 MPa (maximum). From these results, we calculated the span of θ_crit (S_θcrit, unitless) for each soil textural class as S_θcrit = max(θ_crit) − min(θ_crit), which allowed us to compare the effects of varying soil and plant properties to the variance of the observations—that is, the span of the FN and SFN observations per soil textural class. Varying atmospheric conditions were simulated using future projections of potential transpiration rates. Details of the future climate modelling are described in the section ‘Global map’.

The relative importance of VPD and soil moisture in ecosystem water limitation

To evaluate the relative importance of VPD and soil moisture in ecosystem water limitation, we compared the soil-specific simulations with observed ecosystem fluxes for two contrasting soil textures (median of five eddy covariance sites for clay and sand, respectively). For these simulations, we assumed that T_pot = 4 mm d^−1, corresponding to VPD = 1.5 kPa (Fig. 3a–d). Our model predicted that, at VPD = 1.5 kPa, plants could transpire at full stomatal opening (g_cmax) as long as the soil moisture was higher than, or equal to, the simulated θ_crit (that is, moving horizontally in Fig. 3). Rising VPD (that is, moving vertically in Fig. 3) triggered stomatal closure at a critical VPD, which was set by the stress onset limit (yellow line in Extended Data Fig. 1). The critical VPD declined with decreasing soil water content, but it remained relatively constant for θ > θ_crit (particularly in sandy soils; in Fig. 3, this critical VPD is approximately 2 kPa). This critical VPD in wet soil depends on plant hydraulics and T_pot. More precisely, critical VPD depends on the difference between T_pot/K_plant and the critical leaf water potential where K_plant declines (psi_star, Supplementary Information). For instance, transpiration would become water-limited at a low VPD (and high soil moisture) when the plant hydraulic conductance was too low to sustain high transpiration fluxes. In this case, the system becomes water limited even in wet soils due to plant limitations, the key driver being the rising VPD.

For the observed evaporative fractions in clay and sand, we used (half-)hourly fluxes from these eddy covariance towers and filtered them as follows, referring to previous studies^1,2, to remove unmeaningful data for our analysis: only positive latent and sensible heat fluxes; only during daytime; without negative soil moisture and VPD values; sufficient incoming radiation and VPD to drive substantial transpiration (photosynthetic photon flux density of more than 500 μmol m^−2 s^−1, VPD of more than 0.5 kPa); sufficient wind speed (more than 1 m s^−1) to foster vegetation–atmosphere coupling; and without ‘cold’ days limiting plant metabolism (where the median daily temperature was less than 15 °C). Because we aimed to analyse as much of the VPD–soil moisture space as possible, we used the maximum available time resolution (half-hourly to hourly) and all levels of gap-filled data. Binning and visualization of the evaporative fraction along VPD and soil moisture (thereby removing soil moisture values below and above the residual and saturated water content of each soil texture, respectively, as well as cutting off the extreme VPD conditions (VPD of more than 5 kPa) occurring in one sand site) revealed soil-specific responses to the two environmental drivers. Simulations of transpiration rate with soil drying, based on the SOL in each of the two soil textures, were anchored to the VPD axis by the assumption that the maximum transpiration rate being sustained by the underlying soil–plant hydraulic constraints corresponded to the experimentally observed maximum evaporative demand (we took the 99th percentile of the median VPD distribution of the five sites) in wet soil (θ > θ_crit) for each textural class. Inset plots displaying the median relative sensitivity of evaporative fraction to VPD and θ confirmed the stronger relative contribution of soil moisture limitation in coarse-textured soils.

Global map

To test the effects of variable atmospheric conditions and evaluate the impacts of future climate on θ_crit, we simulated θ_crit in all soil textural classes under changing potential transpiration demands. For each soil textural class, we obtained a numerical function of T(θ), as in Supplementary Fig. 9 (note that the slightly different model parameters enabled projections of future transpiration rate up to the maximum increase in T_pot, that is, +65%). An analytical sigmoidal function was fitted to T(θ) to calculate θ_crit and Δθ_crit for the expected changes in potential transpiration rate under future climate (2060–2069) (Fig. 4 and Extended Data Figs. 3 and 4). To estimate the potential transpiration rate, we chose a simple approach based on air temperature and relative humidity, as described by Ivanov’s formula^79,80. Temperature and relative humidity data for the years 2005–2014 (current climate) and 2060–2069 (future climate) were downloaded from World Climate Research Programme Coupled Model Intercomparison Project 6 (SSP2-4.5 scenario) using the EC-Earth3 model with a spatial resolution of 0.7° (refs. ^81,82). Additionally, current and future precipitation data were acquired to classify world regions based on the current and future AI. The AI was calculated on an annual basis as AI = precipitation/T_pot. Next, the global map of Δθ_crit resulting from changing climate (T_pot) was calculated as follows: in a first step, global maps of sand and clay content from SoilGrids⁸³ were used to determine a global map of soil textural classes (Extended Data Fig. 4a). For each pixel, the change in potential transpiration rate (Extended Data Fig. 4c) and the corresponding change in critical water content were then computed and visualized (Fig. 4 and Extended Data Figs. 3 and 4).

Statistical information

Functional relationships between key variables were underpinned by linear regression analyses using ‘lm()’ from the stats package (v.4.3.2) of R statistical software v.4.3.2 (ref. ⁸⁴). The linear regression prerequisites were verified using ‘check_model()’ from the performance package (v.0.11.0) (ref. ⁸⁵) (see ‘Code availability’ for details). As a goodness of linear regression fit, adjusted R²_adj values and two-sided P values of the linear regression slopes were calculated and are indicated in the figures, where meaningful (significant P values are indicated in the plot area using standard significance levels of *P < 0.05, **P < 0.01, ***P < 0.001, not significant (NS) P > 0.05). The main linear regressions, indicated by R²_adj and significance levels in the main text figures, are detailed as follows: Fig. 2a, solid black regression line, R²_adj = 0.34, P < 0.001, y = −0.0017x + 0.23; Fig. 2b, solid black regression lines encompassing all textures, R²_adj = 0.35, P < 0.001, y = 0.61x + 0.05, and excluding clay soils, R²_adj = 0.48, P < 0.001, y = 0.98x + 0.01; Fig. 2c, solid black regression line, R²_adj = 0.4, P < 0.01, y = −0.0019x + 0.23, dotted brown regression line, R²_adj = 0.87, P < 0.001, y = −0.0018x + 0.24, dotted dark green regression line, R²_adj = 0.87, P < 0.001, y = −0.0013x + 0.15; Fig. 2d, solid black regression line, R²_adj = 0.37, P < 0.001, log₁₀ (y) = −0.048x + 2.6 (see ‘Code availability’ for further details). Note that, in Fig. 2b, we tested whether the 95th intervals of the estimated slopes and intercepts of the linear regressions overlapped with the 1:1 line (see ‘Code availability’ for details). The error bars for the sand fraction (in Fig. 2a) show the entire range of sand percentages within each soil textural class⁶⁸. The data distributions are displayed using grouped box plots (the thick solid horizontal line represents the median, the lower and upper hinges correspond to the first and third quartiles, the whiskers extend to the highest or lowest value, respectively, but no further than 1.5 times the interquartile range, and the widths of the boxes scale with the square root of the number of observations in each soil textural class) and individual observations in the groups are displayed as points along the boxes. Note that the grouped (FN and SFN) boxes and points (Fig. 2a,b) are slightly shifted around the true x coordinate (the same for both) for readability. The number of sites in each soil textural class is given in Supplementary Table 2, while the number of θ_crit observations per soil textural class are displayed in Fig. 2a,b.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.