Hydraulic Redistribution Decreases with Precipitation Magnitude and Frequency in a Dryland Ecosystem: A Data-Model Fusion Approach

Hydraulic redistribution (HR), the movement of water via plant root systems that connect soil compartments with different water potential, should influence soil moisture dynamics particularly in dryland ecosystems, where water availability strongly constrains ecosystem function. Realistic representation of HR in ecosystem models is essential to improve the ability of these models to predict ecosystem function in dryland regions. In this study, we integrated HR into the Terrestrial ECOsystem model and employed a Bayesian Markov Chain Monte Carlo technique to optimize soil hydraulic parameters and root conductance using four years of soil moisture observations from a piñon-juniper woodland. We found that (i) integrating HR generally improved model prediction of soil moisture during dry periods, particularly in the top 30 cm of the soil profile, where more than 50 % of root biomass exists; (ii) HR increased surface soil moisture by up to 60 % during dry periods; (iii) HR decreased with increasing precipitation magnitude and frequency, however, the length of dry spells between rainfall events also influenced HR rates; and (iv) upward HR in the top 60 cm soil profile became more pronounced as dry conditions progressed, with rates ranging from 0.10 to 0.50 mm d⁻¹. These findings highlight that HR plays a critical role in sustaining soil moisture during extended dry periods and has a limited effect during precipitation events. Future research should investigate the effect of HR on other ecosystem processes, such as net ecosystem exchange of carbon and evapotranspiration under varying climatic conditions.

1 Introduction

Soil volumetric water content (VWC), defined as the amount of water stored in the unsaturated zone of the soil profile, is a fundamental state variable regulating ecosystem water and energy exchanges, particularly in dryland ecosystems (Seneviratne et al., 2010). Drylands cover over 40 % of Earth's terrestrial surface and support more than 38 % of the global population (Prăvălie, 2016), underscoring the importance of understanding soil moisture dynamics in these regions. While VWC provides a useful measure of soil water status, the movement and availability of this water are governed by soil water potential, particularly its matric component, which reflects the capillary and adsorptive forces binding water to soil particles (Hillel, 2003; Novick et al., 2022). In unsaturated soils, matric potential determines how tightly water is retained and how readily it can move toward plant roots. Because dryland ecosystem functioning is strongly constrained by precipitation variability (Beer et al., 2010; Ukkola et al., 2021), understanding how plants regulate water under fluctuating moisture conditions is essential for predicting ecosystem stability.

One key mechanism underlying this regulation is hydraulic redistribution (HR), the passive movement of water through plant roots, usually at night, from wet to dry regions of the plant rooting volume driven by differences in water potential. This passive process can favor plant survival during droughts by tapping into deep soil layers having relatively higher water potential and redistributing water to the shallow root zone (upward HR) (Nadezhdina et al., 2015; Prieto et al., 2012; Nicola and Ram, 2022). During wet seasons, HR can redistribute water from wet surface soil into deeper, drier soil (downward HR), supplementing the infiltration process in recharging deeper soil layers (Hultine et al., 2003; Scott et al., 2008; Fu et al., 2016; Bleby et al., 2010). Despite its potential role in regulating plant and ecosystem productivity, nutrient cycling and soil microbial activity (Grünzweig et al., 2022; Sardans and Peñuelas, 2014), most current dynamic global vegetation models and Earth system model still lack an explicit representation of HR (Fu et al., 2016).

HR has been observed across diverse ecosystems and plant species (Neumann and Cardon, 2012; Nadezhdina et al., 2010; Yu et al., 2013; Priyadarshini et al., 2016). It is recognized as a structural driver of dryland plant communities, regulating ecosystem productivity, and enhancing resilience to climate extremes (Lee et al., 2018; Barron-Gafford et al., 2021, 2017; Hafner et al., 2020). The dynamics of HR are influenced by various biotic (rooting architecture, plant capacitance, transpiration demand, senescence, and dormancy), abiotic factors (soil hydraulic characteristics, soil moisture status), and climatic conditions (precipitation) (Prieto et al., 2012; Katul and Siqueira, 2010; Wei et al., 2022). While several studies in arid and semi-arid ecosystems have reported upward HR during dry periods (Hao et al., 2013a; Lee et al., 2018; Scott et al., 2008; Yu et al., 2013) and downward HR following precipitation (Hultine et al., 2003), study on the fine-scale temporal variability of HR across multiple soil depths and multi-year timescales remains limited. Moreover, a quantitative understanding of how precipitation magnitude and frequency, key limiting factors in dryland ecosystems, influence HR rates is still lacking.

In this study, we explicitly test two hypotheses: (1) Direction of HR: Upward HR should be the dominant form of HR in dryland ecosystem. This is due to the recharge of deeper soil layers from precipitation which can retain moisture for longer periods, and during dry periods roots facilitate the movement of this retained water to the drier surface soils. (2) HR-precipitation relationship: Upward HR should decline following precipitation events, reaching its maximum rates during prolonged dry periods as drought creates steep water potential gradients between deeper, moist soil layers and the drier surface layers, facilitating the upward redistribution of water.

Meanwhile, soil moisture dynamics are governed by a complex interplay of forces that drive water movement through the soil profile. The primary drivers include matric potential (capillary and adsorptive forces binding water to soil particles), gravitational potential (driving downward drainage), and HR (Caldwell et al., 1998). These forces collectively determine soil water retention, redistribution, and plant water availability (Hillel, 2003). However, isolating their individual contributions from field soil moisture observations is challenging, because these processes operate simultaneously and are strongly influenced by soil properties, root activity, and atmospheric conditions. Consequently, a data-model fusion approach, which integrates process-based models with soil moisture observations, provides a robust framework to isolate and quantify HR, offering a more mechanistic and quantitative understanding of soil-plant water dynamics.

Several modeling studies have incorporated various HR schemes into process-based models to improve understanding of hydrological and ecological processes (Ryel et al., 2002; Amenu and Kumar, 2008; Wu et al., 2020; Fu et al., 2016; Zheng and Wang, 2007; Tang et al., 2015; Lee et al., 2018; Quijano and Kumar, 2015). However, realistic representation and estimation of parameters related to HR remains a challenge, as neither the magnitude of HR nor its associated parameters can be directly observed in the soil (Ryel et al., 2002; Quijano and Kumar, 2015). As a result, most models rely on default HR parameter values derived from Ryel et al. (2002). For example, in a study by Fabian et al. (2010), the maximum soil-root radial conductance (C_RT), a key parameter controlling HR, was assigned as the mean value between C_RT reported by Ryel et al. (2002) for Artemisia tridentata and by Williams et al. (1996) for Quercus-Acer stand. Similarly, Zheng and Wang (2007) and Yan and Dickinson (2014) prescribed a constant C_RT value based on Ryel et al. (2002). Alternatively, some studies estimated parameters during specific periods of time when upward or downward HR is assumed negligible, such as wet or dry season (Fu et al., 2018; Fu et al., 2016). These challenges in direct measurement, the reliance on assumed parameter values, and the parameterization of HR under the assumption of negligible redistribution constitute key gaps in our understanding of HR dynamics.

To address these gaps, we focused on piñon-juniper (PJ) woodlands, the most widespread semi-arid ecosystem in the US. PJ woodlands are spatially widespread, ecologically important, temporally dynamic, and structurally unique dryland ecosystem in the western US, spanning 10 US states and 40 million hectares across the American Southwest (Eastburn et al., 2024; Romme et al., 2009). Despite their importance, HR has not been previously studied in PJ woodlands. However, our continuous root sap flux measurements provided direct evidence of HR in both piñon and juniper roots, indicated by sustained negative root sap flux during nighttime at the study site (Fig. S1).

In this study, we used the process-based Terrestrial ECOsystem (TECO) model to (i) develop and implement a data assimilation approach to integrate HR into the TECO model; (ii) quantify and characterize the magnitude and dynamics of HR across multiple soil depths; and (iii) analyze the temporal patterns of HR and its relationship with precipitation magnitude and frequency. The TECO model is a well-established ecosystem model that integrates ecological processes to simulate carbon, water, and energy fluxes within terrestrial ecosystems (Weng and Luo, 2008). We employed data assimilation to constrain the TECO model including HR using four years of soil moisture data measured at multiple soil depths, encompassing both wet and dry periods.

2 Data and Methods

2.1 Study site and data

Our modeling study utilized data from a PJ woodland plot (Lat. 35.642, Long. −104.607, elevation 1925 m) located in New Mexico, USA, and previously described in Schwinning et al. (2020). The site is a private ranch covering an area of over 6800 ha that was ungrazed from 2012 through the measurement period used for this study and is characterized by a semi-arid climate. Mean annual precipitation of the site is approximately 460 mm, with the majority falling between May and October, and a mean annual temperature of 10.5 °C. The soil texture at the site varies with depth, ranging from loam to clay loam. The vegetation consists of distinct tree clusters dominated by piñon pine (Pinus edulis (Englem.)) and juniper (Juniperus monosperma (Englem.) Sarg.) separated by open areas of bare soil and herbaceous cover.

VWC was continuously monitored using multi-sensor frequency domain capacitance probes (Decagon EC-5) installed at four depths (5, 15, 30 and 60 cm), in four soil pits under the tree canopies. All sensors were monitored every minute by a datalogger (model CR6, Campbell Scientific), and 15 min averages were stored. For model parameterization, we used 15 min VWC records aggregated to daily means. Each sensor was calibrated in the lab before installation for both air and water frequency. Because soil temperature can affect both soil permittivity and the response of capacitance sensors, potentially confounding the small fluctuations in VWC caused by HR, temperature correction factors were applied to the measured VWC at each depth, using the nearest measured temperature, following the method described by Saito et al. (2009). Rather than excluding data below 0 °C, we used this temperature-correction approach to reduce the influence of temperature-driven artifacts on the soil moisture signal. This strategy allows retention of continuous soil moisture records while accounting for the known sensitivity of capacitance sensors to temperature-dependent changes in dielectric permittivity.

2.2 Modeling framework

TECO is a process-based ecosystem model (Hou et al., 2021; Jiang et al., 2018; Weng and Luo, 2008), and has evolved from the TCS model (Luo and Reynolds, 1999). The model consists of major components: canopy photosynthesis, plant growth, soil water dynamics, and soil carbon transfers. The canopy photosynthesis and soil water dynamics submodels run at the hourly time step whereas the plant growth and soil carbon submodels run at the daily time step. The model is driven by seven environmental variables, including precipitation (mm), wind speed (m s⁻¹), solar radiation (W m⁻²), air and soil temperature (°C), relative humidity (%), and vapor pressure deficit (kPa). The detailed description of TECO model is available (Weng and Luo, 2008) and only a brief description of soil water dynamics is provided here.

The soil profile is divided into 10 layers with a total depth of 180 cm, with varying thickness: 5 cm for the first layer, 10, 15, and 30 cm for the second, third, and fourth layers respectively, and 20 cm for each of the fifth through tenth layers. VWC in each layer results from the mass balance between influx and efflux, with changes primarily attributed to vertical unsaturated flow, transpiration, precipitation, runoff, and drainage. Evaporation depletes water from the first two soil layers, while transpiration depletes water from all soil layers containing roots, allocated based on root fraction in each layer (Eq. 8). Given the predominantly arid conditions of the study site, runoff and drainage were found negligible. Thus, water movement between soil layers is simulated as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}{W}_i}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{F}_i}{\mathrm{d}z}}-E_i-T_i\end{array}$$

where W_i is the water storage (cm) in layer i, t is time (h), F_i is net unsaturated flow of water into layer i (cm h⁻¹), z is vertical thickness, E_i and T_i are evaporation and transpiration water loss from layer i (cm h⁻¹).

The unsaturated soil water movement is simulated vertically according to modified form of Buckingham-Darcy's law (Campbell, 1985) (Eq. 2), with Brooks and Corey (1965) equation (Eq. 4) estimating hydraulic conductivity and soil water retention curve (SWRC) to simulate soil water potential (Ψ).

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{d}{F}_i}{\mathrm{d}z}}=K\left({\mathit{\theta }}_i\right)\left({\displaystyle \frac{\mathrm{d}{\mathrm{\Psi }}_i}{\mathrm{d}z}}+\mathrm{1}\right)\end{array}$$

where K(θ_i) is the unsaturated soil hydraulic conductivity (cm h⁻¹) for VWC θ (cm³ cm⁻³) in layer i, Ψ_i is soil water matric potential (MPa) in layer i, and z is the vertical thickness (cm) of the soil.

$$\begin{array}{}\text{(3)}& \mathrm{K}\left({\mathit{\theta }}_i\right)=K_{\mathrm{s}}{\left[{\displaystyle \frac{{\mathit{\theta }}_i-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}\right]}^{(\mathrm{2}m+\mathrm{3})}\end{array}$$

where, K_s is the soil saturated hydraulic conductivity (cm h⁻¹), m is the pore size distribution index, θ_s and θ_r are saturated and residual VWC (cm³ cm⁻³)

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}={\left({\displaystyle \frac{\mathrm{\Psi }}{{\mathrm{\Psi }}_b}}\right)}^{-\mathrm{1}/m}\end{array}$$

Ψ_b is the soil air entry water potential.

To quantify the direction and magnitude of HR, we integrated the HR model by Ryel et al. (2002) into equation 1 of TECO model (Eq. 5). This HR model empirically describes HR flux based on the soil water potential gradient between two soil layers (Eq. 6). HR was assumed to occur only at night, with its occurrence controlled by solar radiation instead of fixed day and night hours. Daytime starts as solar radiation exceeds 10 W m⁻², thereby inhibiting HR since the water potential gradient typically favors water movement from roots to canopy to meet transpiration demand during the day. This pattern is evident in Fig. S1, where under low or zero solar radiation, root sap flux was found to be negative, indicating water movement away from the root zone which is an indicator of occurrence of HR at the study site. Using these assumptions, the net water movement into soil layer i from other soil layers j can be expressed as:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}{W}_i}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{F}_i}{\mathrm{d}z}}-E_i-T_i+H_i\text{(6)}& {\displaystyle }{H}_i=C_{\mathrm{RT}}\sum _j\left({\mathrm{\Psi }}_j-{\mathrm{\Psi }}_i\right)max(c_i,c_j){\displaystyle \frac{R_i{R}_j}{\mathrm{1}-R_x}}{D}_{\mathrm{tran}}\text{(7)}& {\displaystyle }{c}_i={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\left(\frac{{\mathrm{\Psi }}_i}{{\mathrm{\Psi }}_{\mathrm{50}}}\right)}^b}}\text{(8)}& {\displaystyle }{R}_{d,i}={\displaystyle \frac{R_{\mathrm{0}}}{\mathrm{1}+{\left(\frac{d_i}{d_{\mathrm{50}}}\right)}^a}}\end{array}$$

where in Eq. (6), H_i is the net water redistributed by roots into layer i (cm h⁻¹) from other soil layers j, C_RT is the maximum radial soil-root conductance of the entire active root system for water (cm MPa⁻¹ h⁻¹), Ψis soil matric potential (MPa), c_i is a factor reducing soil-root conductance based on Ψ_i, R_i is the fraction of active roots in layer i, R₀ is the average vertically summed root dry mass from the bottom to the root zone to the soil surface, and D_tran is a factor reducing water movement among layers by roots while plant is transpiring and is assumed to be 1 during the night when transpiration is minimal and 0 during day. R_x=R_i when θ_i>θ_j or R_x=R_j when θ_j>θ_i. In Eq. (7), Ψ₅₀ is the soil water potential (MPa) where conductance is reduced by 50 % and b is an empirical constant. In Eq. (8), R_d,i is cumulative fraction of roots above soil depth d_i (cm) for the ith layer, and d₅₀ is the soil depth (cm) at the median of the root distribution and a is the shape parameter (Table 1). The Brooks and Corey (1965) model for SWRC was utilized to simulate soil water potential (Ψ), facilitating the development of soil water potential gradients necessary for HR by tree roots (Eq. 4). Due to lack of site-specific parameters, the default values of b and Ψ₅₀ were used as 3.22 and −1 MPa, respectively (Ryel et al., 2002).

2.3 Data assimilation for parameters estimation

We used Bayesian probabilistic inversion to calibrate parameters associated with soil hydraulics, where posterior probability density functions of parameters are obtained from prior knowledge about the parameters and the error between model and observations. According to ..Mosegaard and Sambridge (2002), Bayesian inversion can be summarized by the following equation:

$$\begin{array}{}\text{(9)}& p\left(c|Z\right)\propto p\left(Z|c\right)p\left(c\right)\end{array}$$

where p(c|Z) is posterior probability density function of model parameters c; p(Z|c) is a likelihood function of parameters c; p(c) is prior probability density function of parameters c. We assumed that the prediction errors were normally distributed and uncorrelated, hence, the likelihood function, p(Z|c), was calculated as follows:

$$\begin{array}{}\text{(10)}& p\left(Z|c\right)\propto \exp \left\{-\sum _{i=\mathrm{1}}^k{\displaystyle \frac{(Z_i-X_i{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_i^{\mathrm{2}}}}\right\}\end{array}$$

where Z_i is observed VWC at ith soil layer, X_i is VWC simulated by TECO at a corresponding soil depth; ${\mathit{\sigma }}_i^{\mathrm{2}}$ is the variance of a measurement at a soil layer; k is the total number of soil layers.

To generate the posterior distributions, we first specified the priors of the parameters to be uniformly distributed over the intervals specified in Table 1. We put constraints on parameters based on literature. The initial set of parameters was randomly selected within the prior parameter ranges. Once we specified parameter ranges, we used the Metropolis–Hastings (M–H) algorithm (Hastings, 1970; Metropolis et al., 1953), a Markov chain Monte Carlo method, to sample from the posterior parameter distribution. To generate a parameter set, we ran M–H algorithm in two steps: proposing step and moving step. In the proposing step, a new parameter set c^new was generated from a previously accepted parameter set c^k−1 through a proposal distribution ($c^{\mathrm{new}}|c^{k-\mathrm{1}}$):

$$\begin{array}{}\text{(11)}& c^{\mathrm{new}}=c^{k-\mathrm{1}}+r\times {\displaystyle \frac{c^{\max }-c^{\min }}{D}}\end{array}$$

The value of $P\left(c^{k-\mathrm{1}}\right|c^{\mathrm{new}}$) was then compared with a random number U from 0 to 1. Parameter set c^new was accepted if $P\left(c^{k-\mathrm{1}}\right|c^{\mathrm{new}})\ge U$, otherwise c^k was set to c^k−1. In the moving step, a probability of acceptance $P\left(c^{k-\mathrm{1}}\right|c^{\mathrm{new}}$) was calculated as in the following (Marshall et al., 2004):

$$\begin{array}{}\text{(12)}& P^{\left(c^{k-\mathrm{1}}|c^{\mathrm{new}}\right)}=min\left\{\mathrm{1},{\displaystyle \frac{p\left(Z\right|c^{\mathrm{new}}\left)p\right(c^{\mathrm{new}})}{p\left(Z\right|c^{k-\mathrm{1}}\left)p\right(c^{k-\mathrm{1}})}}\right\}\end{array}$$

The M–H algorithm was repeated for 50 000 simulations, and then all accepted parameters values were used to generate the probability distribution functions (Xu et al., 2006).

To evaluate the impact of HR on soil moisture dynamics in a PJ woodland, we conducted two multi-year simulations using two configurations of the TECO model: TECO+HR (with HR) and default TECO (HR turned off). To distinguish the influence of HR from soil hydraulic properties, we adopted a data assimilation approach focused on calibrating only the TECO+HR model. Soil moisture observations were available at 5, 15, 30 and 60 cm depths, and data assimilation was therefore applied only to these four soil layers over a four-year period, and data-model comparisons are presented exclusively for these depths. The parameters values deeper than 60 cm were not calibrated and were empirically defined from literature values. This modeling strategy is consistent with other data assimilation studies, in which model calibration is restricted to soil layers with available observations (Fu et al., 2016; Zhang et al., 2016). In total, 21 parameters were optimized including 20 soil hydraulic parameters (θ_s, θ_r, K_s, m, and Ψ_b for each of the four soil depths) across four soil layers and one HR-related parameter (C_RT). The prior range of soil hydraulic parameters were informed by established relationships between soil texture and hydraulic properties (Rawls et al., 1982; Clapp and Hornberger, 1978). The prior range for C_RT was based on values reported in Fu et al. (2016). Within this range, we optimized depth-specific soil hydraulic parameters to achieve a close match between modeled and observed soil moisture (Table 1).

After calibrating the TECO+HR model, we deactivated the HR process and ran simulations with the same optimized parameters to generate the default TECO scenario. Before each model simulation, we performed a 200-year spin-up separately for both model configuration (TECO+HR with HR active and default TECO with HR disabled) to ensure that each model reached stable carbon stocks as initial conditions.

The motivation to calibrate only TECO+HR model, rather than the default TECO is to avoid parameter compensation for unresolved processes (Luo and Schuur, 2020), in which the absence of HR could lead to unrealist adjustments of soil hydraulic parameters to indirectly capture its effects. This approach allowed us to ensure that differences in soil moisture dynamics between TECO+HR and default TECO simulations were attributable solely to the presence or absence of HR.

Furthermore, to evaluate the influence of cumulative precipitation and soil moisture memory on HR, we calculated the Antecedent Precipitation Index (API) for the study period (2018–2021) following Kohler and Linsley (1951). API acts as a proxy for soil moisture status by accounting for the decaying effect of past rainfall events. The daily API (API_t) was calculated using the recursive decay function:

$$\begin{array}{}\text{(13)}& {\mathrm{API}}_t=P_t+(k\times {\mathrm{API}}_{t-\mathrm{1}})\end{array}$$

where P_t is the precipitation on day t (mm), API_t−1 is the index value of the preceding day, and k is a decay constant representing the recession of soil moisture due to evapotranspiration and drainage. We used a decay constant of k=0.90, which falls within the commonly applied range for antecedent precipitation indices and is consistent with optimization analyses indicating optimal decay constants near 0.90 (Li et al., 2021). This metric enables differentiation between short dry intervals following wet conditions and extended dry spells with limited antecedent moisture.

2.4 Statistical analysis

Model performance was assessed by comparing simulated outputs with observed data during full simulation periods (2018–2021), dry, and wet periods, defined as days without and with rainfall events, respectively. During the study period, wet days accounted for 22 % of all days, whereas dry days comprised the remaining 78 %. Model performance was evaluated using two statistical metrics: root mean square error (RMSE) and absolute mean error (MAE). RMSE is widely used to quantify model accuracy, but its squared-error formulation overemphasizes the effects of large deviations (Willmott and Matsuura, 2005). Therefore, MAE was also calculated as a measure of the average magnitude of deviation from observed values (Alfieri et al., 2017). Both metrics were calculated as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n(m_i-o_i{)}^{\mathrm{2}}}\text{(15)}& {\displaystyle }\mathrm{MAE}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n{\left|m_i-o_i\right|}^{\mathrm{2}}\end{array}$$

Where: o_i represents observed values, m_i represents modeled values, and n represents the number of data points.

Table 1 Parameters constrained using data assimilation in the TECO model from soil moisture data from 2018 to 2021.

Note: Four values represent parameters in the four modeled VWC at the depths of 5, 15, 30, and 60 cm, respectively.

3 Results

3.1 Parameter estimation via data assimilation and water mass balance

The data assimilation approach, using VWC data to constrain the model, yielded well-constrained soil hydraulic parameters (Table 1; Figs. S2 and S3). The resulting posterior probability density functions, characterized by sharp peaks, narrow spread, and consistency across soil depth support the reliability and accuracy of these calibrated parameter values. Additionally, soil water mass balance of soil profile was conserved before and after incorporating the HR process into the TECO model (Fig. S4). The key components of the water budget: precipitation, evapotranspiration, and changes in soil water content remained balanced, ensuring that the model accounted for all water fluxes. Furthermore, the sum of HR across all soil layers (10 layers) was consistently equal to zero, further ensuring that no water was artificially introduced or lost from the system.

3.2 Observed and simulated soil moisture

The data assimilation-constrained models generally captured both the magnitude and dynamics of observational data, reproducing seasonal variations in soil moisture across four soil depths. Minor mismatches at the topsoil (5 cm) likely reflect the complexity of near-surface processes that are not fully represented in simplified models and potential sensor lag due to imperfect soil contact, as similar discrepancies are not observed at deeper layers. In addition, following prolonged dry periods, initial rainfall can be largely offset by evaporation, resulting in muted surface soil moisture responses (Miele et al., 2023; Cattry et al., 2025; Asadollahi et al., 2022). While TECO+HR simulation showed an improvement in the overall model performance, the impact of HR was mostly pronounced during dry periods (Figs. 1 and 2). We further examined diurnal soil moisture fluctuations (Fig. S5) and found that TECO+HR closely tracked the observed diurnal cycles, whereas the default TECO failed to capture this pattern, suggesting that the observed diurnal variability was likely driven by HR. Additionally, we compared min–max normalized soil matric potential at 15, 30, and 60 cm with simulations derived from Eq. (4) (Fig. S6). Both models reproduced the general trends of the observations, suggesting that the simulated soil water potential gradients were consistent with measurement.

Moreover, during periods of limited precipitation, the TECO+HR (blue lines) consistently maintained higher soil moisture compared to default TECO (red lines), aligning closer to observation particularly in the topsoil layers (Fig. 1a–c). Following precipitation events, the default TECO and TECO+HR simulations converged, suggesting the minimal influence of HR under wet conditions at the study site. However, as surface soil moisture decreased following precipitation, the two simulations diverged again, with TECO+HR maintaining higher moisture levels in the topsoil layers, highlighting the role of HR in maintaining soil moisture during prolonged drought.

Figure 1 (a–d) Observed and simulated soil volumetric water content (from January 2018 to December 2021) at soil depths of 5 cm (a), 15 cm (b), 30 cm (c), and 60 cm (d). Black lines indicate observations, cyan lines indicate TECO+HR, and magenta lines indicate default TECO. Vertical blue bars indicate daily precipitation (right axis).

The incorporation of HR into TECO resulted in reductions in model errors. During dry periods, the RMSE decreased by 25 %, 43 %, and 52 % at 5, 15, and 30 cm soil depths, respectively. However, limited improvement was observed at 60 cm soil depth. Correspondingly, the MAE was reduced by 30 %, 53 %, and 60 % at 5, 15, and 30 cm, respectively. Over the entire study period, RMSE decreased by 24 %, 25 %, and 47 % at 5, 15, and 30 cm, with MAE reductions were 29 %, 34 %, and 55 % at the same depths (Fig. 2a–d). Overall soil profile performance improved as well, with RMSE and MAE reductions over 40 % for both the four-year simulation and dry periods. These improvements during dry periods are especially important, as roots are most vulnerable to drought. By mitigating soil water deficits in surface layers, HR could reduce the risk of hydraulic failure, thereby supporting plant species survival and it could enable better prediction of ecosystem responses to water stress, such as carbon uptake (Domec et al., 2010), and evapotranspiration (Zhu et al., 2017). In contrast, during wet periods, HR had minimal influence on soil moisture (Fig. 2e, f).

Figure 2 Model performance for soil moisture across different depths (5, 15, 30, and 60 cm, and 0–60 cm integrated soil profile), considering temporal variations in soil moisture conditions. Root Mean Square Error (RMSE, %) and Mean Absolute Error (MAE, %) are presented for the complete time series (a, b), dry periods (c, d), and wet periods (e, f). Lower values of RMSE and MAE indicate better model performance.

3.3 Effect of HR on soil moisture

The direct impact of HR on hydrological processes should be evident in the soil profile water content. We tested this by comparing VWC model simulations in TECO with and without HR processes, to observed VWC time series at four depths (Fig. 3). We found that cumulative effects of HR on soil moisture vary with depth, primarily due to the non-uniform root biomass distribution throughout the soil profile (Fig. S8). The most pronounced effects of HR were observed in the topsoil layers (5, 15, and 30 cm), where average daily water content increased by up to 60 % compared to simulation without HR. This increase was driven by upward HR, especially during dry-down periods (Fig. 4b).

Figure 3 Relative change in modeled volumetric water content (VWC, %) across four soil depths (5, 15, 30, and 60 cm). The color gradient represents the relative change (%) in modeled VWC, calculated as (HR − No HR)$/$No HR × 100, with HR and No HR indicating simulations with and without hydraulic redistribution, respectively. The relative change is overlaid on observed VWC.

Figure 4 Temporal dynamics of hydraulic redistribution (HR). (a) diurnal pattern of modeled HR across soil depths from 22–24 July 2018. The graph illustrates HR patterns during a dry period followed by a precipitation event (right y-axis). Colored lines represent different soil depths. (b) long-term daily HR trends and precipitation from January 2018 to December 2021. The blue vertical bars represent precipitation (right y-axis).

3.4 HR simulations

Model simulation revealed distinct patterns of HR dynamics across soil depths and temporal scales. Figure 4 illustrates these patterns over two timescales: a short-term, diurnal pattern (Fig. 4a), and a long-term perspective from 2018 to 2021 (Fig. 4b). HR is a process with both a source and a sink for water movement. In Fig. 4, positive HR suggests that a soil layer is gaining water (sink), whereas negative HR values suggest that the layer is losing water (source).

The short-term modeling analysis highlights diurnal pattern of HR during dry conditions and a precipitation event (Fig. 4a). For instance, on 23 July 2018, during a dry period, upward HR occurred, moving water from deeper (>100 cm) to shallower (0–30 cm) soil layers. However, following a precipitation event on 24 July 2018 (12 mm), this pattern shifted. The top 5 and 15 cm layers showed negative HR and a decrease in the upward HR rate, respectively, acting as a water source for deeper layers. At the same time, deeper soil layers showed a decline in negative HR rates, suggesting signs of receiving water likely from the topsoil layers. The sum of HR across all soil layers remained zero, confirming that HR redistributed water rather than adding to the system. Consequently, downward HR from the topsoil supplemented infiltration, enhancing water movement into deeper soil layers, reflected by a decrease in the negative HR rates at depths and an increase in the positive HR rate at 30 cm (Fig. 4a).

While our model simulates HR across 10 soil layers, we present long-term results for only the top four soil layers (5, 15, 30, and 60 cm) to enable direct comparison with the available observed soil moisture data. A clear seasonal pattern emerged, with HR generally intensifying during dry periods (Fig. 4b).

Our model showed that upward HR was predominantly occurring in up to top 30 cm of soil profile, with values ranging from −0.066 to 0.29 mm d⁻¹ in each soil layer and an average of 0.30 mm d⁻¹ across the top 30 throughout the study period. Downward HR, while less pronounced, moved water only from the 5 cm soil layer during monsoon seasons and large precipitation events (e.g., July 2018, 2019, 2020, and 2021; Fig. 4b). In contrast, 60 cm soil layer typically exhibits a negative HR during dry periods, acting as a water source for upper layers, and positive HR during wet periods, suggesting occasional water input from surface layers ranging from −0.096 to 0.059 mm d⁻¹ (mean 0.0015 mm d⁻¹). Moreover, integrated soil profile (top 60 cm of soil profile), showed that upward HR was the dominant form of HR throughout the year, ranging from 0.10 to 0.53 mm d⁻¹ with a mean value 0.31 mm d⁻¹ (Fig. S7).

3.5 Precipitation influences on HR

The model results showed a significant linear relationship between weekly HR and mean weekly API (mm) (Fig. 5a–b). In the topsoil layers (5 and 15 cm), HR exhibited negative relationships with API (R²=0.31 and 0.30 respectively, both p-values < 0.001), indicating that HR activity decreased as antecedent moisture conditions became wetter. In contrast, positive correlation was observed at 30 and 60 cm depths (R²=0.13 and 0.14 respectively, both p-values < 0.001), suggesting downward redistribution under wetter antecedent conditions.

Figure 5 Relationships between hydraulic redistribution (HR) and Antecedent Precipitation Index (API). (a) Weekly HR rates versus mean weekly API at different soil depths (5, 15, 30, and 60 cm). For each soil depth, the trend lines, R², and corresponding p-values are shown. (b) Depth-integrated weekly HR across 0–60 cm soil profile versus mean weekly API, with trend line, R², and p-value.

Additionally, when HR was integrated across the soil profile (0–60 cm), a significant negative relationship with mean weekly API was observed (R²=0.45, p-values < 0.001, Fig. 5b). This suggests that overall HR activity was greatest under drier antecedent conditions and declined as cumulative moisture availability increased, highlighting the stronger role of HR in regulating soil water dynamics during prolonged dry periods.

Figure 6 Relationship between weekly mean hydraulic redistribution (mm d⁻¹), dry spell length between two rainfall events (log₁₀(days + 1)). The color scale indicates the number of rainfall events per week, while marker size represents the weekly precipitation amount (mm per week). Dry spell length denotes the number of rainless days between two consecutive precipitation events. The x-axis shows dry spell length transformed as log₁₀(days + 1) to allow inclusion of zero-length dry spells.

Despite these clear trends, variability in HR was observed across the full range of mean weekly API values (Fig. 5a–b). This variability could be attributed to the rainfall frequency, event size, and the duration of dry periods between rainfall events (Fig. 6), factors that are implicitly captured by API. Figure 6 provides a more detailed view of these controls by relating HR directly to dry spell length and rainfall frequency, suggesting that HR may increase not only with reduced precipitation frequency but also as the interval between consecutive rainfall events lengths. HR. HR was lowest under conditions of high rainfall frequency and shorter dry spells, progressively increasing to its peak in the absence of rainfall. However, as the drought period extended beyond 30 d, HR declined, suggesting potential limitation on availability of deeper water to sustain HR. This variability is further illustrated through three scenarios (Fig. S7): (1) Following a rainfall event (28 mm on 5 July 2018), HR in the top 60 cm of soil profile was minimal at 0.13 mm d⁻¹, indicating limited driving force for water redistribution when soil moisture was abundant. (2) During a transition period between rainfall events (5–10 July 2018), HR gradually increased but remained moderate, ranging from 0.13 to 0.20 mm d⁻¹, suggesting a progressive activation of the redistribution process as soil began to dry. (3) During a prolonged dry period (23–30 November 2018), HR peaked at 0.20–0.52 mm d⁻¹, demonstrating enhanced redistribution activity in response to the development of soil moisture gradients.

4 Discussion

4.1 Patterns of hydraulic redistribution

Our findings support the hypothesis that upward HR is the dominant form of HR in dryland ecosystems due to limited precipitation amount and sporadic rainfall events (Fig. S7). This prevalence of upward water movement is characteristic of semi-arid regions, where deep-rooted plants often redistribute water from moist deeper layers to drier surface soils during periods of water stress (Caldwell et al., 1998; Ryel et al., 2002). Notably, the most pronounced HR occurred in the topsoil layer (5, 15, and 30 cm), (Fig. 4b), which can be attributed to vertical root distribution, with over 50 % of root biomass concentrated in the top 30 cm (D₅₀=25 cm) of the soil profile (Fig. S8). Similar relationships between root distribution and HR intensity have been reported in previous studies, where deeper root systems extend redistribution to deeper soil layers, whereas shallow root systems amplify HR effects in surface soil due to higher root density and activity (Hao et al., 2013b).

The magnitude of HR simulated in this study (0.10–0.53 mm d⁻¹ for top 60 cm soil) falls within the range reported in previous studies. Estimated HR rates for the topsoil were comparable to values synthesized in the global review by Neumann and Cardon (2012), which reported HR magnitudes ranging from 0.04 to 3.2 mm d⁻¹ across ecosystem. However, our estimates are slightly higher than the upper range reported by Yang et al. (2022) for desert or sparsely vegetated ecosystems (0.014–0.475 mm d⁻¹). These differences likely reflect variations in vegetation structure, rooting depth, soil texture, and water availability among ecosystems.

4.2 Effects of precipitation variability on HR

Our findings support the hypothesis that the precipitation pattern significantly influences the magnitude and variability of HR (Figs. 4a–b, and S7). The rate of HR in the topsoil profile (<60 cm) exhibited a consistent response to precipitation events, characterized by sharp declines following large rainfall, and a gradual recovery to pre-rain levels during subsequent dry periods. Similar responses have been reported in previous study, where precipitation temporarily suppresses HR by reducing water potential gradients between shallow and deeper soil layers (Hao et al., 2013b). However, as water redistributes through the soil profile, new hydraulic gradients re-establish, leading to enhanced HR activity. In this phase, roots actively redistribute water from newly moistened deep layers to drier shallow layers (Yu and D'Odorico, 2014; Ryel et al., 2002).

Our model predicted that HR rates were generally higher during rainless periods compared to rainfall periods within a given year. For instance, during the prolonged dry period in 2020 (driest year), HR rates remained consistently high, 0.17–0.40 mm d⁻¹, with minimal fluctuations. The consistent high HR rates, likely arises from more pronounced soil water potential gradients derived from sustained plant water demand and surface evaporation in the absence of frequent precipitation (Fu et al., 2016; Meinzer et al., 2004).

Seasonal patterns include higher HR rates (0.12–0.53 mm d⁻¹) during the drier periods (typically from November to May) and lower rates (0.10–0.30 mm d⁻¹) during the monsoon season (usually from June to October) (Fig. S7). This seasonality underscores the influence of both precipitation patterns and potential evapotranspiration on HR dynamics, highlighting that HR is likely more pronounced during drier seasons when soil moisture gradients are likely to be more substantial due to reduced precipitation and potentially higher evaporative demand (Scott et al., 2008; Fu et al., 2016; Yu and D'Odorico, 2014).

4.3 Limitation and future perspectives

While our modeling study provides valuable insights into HR dynamics in PJ woodlands, several limitations should be noted. (1) The model does not account for interannual changes in vegetation cover or species composition. Variations in plant functional types and leaf area index may influence soil moisture and HR, and incorporating these dynamics could improve long-term simulations. (2) Our analysis focused on the dominant tree species at study site; however, other plant species may also benefit from water redistributed by these trees, potentially influencing ecosystem water dynamics. (3) We did not include stem water refilling or nighttime transpiration reported by Howard et al. (2009) and Neumann et al. (2014), which may influence the magnitude of HR. (4) Finally, future studies should address these limitations and further investigate the role of HR in regulating ecosystem functions, such as carbon exchange and evapotranspiration.

5 Conclusions

This study demonstrates the role of hydraulic redistribution (HR) in soil water dynamics in piñon-juniper woodlands. By integrating HR processes and observations into the Terrestrial Ecosystem Model (TECO) via data assimilation, we successfully constrained model soil hydraulics parameters and improved simulations of soil water content across multiple depths, particularly in shallow soil layers (0–30 cm) and during dry periods. Our model results indicate that HR rates vary in response to the duration of dry spells between rainfall events. Generally, HR rates tend to increase as soil becomes drier and decreases with increasing precipitation magnitude and frequency. Across the wet to dry transition, HR rates exhibit a range of 0.10 to 0.50 mm d⁻¹. Consequently, HR increased soil moisture in topsoil layers by up to 60 % during dry periods, with upward HR emerging as the dominant flux, especially in the top 30 cm. These findings underscore the potential influence of HR during dry periods and highlight its role in sustaining soil water availability for vegetation. Future research should explore how HR-mediated water redistribution affects ecosystem functions including carbon exchange, and evapotranspiration.