The physical basis for ecohydrologic separation: the roles of soil hydraulics and climate

The degree of water mixing in the critical zone is under intense debate. Field measurements of isotope ratios indicate varying degrees of separation between pools of water that supply streams and vegetation. The exact physical mechanisms behind ecohydrologic separation are unknown, but local conditions such as soil heterogeneities likely inuence the extent of mixing and separation of subsurface water pools. Using a well-established soil physics model, we simulated if isotopic separations occur within 650 distinct congurations of soil properties, climatologies, and mobile/immobile soil-water domains. Simulations demonstrated separations in isotope ratios between storage and drainage waters during periods of high precipitation, soil water content, and drainage. Separations grew with larger immobile domains and, to a lesser extent, higher mobile-immobile transfer rates. Across soil types and climates, lower saturated hydraulic conductivity and higher rainfall rates amplied isotopic differences, illustrating how mobile and immobile domains interact with local conditions to physically result in subsurface separations. These results show how different critical-zone solute uxes can be generated by representing contrasting transport dynamics in distinct domains across a range of soils and climate conditions. water immobile head gradient and a mass transfer coecient ( ) 26 . We investigated ve model congurations characterizing different heterogeneous mixing and transport: 1) a model representing a single-pore domain with no f and the ω equaling 1 (0 f 1 ω ), 2) a model with a high f representing 40% of the total saturated water content (VWC sat ) and high ω equaling 0.75 (H f H ω ), 3) a model with a low f representing 20% of the VWC sat and high ω equaling 0.75 (L f H ω ), and 4) a model with a high f representing 40% of the VWC sat and low ω equaling 0.25 (H f L ω ), and 5) a model with a low f representing 20% of the VWC sat and low ω equaling 0.25 (L f L ω ). These ve model congurations were tested across a range of different precipitation inputs and soil hydraulic properties, which were varied in relation to the original soil column congured from observed datasets from Watershed 10, H.J. Andrews Experimental Forest. In total, we present 650 of HYDRUS-1D with varying f, ω, VWC sat , K s , total accumulated precipitation, and amount effects describing the correlation between precipitation amount and its isotopic composition. When varying soil hydraulic parameters across model congurations, all other parameters were held constant and VWC sat or K s was altered. We increased or decreased the total accumulated precipitation by multiplying each precipitation event in the input time series by a specic percentage. New precipitation isotope ratios were generated (refer to the Supplemental Material) with stronger or weaker negative correlations between precipitation depth and its isotopic ratio. All model simulations had relatively low mass balance errors (average relative error = 0.09 %) and solute balance errors (average relative error = 0.001 ‰ calculated

zone and the spatiotemporal variability of water and solute transport. Sprenger et al. 24 investigated separations between soil water pools using a one-dimensional ow model at three long-term northern latitude research sites and found modeling with two pore domains improved the representation of soil water isotope dynamics when compared to eld measurements. We expand on this work by utilizing a fully mechanistic modeling approach to examine how separation phenomena arise across a wide range of soil hydraulic parameters or climate conditions to address which environmental conditions do (or do not) lead to the separation phenomena that have been observed 1,16,25 to be intrinsic to soil-water ow.
We con gured 650 geochemically enabled soil physics models (HYDRUS-1D 26,27 ) to evaluate soil water transport and mixing across a range of soil types and with simulated precipitation from different climates. This approach explicitly allowed for heterogeneous mixing of solutes at depth via a dualporosity method. Speci cally, our research objective was to test how pore heterogeneity alone (i.e., without confounding evaporative or instrumental effects 14 ) can give rise to separations in soils across a range of soil types and climates. While the subsurface processes conceptualized in advanced models, such as HYDRUS-1D, represents current mechanistic knowledge within the soil-physics community, the high degree of numerical discretization coupled with the large number of non-linear functions programed into these models, requires an evaluation of these models across an ensemble of con gurations to deduce those emergent priorities inherently arising from a chosen set of soil physics. We hypothesized (1) isotopic separations should arise between water pools if a soil has heterogeneous porosities and (2) the degree of separation between water pools is controlled by soil hydraulic properties and climate conditions. This work extends beyond previous studies by identifying how interactions between climate and porosity heterogeneity within a soil pro le drive incomplete mixing in soils. These results help frame how we understand and represent solute transport and mixing in the subsurface.

UNDERSTANDING DUAL-POROSITY ISOTOPE SEPARATION IN SOILS
Isotopic separations were simulated across different hypothetical conceptualizations of water transport for a single soil type with high, low, and zero immobile water fractions and high and low mass-transfer rates between the mobile and immobile domains (Fig. 1). Across the range of fractions (f) of pore space in the immobile domain and mass transfer coe cients (ω) between the mobile and immobile domains, soils with immobile fractions exhibited distinct transport processes that were not apparent with a single porosity soil con guration (Fig. 1c-e). Given that these models were driven with the same inputs and all other model parameters were held constant, the temporal patterns which arose were attributable to differences in simulated pore heterogeneities. Immobile soil water has limited d 2 H variability regardless of precipitation driven changes in volumetric water content (VWC) or drainage rates (Fig. 1d).
This contrasts with mobile-soil and drainage waters, of which d 2 H varied considerably during wetter and drier simulations periods. When parameterizations included large fractions of immobile soil water, vertical water transport reduced, matrix holding capacities increased, and mobile water fractions decreased; however, this manifested in more rapid breakthrough curves that preceded the low-f breakthrough curves (see days 250-260, Fig. 1e). These differences in tracer breakthrough curves occurred because vertical movement within the soil column and mixing of precipitation inputs was limited to the mobile domain, and thus smaller f values implied more mixing within the column. Varying ω did not in uence the breakthrough curve response. While the separations between pools were distinct, overall differences in δ 2 H from each domain and model con guration were relatively small (all δ 2 H ratios ranged from -70.5 to -65.3 ‰).
Consistently across soil con gurations, as the VWC bulk (mobile + immobile VWC) and drainage increased, the separation between water pools shifted so the drainage water more closely re ected the mobile and immobile water (i.e., the difference between drainage and mobile or mobile water 0 ‰; Fig. 2b,c). During drier periods, stronger separations were observed between soil water pools and the drainage water (Fig.  2b,c). This difference was likely a function of the incomplete mixing of the soil column with new incoming precipitation and the time-lag between the drained and mobile soil water pools. As VWC bulk and drainage increased, models with a higher immobile fraction had the largest increase in isotopic separation (i.e. more positive slopes in Fig. 2b,c; Supplemental Tables S1,S3). All dual-porosity models had positive slopes (most had signi cant p-values < 0.05; Supplemental Table S1) and models with higher immobile fractions had steeper slopes. Interestingly, the separation between the soil and drainage water simulated by a single porosity model had a negative slope as VWC bulk increased , this was distinctly different from all other models (Fig. 2b). Within model variability in these positive and negative slopes arises as a function of the ten stochastically simulated isotope precipitation time series which served as input for each of the modeled porosity heterogeneities. Changing the incoming precipitation would likely change the slopes, however differences among modeled separations should still be apparent across soils with smaller or larger immobile fractions. Supplemental Tables S1, S2, and S3 summarize the linear regression, p-value, Pearson correlation coe cient, and Spearman correlation coe cient calculated between isotopic separations and VWC bulk , daily accumulated precipitation, and drainage. In summary, periods with drier soil water contents will have larger separations between soil water pools and this relationship is ampli ed for soils with larger immobile fractions.

DUAL-POROSITY ISOTOPE SEPARATION ACROSS SOILS AND CLIMATES
Moving beyond a representation of one soil type, we explored a wide range of soil and climate characteristics with varying immobile fractions and transfer rates. Changing a soil's total saturated water content (VWC sat ) and saturated hydraulic conductivity (K s ) in uenced the transport and mixing of soil water, thus producing different degrees of separation between mobile, immobile and drainage pools ( Fig.  3a-d). Separations between the mobile and immobile pools were more pronounced with lower VWC sat and K s (Supplemental Fig. S1). Lowering the VWC sat produced the largest separations in soils with high f and ω (Supplemental Fig. S1a), accordingly large tightly bound water pools and low VWC sat drive these larger separations. Furthermore, soils with large tightly bound water pools and low VWC sat demonstrated decreases in the degree of separation between drainage and mobile water, and this decrease was not observed for models with smaller tightly bound pools (Fig. 3a). Lowering VWC sat and increasing the strength of transfer (ω) between the mobile and immobile domains increased the separation between drainage and immobile water (Fig. 3b). Changing K s in uenced the separation between soil water pools as well, wherein soils with lower K s exhibited larger separations between drainage and tightly bound water (Fig. 3d). Varying the K s of the soil had a minimal effect on the difference between drainage and mobile water (Fig. 3c). Based on the 2-sample t-test, none of the dual-porosity models were statistically different compared to the single pore model (p-values > 0.5; Fig. 3). It is worth mentioning that all simulated differences in Fig. 3a-c varied by approximately 0.5 ‰ and this difference would be di cult to observe from eld observations given typical instrumental accuracy. Nevertheless, we demonstrated that appropriately representing mobile and immobile pools is especially important for soils with low VWC sat and K s (such as with clays) when representing subsurface transport processes and drainage.
Next, we showed that precipitation input amounts can drive the degree of separation between mobile, immobile, and drainage pools (Fig. 4). We both increased and decreased the total input precipitation in simulations, since precipitation drove increases in VWC and drainage in the previous H.J. Andrews soil results. Soils receiving less precipitation yielded lower variances in the separations between drainage, mobile, and immobile pools (Fig. 4a,b, Supplemental Fig. S2a). For these simulations, the soil had a high K s (65.64 cm d -1 ) and therefore smaller precipitation events were less likely to saturate the soil column.
This in uenced the rate of mass transfer between mobile and immobile regions, which was driven by a pressure head gradient, thus decreasing the separations simulated at lower precipitation rates.
Given that precipitation amount controls both separation (as discussed in the previous paragraph) and stable isotope ratios 4,17,28,29 , we posit that the interaction between these two relationships could be a signi cant factor contributing to the ecohydrologic separation phenomena. We tested if precipitation driven separations were disparate between single and heterogeneous porosity landscapes. We show that varying the strength of the amount effect on precipitation isotopic concentrations, which describes the often negative correlation (ρ) between precipitation amount and its isotopic concentration, in uences separations between water pools (Fig. 4). A large amount effect (more negative ρ) means higher intensity rainfall will have much lower δ 2 H ratios than low intensity rainfall. ρ can vary depending on source precipitation (e.g., oceanic versus continental precipitation), latitude, and seasonal temperatures or precipitation patterns 30 . Observed daily values of ρ can range from +0.24 to -0.47 31 and an amount effect of zero removes any correlation between precipitation intensity and its isotopic composition. The most negative precipitation events were simulated when the amount effect was -0.4 and the input precipitation amount was large. Consequently, when ρ = -0.4, we simulated the largest separations between mobile and immobile soil water (Supplemental Fig. S2b) and the drainage and mobile soil water (Fig. 3c). The larger separations can be attributed to the large, negative precipitation events entering the mobile soil domain later in the time series when precipitation intensity was the largest. These large, negative precipitation events also shifted the degree of separation between drainage and immobile pools (Fig. 4d).
Thus, regions with large storm events and large amount effects (i.e., rain-out effects), seasonal precipitation which transitions from snow to rain, or intense temperature changes between seasons are more likely to exhibit larger separations between pools and uxes.

ECOHYDROLOGIC IMPLICATIONS
Many research studies have studied the effects of ecohydrologic separations [1][2][3][4]6,15,[20][21][22] and suggest different storage pools supply plant transpiration (immobile) and groundwater recharge and stream ow (mobile) 1,25 . Here, we demonstrate how soil properties such as a low VWC sat and high f can drive larger separations between pools supplying plants and other water pools. Our results showed the largest separations between drainage and mobile or immobile water were ~4‰ in δ 2 H (Fig. 2a) and when averaged over time, were ~1-2‰ in δ 2 H (Fig. 3a-d, 4a-d). We acknowledge only a few soil physical properties were investigated, and other factors which were not simulated in this study could also drive separations. Literature has cited varying degrees of separation between soils and streams, often far exceeding 4 ‰ 1,32 . However, recent literature has shown that these separations are likely overestimated 7,12 (also refer to Allen and Kirchner 33 ). Moreover, our simulated values cannot be directly compared to previous studies since stream ow is largely composed of temporally integrated groundwater at any given time and there are much longer time lags involved in stream ow generation 16,33 . Regardless, our results show that soils represented with heterogeneous porosities are likely to produce larger separations between soi1 and drainage water, and thus such soils are likely to manifest in previously described isotopic differences between soil water and stream ow.
The implications of varying of the size and in uence of f and ω for tracer transport and mixing have not been explicitly investigated. Prior work by Hu et al. 23 used the simple assumption that half of the total storage was mobile water and no mass transfer occurred between regions, nevertheless they found improved transit time estimates when compared to a complete mixing model. Here we clarify how increases in f increased the proportion of the soil matrix characterized by porous ow media, thus decreasing the amount of incoming precipitation likely to bypass the matrix and decreasing a soil's hydrologic connectivity. The largest separations between drainage and mobile or immobile water were simulated with large f values (Fig. 2b,c). Lowering ω decreases the exchange between the matrix and preferential ow, further limiting mixing and transfer between domains especially under dry soil conditions. A high ω often increased the separation between drainage and immobile pools (Fig. 3b,d). Based on our ndings and previous work, we hypothesize that soil representations accounting for both f and ω will better characterize transport and mixing, as well as prescribing some subsurface water to reside longer (e.g., increase in f will increase the time water travels through the matrix), while other con gurations will result in more water transported through preferential ow paths.
While not explored here, other factors could drive separations quoted in previous studies. These include accounting for evaporation and the corresponding isotopic fractionation mechanisms 34 or root water uptake from the mobile or immobile domains. It also should be noted, plant and soil water extraction techniques can introduce biases and uncertainties in stable water isotope analyses 7,12,33,35,36 , and many reported simulated differences are 1-2 ‰, which is approximately the precision of many laser-based systems when measuring water samples for δ 2 H concentration. Consequently, the simulated differences between domains would not be able to be detected in nature.
We provided the rst comprehensive evaluation of transport in the subsurface under different soil porosity heterogeneities across a range of ne to coarse soils and wet, dry, and seasonally varying climate conditions. Future research should build on this analysis and investigate other mechanisms (e.g., evaporation, root water uptake) that might drive larger separations between soils, plants, and stream ow. Regardless, we have demonstrated that separations can be explained by heterogeneous mixing in soils alone and how these separations are expected to vary with local ecohydrologic characteristics. Thus, any models aiming to realistically represent transport processes, especially those characterized through tracer observations, must represent the heterogeneity of soil pores within a soil column and the exchanges among them. This analysis demonstrated that the coexistence of both ner and coarser pores within a single soil pro le by itself can manifest in real, complex tracer phenomena that have otherwise been attributed to myriad ecophysiological and hydrological processes.

Methods
The one-dimensional HYDRUS-1D dual-porosity numerical model 26 with modi cations for stable water isotope transport 27 was con gured to represent 650 different porosity heterogenies across soils with different hydraulic properties and climate conditions. HYDRUS 26 is one of the most widely used numerical models for simulating the movement of water, heat, and solutes in various soil conditions and has been used in other isotope studies in a single porosity con guration 14,20,22,27,33,37 . We rst con gured a dual-porosity HYDRUS-1D model with parameters based on observed soil hydraulic properties from Watershed 10, H.J. Andrews Experimental Forest, Oregon, USA 1 and simulations of VWC bulk were compared against observed VWC bulk from September 14, 2006 to December 23, 2006. These 100 days were selected for our analysis because precipitation events during this period were sampled for their isotopic concentration in 5 mm increments, making this a unique dataset on which to test our hypothesis since many other sampling approaches integrate precipitation sampling over a week or longer time scales. The observed precipitation samples were used within a statistical downscaling method 31 to simulate realizations of possible precipitation inputs corresponding with observed precipitation amounts. The model represented the top 100 cm of the soil pro le, no evaporative effects were considered from the surface, and no root water uptake was simulated within the column. All pore water heterogeneities were driven by model parameterization of the soil physical properties. Refer to the Supplemental Material for further details on model con guration, parameterization, and initial conditions and the simulated input precipitation datasets.
The modeled "mobile" soil water domain represented regions of the soil matrix such as preferential ow paths or large pore-spaces and the modeled "immobile" soil water domain represented tightly bound water held at water potentials below what can drain by gravity. Water was only transported vertically in the mobile region and water movement into or out of the immobile region was controlled by the pressure head gradient and a mass transfer coe cient ( ) 26 . We investigated ve model con gurations characterizing different heterogeneous mixing and transport: 1) a model representing a single-pore domain with no f and the ω equaling 1 (0 f 1 ω ), 2) a model with a high f representing 40% of the total saturated water content (VWC sat ) and high ω equaling 0.75 (H f H ω ), 3) a model with a low f representing 20% of the VWC sat and high ω equaling 0.75 (L f H ω ), and 4) a model with a high f representing 40% of the VWC sat and low ω equaling 0.25 (H f L ω ), and 5) a model with a low f representing 20% of the VWC sat and low ω equaling 0.25 (L f L ω ). These ve model con gurations were tested across a range of different precipitation inputs and soil hydraulic properties, which were varied in relation to the original soil column con gured from observed datasets from Watershed 10, H.J. Andrews Experimental Forest.
In total, we present 650 simulations of HYDRUS-1D with varying f, ω, VWC sat , K s , total accumulated precipitation, and amount effects describing the correlation between precipitation amount and its isotopic composition. When varying soil hydraulic parameters across model con gurations, all other parameters were held constant and VWC sat or K s was altered. We increased or decreased the total accumulated precipitation by multiplying each precipitation event in the input time series by a speci c percentage. New precipitation isotope ratios were generated (refer to the Supplemental Material) with stronger or weaker negative correlations between precipitation depth and its isotopic ratio. All model simulations had relatively low mass balance errors (average relative error = 0.09 %) and solute balance errors (average relative error = 0.001‰) calculated by HYDRUS-1D during the numerical computations 26 .
This study represented 100 days in a Paci c Northwest winter wet season 1,32 and not a full year with dry periods, large evaporative effects, or high transpiration rates. The 100-day time series was repeated three times to remove the effects of the initial condition of the soil's stable water isotope signature, and the nal 100 days (days 200-300) were analyzed in the presented results. We chose this approach to reduce the confounding effects of evaporation and we were constrained by the computation limitations imposed by running HYDRUS-1D with many con gurations. The observation dataset from Brooks et al. 1 was used to calibrate our modeled soil properties, however our objective was not to exactly simulate the observed water isotope datasets. The simulated stable water isotopes in precipitation are not equal to the actual precipitation during the period and we did not know the initial soil water's isotopic concentration. Our focus was to test the dual-porosity approach across 650 representations of transport dynamics across soils and climates.
Declarations expressed in this paper are those of the author(s) and do not necessarily re ect the views or policies of the U.S. Environmental Protection Agency. Mention of trade names or commercial products does not constitute endorsement or recommendation for use.

AUTHOR CONTRIBUTIONS
CEF, SPG, and JRB designed the study and developed the methodology. JRB provided observation datasets. CEF performed the analysis and wrote the manuscript. SPG, JRB, STA, and SS provided guidance and reviewed the manuscript. Figure 1 Time series of simulated precipitation, soil water pools, and drainage isotope ratios. (a) Averaged (+/-std from the 10 simulations) downscaled precipitation's stable water isotope ux-weighted concentration (δP_input) and the daily input precipitation. Averaged time series of the 10 simulated (b) volumetric water contents (VWC) of a single porosity soil (VWCbulk), the mobile (solid lines) and immobile (dashed lines) soil domains for high (Hf, dark grey, 40% of VWCsat) and low (Lf, light grey, 20% of VWCsat) immobile fractions, and the drainage from the column (green), (c) isotope ratio of the mobile soil water domain, (d) isotope ratio of the immobile soil water domain, and (e) isotope ratio of the column's drainage for the soils simulated with a single pore domain model (0f1ω, black), high immobile fraction, f, models with high(0.75)/low(0.25) transfer rates, ω (HfHω, red; HfLω, orange) and low immobile fractions with high/low transfer rates (LfHω, purple; LfLω, blue).

Figure 2
Instantaneous δ2H differences between soil pools and uxes during different soil moisture conditions. The separation for 10 simulations in (a) isotopic concentration of the mobile and immobile soil water at different volumetric water contents (VWCBulk) for models with high and low immobile fractions (Hf and Lf) and transfer rates (Hω and Lω). The separation in isotopic concentration of the drainage from the column and (b) mobile soil water and (c) immobile water at different water contents. Refer to Figure S1 in the Supplemental Information for (a-c) evaluated against depth of precipitation into and drainage from the soil column.

Figure 3
The in uence of modeled soil hydraulic parameters on average δ2H difference between drainage and (a, c) mobile or (b, d) immobile soil water. Each boxplot represents ux and volumetrically weighted averaged differences, the black triangles indicate the mean calculated from 10 simulations, and the diamonds are outliers. Missing boxplots at 0.35 cm3cm-3 (a-b), 0.40 cm3cm-3 (a-b), and 40 cm day-1 (cd) indicate where Hydrus model con gurations failed to converge. Refer to Fig. S1 in the Supplemental Information for the average isotopic difference between mobile and immobile soil water evaluated against the soil hydraulic parameters in (a-d).