The study site is located in the southern Sacramento Mountains near the village of Cloudcroft, NM. Elevations are 2,127 m in the valley site and 2,748 m on the ridge (Figures 1A–C). The climate is considered warm-dry (Fulé et al., 2009), with a mean monthly temperature of 22 °C in July and 5 °C in January, and a mean annual precipitation of 752 mm (1987–2016; WRCC). Topography is steep, with most slopes ranging from 20 to 40%, and shallow soils approximately 35 cm deep, consisting of clay and clay loam textures (Garduño et al., 2015). The average wilting point (Wp), field capacity (Fc), and hydraulic conductivity (Kh) are shown in Table 1. The soil type is Typic Argiborolls-Aquic Haploborolls. The experimental watershed lies above the Permian San Andres and upper Yeso formations. A complex hydrological system is influenced by exposed geological features and regional fracture systems. This system enables rapid movement of much of the high-altitude precipitation through interconnected shallow carbonate aquifers perched and linked by fracture networks and surface water drainages (Newton et al., 2015).

The overstory community is dominated by Douglas fir (Pseudotsuga menziesii Mirb.) and ponderosa pine (Pinus ponderosa C. Lawson). The understory herbaceous community is dominated by Stipa and Sporobolus species, Alder leaf mountain-mahogany (Cercocarpus montanus Raf.), and Longflower snowberry (Symphoricarpos longiflorus A. Gray) (Garduño et al., 2015). The height of trees, diameter at breast height (DBH), and density ranged from 6 to 13 m, 10 to 25 cm, and 80 to 159 trees ha⁻¹, respectively (Table 2).

In mid-August 2009, one randomly selected plot was thinned from each of four experimental plot pairs (8,100 m² or a dimension of 90 m × 90 m each plot) (Figure 2). The thinning was performed using a boom-mounted masticator attached to a Caterpillar 320D LLR Hydraulic excavator (23,700 kg weight), which shredded trees from the top to the base. The resulting slash (i.e., bark and branches) was left on site (Supplementary Figure 1A). The other plot within a pair was left as an untreated control plot. Within plot pairs (i.e., blocks), the treatment and the control plot had similar aspects (north or south). Two plot pairs were installed on the ridge and two on the valley (Supplementary Figure 1A). This blocking design with four replicates (two in the ridge and two in the valley) is robust for detecting treatment effects. However, the limited number of blocks may constrain the spatial extrapolation of the findings to landscapes with greater topographic complexity.

On control plots, pine needles and litter covered the soil, ranging from 46.81 to 95.29% with a depth of 1.27 to 4.43 cm. Pre-treatment canopy cover (i.e., calculated with the Normalized Difference Vegetation Index) on selected thinned plots averaged 85 and 77 % in control plots. Post-treatment canopy cover averaged 30% with slash depth ranging from 1 to 7 cm (Garduño et al., 2015).

The soil water retention curve was used to calculate the soil volumetric water content at field capacity, θ_fc, and permanent wilting point, θ_pwp (Supplementary Table 1). Mean daily values of ψ and θ were used to generate SWRC for each depth and site. Soil ψ and θ were screened for outlying values, which were removed from the subsequent curve analyses. After the screening process, ψ was converted to cm of H₂O by multiplying MPa by 10197. Then, the soil particle size distribution data, ψ and θ were used as input data to calculate the SWRC for each plot at each depth. The Retention Curve program (RETC, Version 6.02) was used to fit the SWRC using the van Genuchten model described by

Where s is the effective degree of saturation, α, k (equivalent to the van Genuchten shape parameter n), and m are empirical constants affecting the shape of the retention curve (Supplementary Table 1). For further information, (see van Genuchten, 1980). Equation 1 was selected for fitting the SWRC because it performs best for many soils (van Genuchten et al., 1991). Unsaturated hydraulic conductivity (K) was calculated for each depth using Gardner's equation (Gardner, 1958) described by

Where n is the porosity, a and b are constants estimated by hydraulic properties (Campbell, 2008) described by

Where lnθ_pwp is the natural logarithm of the wilting point, lnθ_fc is the natural logarithm at field capacity of the soil volumetric water content (cm³ cm⁻³) for each soil at each depth determined with the RETC program. Although hysteresis can significantly influence water flux, it was not considered in describing water flux, as no unique functional relationship can be easily assumed (Weiler and McDonnell, 2004).

Soil temperature, soil volumetric water content, and water flux were analyzed by depth using a mixed model with fixed effects for treatment, month, and their interaction, but focused on two-tailed pre-planned comparisons of control-treatment differences within specific seasons. To allow considering a relatively flexible set of parsimonious candidate covariance structures, we fitted covariance structures with random pair × month effects, and different combinations of random plot effects (uniquely identified as pair × treatment) and AR (1) (autoregressive) variance structures fit pooling across treatment or separately by treatment. We based formal inference on the covariance structure with the lowest AICc (Akaike Information Criteria). We conducted 12 pre-planned comparisons of treatment to control by averaging across the months in the indicated seasons. Seasons fell into the following five categories, with all but the pre-treatment season corresponding to multiple pre-planned comparisons for each incidence of that season over the 2009-2011 study years. (1) The pre-treatment season (May-August 2009); (2) The dry season fall (October) and spring (May); (3) The snow season (November–January); (4) The snowmelt season (March–April); and (5) The rainy season (June–September). We assessed the sensitivity of the findings to the specific model chosen and to outliers by fitting the two models with the lowest two AICc values with and without extreme outliers, defined as those observations with marginal studentized residual values of magnitudes exceeding 4. Reported results are based on the specific model chosen, noting when findings were sensitive to either the alternative model or outliers. We defined significance as α ≤ 0.05 and conducted analyses using SAS (version 9.4−2006 ^® SAS Institute Inc).

Based on the Cooperative Climatological Data Summaries of precipitation data covering 30 years, the highest precipitation was in July, followed by August, while the driest month was April, followed by May (Table 3). For 2009, 2010, and 2011, the total annual precipitation recorded was 238, 408, and 159 mm, respectively. The lowest monthly precipitation of 1 mm was recorded in November 2010, March, and April 2011, while the highest monthly precipitation of 140 mm occurred in July 2010 (Table 3). The rainy season was marked by a different precipitation regime throughout the study period (2009–2011). The lowest precipitation, during the rainy season, was recorded in 2011 (140 mm), while the highest precipitation of 299 mm was recorded in 2010; 2009 was in the middle of this precipitation regime with 169 mm. Examining the seasonal distribution of precipitation (1987–2016), it was observed that the study period (2009–2011) coincided with a drought phase. This period was characterized by a reduction in monthly precipitation amount, evident in the overall seasonal distribution of precipitation.

The variability of slash cover temperature and air temperature followed a seasonal trend characterized by a gradual rise and decline according to the seasons. The highest slash cover temperature and air temperature were recorded in June 2011 (22 °C and 21 °C, respectively). The lowest slash cover temperature and air temperature below zero were obtained from December 2009 through February 2010, except for January, with an air temperature above zero (Figure 4).

The soil temperature at 7 cm depth was higher in the thinned plots as compared to the control plots. The difference was statistically significant from May 2010 (dry season) through the rainy season of 2010, with differences of 2.65 ± 1.07 °C and 2.15 ± 0.96 °C, respectively (Figure 5). At 20 cm depth, temperatures in the thinned plots were higher than those observed for the control plots from snowmelt until October (fall dry season). The difference was statistically significant in May 2011 (dry season), rainy 2011, and October 2011 (fall dry season), with the estimated difference between treated and control plots being 2.50 ± 1.09 °C, 2.62 ± 0.99 °C, and 2.49 ± 1.09 °C, respectively (Figure 5). At 35 cm depth, no significant differences were observed between thinned and control plots (Figure 5).

At the 7 cm depth, θ was similar between the thinned and control plots (Figure 6). Both treatments showed a sharp decrease in θ approaching to wilting point of 0.07 cm³ cm⁻³ (thinned plots) and 0.11 cm³ cm⁻³ (control plots). On the other hand, thinned plots approached field capacity (0.32 cm³ cm⁻³) during the snowmelt season in March 2010.

At 20 cm depth, the thinned plots had a 0.12 ± 0.05 cm³ cm⁻³ higher than the control plots in October 2009 and 2010 (fall dry season). The significance of this difference was model-dependent (Figure 6). In control plots, the field capacity (0.26 cm³ cm⁻³) was attained in February and March 2010. In contrast, the low precipitation registered from October 2010 through October 2011 led the soil to the wilting point (0.10 cm³ cm⁻³) (Figure 6). On the other hand, thinned plots had a θ in the range of 0.20 to 0.30 cm³ cm⁻³. Although the wilting point and field capacity were at 0.13 and 0.45 cm³ cm⁻³, respectively, the θ was at the midpoint for the study period.

At 35 cm depth, the thinned plots had a consistently higher θ than the control plots, starting from May 2010 (spring dry season) through October 2011 (fall dry season) (Figure 6). The snowmelt of 2010 led to similar soil volumetric water content in thinned plots, throughout the soil profile (i.e., three different depths). After the dry season of 2010 (May), the precipitation of 140 mm (July) contributed to replenishing the soil (i.e., close to the field capacity of 0.41 cm³ cm⁻³), followed by a steady decrease until October 2011. Conversely, control plots had a low soil volumetric water content at 35 cm depth that followed a narrow range through the study period, likely attributed to tree root uptake.

In control plots, the total hydraulic head data showed depth-specific patterns (Figure 7A). During the study period, at all depths, the total head increased as a response to low precipitation events. At 7 cm depth, the total head was the most influenced by precipitation and dry spells, even reaching −1,200 cm in June 2011, while at 20 and 35 cm depths reached −700 and −800 cm, respectively.

In thinned plots, water flowed from the deeper layer (35 cm) to the shallower layer (7 cm) (Figure 7B). After the thinning treatment, all depths diminished in total head until September 2010. After this period, it increased (more negative), although at 35 cm depth, there was a less pronounced increase. Thus, at 35 cm depth, the total hydraulic head was in a range of −200 to −590 cm even with prolonged low precipitation (February and August 2011). The response and direction of water flow occurred after the thinning treatment, mainly due to the lack of tree influence, which likely diminished water uptake at deeper soil layers, thereby altering the total hydraulic head.

In control plots, during the study period, the hydraulic gradient was mainly positive at 7–20 cm depth with a sharp increase from January 2010 through December 2010. At this depth, the hydraulic gradient showed multiple peaks and fluctuations, mainly after precipitation (Figure 8A). At 20–35 cm depth, the hydraulic gradient was mostly negative. At this depth, the hydraulic gradient remained steady around zero, indicating minimal variation. The opposite direction at both depths was mainly due to water uptake by trees (probes were installed between three pines). The total hydraulic head corroborates the direction of water flow shown on the hydraulic gradients during low precipitation (April through July 2011).

In thinned plots, after the thinning treatment, the hydraulic gradient was positive at both depths (7–20 cm and 20–35 cm) even during low precipitation (February through June 2010 and from October 2010 to the first week of July 2011) (Figure 8B). At 7–20 cm depth, the hydraulic gradient increased steadily from October 2009 through May 2010. Henceforward, the hydraulic gradient had a sharp decrease until August 2011. At 20–35 cm depth, it remained relatively stable with small variation. During the described period, at both depths, water flow was from the high total head to the low total head as confirmed by the hydraulic head results (Figure 7B).

Following thinning, at 7 to 20 cm depth, seasonal estimates of soil water flux in control plots ranged from −0.002 (dry season 2009) to 0.014 mm day⁻¹ (snowmelt 2011), while for thinned plots, the mean q ranged from−0.060 (pre-treatment) to 0.051 mm day⁻¹ (dry season 2011). Both plots showed similar trends in q with upward (+) and downward (–) directions, although the control plots showed slight upward and downward q movement (Figure 9A). The pre-treatment period showed a downward direction for both plots (i.e., control and thinned); none of these differences were significant (Table 4).

Following thinning, at 20 to 35 cm depth, control plots differed significantly from thinned plots beginning with the snow season 2009–2010 (November 2009–January 2010) through the rainy season 2010 (June–September 2010; Figure 9B). During this period, in control plots, q was in downward direction with the lowest q of −0.040 ± 0.026 mm day⁻¹ during the rainy season 2010 (June–September). In thinned plots, q was predominantly in upward direction with the highest q of 0.134 ± 0.054 mm day⁻¹ during the dry season in May 2010 (Table 4).

Our results suggest that thinning increases soil volumetric water content at 35 cm depth, likely stored at the bottom of the soil column on a rock surface. Following the thinning treatment, the difference in soil volumetric water content between control and thinned plots increased from May 2010 through October 2011, driven by changes in forest structure. Thinning reduced the canopy cover and tree density, leaving open spaces covered by scattered slash. These changes aimed to maintain higher soil volumetric water content for longer periods compared to control plots, even during the lowest precipitation recorded (i.e., the rainy season) in 2011 (Figure 6). The weathered bedrock vadose zone, referred to as “rock moisture” (Hahm et al., 2022; Luo et al., 2024), mitigates drought impacts by storing water in wet years and releasing it in dry periods (Callahan et al., 2022). Forest structure is a key factor influencing soil moisture at the upper 50 cm depth, especially during drought periods (Belmonte et al., 2022). Thinning creates forest gaps, enhancing soil moisture through three key processes: (1) reduced canopy interception (Wang et al., 2015), (2) reduced water uptake (Marthews et al., 2008), and (3) lower transpiration rates (Ritter et al., 2005). Our results agree with Lucas-Borja et al. (2018), who reported warmer soil surface temperature and higher soil volumetric water content in straw mulch plots compared to control plots in a Mediterranean forest. Similar to our study, thinning a semi-arid ponderosa forest increased soil moisture in the first 100 cm vs. unthinned forests (Sankey and Tatum, 2022). Another study reported temporal effects on soil moisture response in thinned stands. The onset of soil moisture depletion was delayed in thinned vs. control stands. However, in years with above-average precipitation, thinned and control stands reached their maximum water storage capacity (Hardage et al., 2022).

Our results revealed a more nuanced effect of forest thinning on soil water flux than initially hypothesized. Specifically, the hypothesis, “In treated plots, water flows in downward direction during the snowmelt and rainy season,” was partially confirmed. It was not supported during the snowmelt period (upward water flow), but was supported during the rainfall period in 2010 (downward water flow).

During the snowmelt season (2010), differences in soil and slash cover temperatures influenced the upward water flux. This was confirmed by the hydraulic head gradients, with a higher water potential at 35 cm and a lower water potential at 7 cm depth. The differences in gradients were likely influenced by high evaporative demand at the surface. In thinned plots, solar radiation reached the slash cover directly, increasing the soil temperature and slash cover, which provided the latent energy for evaporation (Heck et al., 2020). Concurrently, wind removed water vapor from the surface, sustaining a high vapor pressure deficit (Davarzani et al., 2014). Both variables influenced the upward flow observed during this period. Soil temperature gradients influence soil water flux in unsaturated soils; however, the magnitude of the temperature effect depends on the volumetric water content (Bach, 1992). Gierke et al. (2016) characterized the snowmelt in the same watershed as slow, uniform infiltration that primarily fills soil micropores, creating a tightly-bound water source for water uptake during the dry spring. The upward flux we reported was a key mechanism for redistributing moisture in the shallow root zone, although the majority of the snowmelt water likely infiltrated and was stored at the weathered bedrock. Contrary to our results and a different setup, Iwata et al. (2010) reported downward water flux as the snow melted in above-zero temperatures in an agricultural field devoid of vegetation.

During the rainy season (2010), we observed a complex pattern of water flux. It was downward from 7 to 20 cm and upward from 35 to 20 cm. This pattern likely occurred due to the water uptake by pine seedlings known as hydraulic redistribution. This process occurs when water moves from deeper soil layers (moist soil) upward to shallow soil layers (drier soil) through the root system, usually during the night when transpiration has diminished (Richards and Caldwell, 1987; Brooks et al., 2002). This is a plausible process on our site, where ponderosa pine seedlings can extract water from shallow, rocky soils (Rose et al., 1997) despite root growth being influenced by soil moisture and bulk density (Siegel-Issem et al., 2005). According to the tree density assessment conducted in March 2010, pine seedlings averaged 163 per ha in thinned plots and 128 per ha in control plots. Fornwalt et al. (2017) reported similar total understory plant cover on mulching plots (i.e., similar to our treatment) vs. control plots after 2–4 years of the treatment in a mixed conifer forest. Warren et al. (2011) reported that root water extraction was mainly in the upper soil horizon (20 cm) in old-growth ponderosa pine. Isotope analysis from the same watershed showed that trees take water from shallow and deep layers during the rainy season (Gierke et al., 2016). Their results detailed a seasonal pattern where trees relied exclusively on the shallow water pool during the dry spring, then shifted to a mixture of both shallow and deep water sources during the monsoon, and continued to access the deeper water for months after the rains ended.

Our hypothesis, “During the low water input (dry season), treated plots may have upward water flux,” was confirmed. The upward water flux likely resulted from a high evaporative demand (similar to the observed snowmelt season), which created a hydraulic gradient that moved water upward from storage at the base of the rock column. This process reflected water redistribution from the weathered bedrock interface toward a drier surface layer. The upward water flux was observed right after the thinning treatment from October 2009 (dry season) through the end of the study (October 2011), with noticeable differences in the hydraulic gradients in the shallow soil matrix. The scattered slash covered the topsoil and kept the soil warmer compared to the soil under the canopy in control plots. The slash absorbed and re-radiated solar radiation throughout the year, creating a warmer microclimate at the soil surface. This increase in temperature likely enhanced water loss via evaporation. Solar radiation has been widely studied, resulting in warmer soil temperatures in treated than in control plots due to the lack of canopy cover (Everett and Sharrow, 1985; Breshears et al., 1998; Tang et al., 2005; Moroni et al., 2009; Garduño et al., 2010). In addition, forest floor residues (e.g., slash) have been shown to influence soil temperature patterns (Liechty et al., 1992; Zabowski et al., 2000). Similar to our results but different thinning setup, Sun et al. (2016) reported increased soil evaporation after strip thinning in a Japanese cypress forest. They reported that radiation was highly correlated, followed by vapor pressure deficit during this process. Tatum et al. (2025) measured the soil water potential in a thinned Ponderosa forest. The vertical gradient in soil water potential observed in the thinned treatment suggests upward water movement from deep soil layers (less negative at 100 cm) to the shallow root zone (more negative at 25 cm) during the dry season.

The interpretation of our findings on soil water dynamics should be considered within the inherent context of hydrological data uncertainty. Hydrological data are subject to multiple sources of uncertainty, with typical magnitudes of 10 to 40% (McMillan et al., 2018). First, the Darcy-Richards equation relies on key assumptions: (1) water flow is essentially one-dimensional in the vertical direction, (2) the soil is isotropic with uniform hydraulic properties horizontally, and (3) steady-state conditions are approximated over the measurement interval (Warrick, 2003). This approach offers a reliable estimate of shallow soil water fluxes, even though natural soil systems may differ from these ideal circumstances, introducing potential error.

Second, our study considered the uncertainty associated with the scale-up process (using a limited number of experimental blocks, n = 4). Although this design is common in large-scale experiments, it introduces uncertainty when extrapolating results across broader and heterogeneous landscapes. The statistically significant treatment effects observed in soil volumetric content and shallow soil water flux were detectable across this underlying uncertainty and represented a relevant hydrological response to thinning.

The hydrological response to thinning observed in our study is specific to the conditions of a mixed conifer forest overlying fractured limestone bedrock. The direct application of our findings is most relevant to forest and water managers operating in similar climate and shallow, rocky soil environments. Our findings provide critical-grounded guidance for managers aiming to enhance forest resilience and water storage during drought. Despite an increase in soil water storage, thinning had nuanced effects on water flux, mainly in upward direction. This process translates into more stored water and less water lost via evaporation under soils with scattered slash. For managers, this implies that strategic thinning can be a proactive strategy to mitigate drought stress in conifer forests. By reducing stand density, managers can effectively delay the onset of soil water deficit, thereby reducing tree mortality and susceptibility to pests, which are often exacerbated by drought.