Testing the ‘two water worlds’ hypothesis under variable preferential flow conditions

Widespread observations of ecohydrological separation are interpreted by suggesting that water flowing through highly conductive soil pores resists mixing with matrix storage over periods of days to months (i.e., two ‘water worlds’ exist). These interpretations imply that heterogeneous flow can produce ecohydrological separation in soils, yet little mechanistic evidence exists to explain this phenomenon. We quantified the separation between mobile water moving through preferential flow paths versus less mobile water remaining in the soil matrix after free‐drainage to identify the amount of preferential flow necessary to maintain a two water world's scenario. Soil columns of varying macropore structure were subjected to simulated rainfall of increasing rainfall intensity (26 mm h−1, 60 mm h−1, and 110 mm h−1) whose stable isotope signatures oscillated around known baseline values. Prior to rainfall, soil matrix water δ2H nearly matched the known value used to initially wet the pore space whereas soil δ18O deviated from this value by up to 3.4‰, suggesting that soils may strongly fractionate 18O. All treatments had up to 100% mixing between rain and matrix water under the lowest (26 mm h−1) and medium (60 mm h−1) rainfall intensities. The highest rainfall intensity (110 mm h−1), however, reduced mixing of rain and matrix water for all treatments and produced significantly different preferential flow estimates between columns with intact soil structure compared to columns with reduced soil structure. Further, artificially limiting exchange between preferential flow paths and matrix water reduced bypass flow under the most intense rainfall. We show that (1) precipitation offset metrics such as lc‐excess and d‐excess may yield questionable interpretations when used to identify ecohydrological separation, (2) distinct domain separation may require extreme rainfall intensities and (3) domain exchange is an important component of macropore flow.


| INTRODUCTION
The 'two water worlds' hypothesis (McDonnell, 2014) has recently been proposed as a way to explain distinct isotopic signatures of mobile soil and stream water versus tightly bound soil water used preferentially by plants (Goldsmith et al., 2012;Hervé-Fern andez et al., 2016;McDonnell, 2014). This hypothesis, described alternatively as ecohydrological separation (Brooks et al., 2010;McCutcheon et al., 2017;McDonnell, 2014;Pit, 2011), has now been observed across climates and regions (Evaristo et al., 2015;Goldsmith et al., 2012;Good et al., 2015), and at a fundamental level suggests that water in mobile storage can resist mixing with less mobile water over periods of days to months. This conceptual framework challenges the conventional understanding that new precipitation replaces older storage in soils and groundwater en route to streams (Hewlett & Hibbert, 1967;Horton & Hawkins, 1965). This mechanism of translatory flow displacement is often accompanied with idealized assumptions of mixing belowground (Buttle, 2020). Likewise, many hydrological models assume that soil waters represent a well-mixed reservoir whose residence time increases with depth (Sprenger, Stumpp, et al., 2019) and that recent rainfall supplies transpiration via passive uptake of infiltrating water (Dubbert et al., 2019;Sprenger, Stumpp, et al., 2019) or groundwater (Barbeta & Peñuelas, 2017;Fahle & Dietrich, 2014;White, 1932). However, recent global analyses of water isotope signatures have suggested that soil and planttranspired water pools are disconnected from mobile waters that supply streamflow (Bowen, 2015;Evaristo et al., 2015;Good et al., 2015). These studies ultimately argue that some fraction of soil pores conduct water through the vadose zone at velocities that kinetically restrict mixing with more tightly bound storage (Berry et al., 2017).
Instead, recent studies that examined stable water isotopes in soils have suggested that mobile fractions of water obtained from suction lysimeters may not mix with tightly bound soil matrix water (Sprenger, Llorens, et al., 2019;Zhao et al., 2013) or plant-extracted water  over periods of multiple months. As field observations of domain separation become more numerous, the need for mechanistic explanations grows (Sprenger et al., 2016;Sprenger & Allen, 2020a;Sprenger & Allen, 2020b).
One concern is that common measurement artefacts can bias results. The conceptual constraints posed by the two water worlds hypothesis have thus elicited a re-evaluation of isotope methodology for sampling and quantifying water sources in the vadose zone (Berry et al., 2017;Dubbert et al., 2019;Gaj et al., 2017). For example, typical methods used to sample soil water, such as cryogenic extraction, are often unable to recover water with the same isotopic composition as the labelled reference (Orlowski, Breuer, & McDonnell, 2016;Orlowski, Pratt, & McDonnell, 2016). Further, extracting soil water under relatively low tension via suction lysimeters (e.g., < 60 kPa) versus high tension cryogenic distillation methods (>1100 kPa) may produce distinct isotopic compositions (Landon et al., 1999;Penna et al., 2018;Thoma et al., 2018), which are interpreted as evidence of limited mixing between pores of contrasting size. Isotope labelling studies have even suggested that the soil itself can kinetically fractionate water molecules (Newberry, Nelson et al., 2017;Newberry, Prechsl, et al., 2017) and deplete the soil solution of heavy isotopologues (Aragu as-Aragu as et al., 1995), as inner-sphere mineral complexes preferentially bind 18 O versus 16 O. However, global analyses of ecohydrological separation ignore these small-scale mechanisms, and instead rely exclusively on the premise that water sources can be distinguished by their deviation from global (d-excess) or local (lc-excess) precipitation via evaporative fractionation (Evaristo et al., 2015). If surface soils are additionally altering the stable isotopic signature of precipitation, this challenges the use of dual isotope approaches as ideal tracers in the soil-plant-air continuum and, at the very least, requires more detailed insight using pool-weighted approaches Dubbert et al., 2019;Sprenger et al., 2016).
Another explanation is that non-equilibrium or preferential flow may induce different isotopic signatures between fast-moving, mobile water and slow-moving, immobile water. Water within these preferential flow pathways, being held under relatively low tensions, may become over-represented in suction cup lysimeter and similar measurements (Berry et al., 2017;Goldsmith et al., 2012;Sprenger & Allen, 2020b;Sprenger, Llorens, et al., 2019). Subsurface pore topology (e.g., soil structure, macroporosity) is considered to be a primary source of flow heterogeneity in the vadose zone, with rainfall intensity acting as an important factor affecting the proportion of water moving through more conductive pores (Flury et al., 1994;Stewart, 2019). Additionally, as soils approach saturation, the two flow domains can become more distinct (Cueto-Felgueroso et al., 2020;Scaini et al., 2019) and drainage can select for younger water (Rodriguez et al., 2018;Rodriguez et al., 2020;Sprenger, Stumpp, et al., 2019). Under such conditions, it can be expected that structured soils with high macroporosity (e.g., intact soil columns) will exhibit faster dual domain flow compared to soils with limited structural pore space (e.g., repacked soil material), dominated by slower single domain flow through the soil matrix (Gerke, 2006;Köhne et al., 2009a;Köhne et al., 2009b). Despite our knowledge of nonequilibrium flow processes in soils, the role of preferential flow in maintaining ecohydrological separation has not been considered explicitly through experimentation nor has a responsible mechanism been clearly identified. Assuming that (1) ecohydrological separation stems from domain separation (Sprenger & Allen, 2020b) and (2) domain separation arises from a violation of storage homogeneity assumptions (Pfister & Kirchner, 2017), a 'two water worlds' scenario should be replicable under controlled conditions where all sources are quantified and mixing is measured directly.
The main objective of this study was to test the occurrence of a two water worlds scenario in macroporous soils. Here we quantified the separation between mobile water moving through preferential flow paths versus less mobile water remaining in the soil matrix after free gravity drainage. We applied simulated rainfall with varying intensity and isotopic labels to produce and quantify preferential flow under distinct scenarios. We reasoned that a two water worlds scenario would require (1) incomplete mixing between rainwater and matrix water and (2) a significant difference in the amount of preferential flow between a dual-domain system with intact soil structure and a single dominant system with repacked soil material. The key innovation here is our ability to assess the plausibility of persistent and isotopically-distinct mobile and less mobile pore volumes in soils under varying rainfall conditions and flow heterogeneity. Then, by identifying the contribution of preferential flow to apparent domain separation, we aim to better explain or refute the existence of ecohydrological separation.

| METHODS
All soil material was taken from a pasture in Whitethorne, Virginia.
The pasture had a 9-11% slope and was underlain by Shottower series loam (Typic Paleudult). This soil has moderate to strong soil structural features and is particularly prone to preferential transport of water and solutes (Radolinski et al., 2018). Before the experiment, intact cylindrical cores (5 cm diameter Â 5 cm height; n = 6 per horizon) and unconsolidated soil samples were taken from both the A p (0-20 cm) and B t (20-60 cm) soil horizons (n = 6 samples per horizon).
Cores were used to determine soil bulk density [M L À3 ], porosity [L 3 L À3 ], and saturated hydraulic conductivity (K s ) [L T À1 ]. Corederived K s was measured using the falling head method with a UMS KSAT Benchtop Saturated Hydraulic Conductivity Instrument (METER Group, Pullman, WA). Unconsolidated soil samples were air dried, sieved to 2 mm, and analysed for cation exchange capacity (CEC), pH, total organic carbon (TOC), and texture. CEC was measured colorimetrically via ammonium acetate at pH 7 using a Lachat Quickchem 8500 autoanalyser (Lachat, Loveland), soil pH was measured in a 1:1 slurry (soil: CaCl 2 ), TOC was quantified by dry combustion using a Vario MAX CNS macro elemental analyser (Elementar, Hanau, Germany), and textural analysis was conducted via the pipet method (Day, 1965). Soil physiochemical and hydraulic properties are shown in Table 1.

| Column preparation
We conducted a soil column experiment to create hydrological domain separation using simulated rainfall with an isotopic label. We sought to obtain a range of preferential flow conditions from more homogenous, single-domain flow to distinct dual-domain water transport. To this end we (1) manipulated the physical structure of the medium into three distinct conditions and (2) exposed the columns to three rainfall intensities.
Using the same air-dried Shottower loam material for all columns, we imposed three treatments in 60 cm long by 20 cm inner diameter Schedule 40 PVC. One treatment contained minimal soil structure (reduced structure), the second treatment included an intact Bt horizon (intact structure) and the third treatment used artificial macropores to suppress water exchange between domains (noexchange). The reduced structure treatment (n = 6) was comprised of Shottower loam Ap (0-20 cm) and B t (20-60 cm) soil sieved to 2 mm and packed to native bulk densities. The intact structure columns (n = 7) contained an intact Bt soil layer with moderate subangular blocky structure and homogenized Ap material packed to the surface according to packing methods in Radolinski et al. (2018). The noexchange columns (n = 6) were constructed similar to the reduced structure columns, except that we installed a series of artificial macropores using three 120 cm segments of vinyl tubing (0.3 cm OD).
The segments were heat cured at 90 C for 10 min into a 60 cm long, resting spiral conformation ( Figure S1), which was then wrapped with silica yarn (Fibre Glast Developments Corp., Brookville, OH) at a recommended density of 0.3 g cm À3 (Mori et al., 2014;Mori & Hirai, 2013). The chemically inert silica material was chosen to maintain relatively high hydraulic conductivities, promote enhanced vertical infiltration, and resemble natural biopores linings. The macropore tubes were installed into columns by weaving the three tube segments together at a rate of 1 crossing per 10 cm of column length and spanning the repacked Ap and Bt horizons. While packing columns, the uppermost 5 cm of macropores were severed (i.e., no macropores in the 0-5 cm section of the column) to provide near-surface contact with the soil matrix and allow preferential flow to continue after surface ponding subsided. Figure 1 shows a simple schematic of the three treatment designs.
Following packing, all soil columns were fitted with a flat PVC end cap containing a 1.3 cm thread to 0.64 cm barbed hose fitting.
One cm of mixed fine-and coarse-grained sand was placed between the cap and the bottom soil surface to maintain hydraulic connection between the soil and outlet drain. One g of glass wool was added inside the hose fitting and a 250 μm steel mesh screen was sealed above the fitting to prevent fine soil particles from leaching. PTFE hosing connected the barb fitting to a sealed 800 mL glass container, where all drainage was collected.
A series of 50-g air-dried soil column samples from the A (n = 10) and Bt horizons (n = 20) were further dried in an oven at 105 C for 48 h to determine the amount of water remaining in soils before labeling (Table S1). The air-dried soil columns were saturated with deuterium-spiked tap water, δ o (2.90 ‰ δ 2 H and À 7.53 ‰ δ 18 O), and allowed to drain to field capacity for 48 h, with drainage from each column sampled 1-2 times for stable isotope composition. The columns were continuously covered with PVC endcaps and headspace was filled with a sealed 240 mL ice pack every 8-12 h to prevent evaporative losses, taking care to remove any condensate while transferring icepacks. The average water content of air-dried soil columns was determined to be 0.02 m 3 m À3 and the mean absolute differences between drainage samples and the saturation label were 0.73 ‰ δ 2 H and 0.71 ‰ δ 18 O. These low water content and drainage values suggest negligible influence from any residual water present before the pre-experiment.

| Column rainfall experiments and water sampling
We performed three consecutive 38-mm rainfall simulations using a Corning Sprinkle Infiltrometer (Ogden et al., 1997). For 72-h after each event, we collected 2-3 water samples as drainage volume approached the storage capacity of the glass containers. In order to avoid any potential effects of evaporative enrichment on drainage water signature, we removed any samples <40 mL (Prechsl et al., 2015) from subsequent analyses of preferential flow. In addition, the drainage flux through the lower boundary was recorded from 9 of the columns (n = 3 per treatment) by photographing the digital display of load cells at 15-minute intervals for the full experimental period.
F I G U R E 1 Illustration of the three soil structural treatments prepared in this study F I G U R E 2 Schematic of column study rain inputs and their relative stable isotope signatures which oscillate around known saturation values (δ o ) in order to detect matrix water versus preferential water contributions to drainage To best quantify soil matrix water composition, we took soil cores of the A (5-10 cm) and B horizons (20-25 cm and 45-50 cm) for all treatments. Due to the limited number of replicates, we destructively sampled one column per treatment following initial saturation (time 0) and after the 72-h drainage period after Events 1 and 2. We then sampled the remaining columns at the conclusion of drainage after Event 3. This procedure meant that the sample size decreased by one for each rainfall event, from n = 5-6 in the first event to n = 3-4 for the final event. Each column was destructively sampled by aggregating four cores (2 cm diameter) per depth into one aggregate sample. Soil samples were stored in 470 mL mason jars, sealed, and left to equilibrate at 4 C for up to 1 month before analysis.

| Partitioning of outflow
We modelled flow as being partitioned into two distinct hydrological domains: rapid flow through preferential pathways (e.g., root channels and macropores) and slower flow through the soil matrix. Following the conceptual framework given by Stumpp et al. (2007) the water mass balance is: and the isotope mass balance is: where preferential flow, Q PF (t) [L 3 T À1 ], and matrix flow, Q MF (t) [L 3 T À1 ], sum up to total flow Q t (t) [L 3 T À1 ], and C PF (t), C MF (t), and C t (t) correspond to the isotope composition in each flow component.
Under the assumption of no exchange between the preferential flow pathways and the matrix during the duration of a precipitation event, we consider the irrigation isotope signal to be equivalent to the preferential flow signal in the outlet: Thus, the fractional contribution of preferential flow to the outlet signal is calculated as: The per-event fraction of preferential flow was calculated via Equation (4). We considered the matrix signature for each rain event and treatment, C MF (t), to be the depth-weighted soil isotope value measured after the previous event ( We quantified the uncertainty in f PF by applying a Gaussian error propagation method proposed by Genereux (1998) for two component endmember mixing in differential form as: and in analytical form as: where C and U are the volume-weighted (drainage) or depth-weighted (soil water) average tracer concentration and associated uncertainty and for each mixing component per event and treatment at time t. We considered U PF (t) to be equal to the analytical precision of liquid isotope analysis for each respective tracer, based on one rainfall sample measured per event. U t (t) was equated to the standard deviation of volume-weighted drainage samples for each treatment and the T A B L E 2 Details of the three simulated rainfall events, along with assumptions of pre-event soil matrix and preferential flow signatures used to quantify preferential flow  Table S2.
Some of the f PF calculations resulted in negative values, due to C t values being smaller than C MF values. Nonetheless, these negative estimates had a relatively small mean f PF value of À0.06 and fell within the general range of uncertainty for the experiment (Table S2), so we assumed f PF = 0 for subsequent analyses. LGR analyses. The average differences between four external standards run on both spectrometers were <5% for δ 2 H and <1% for δ 18 O; thus, we considered any analytical bias to be negligible between the two instruments.

| Stable isotope analysis
The isotopic signature of soil matrix water was determined via direct vapour equilibration using the LGR spectrometer. Each mason jar was fitted with a sealed cap onto which two Swagelok bulkhead fittings were mounted with air tight gaskets. HDPE tubing connected one bulkhead fitting to the intake port on the LGR via Swagelok compression fitting, while the other bulk head was connected with similar components to the air-discharge port, ( Figure S1) forming a closed system. Water vapour concentration was monitored until the system we used the standard equations of Horita and Wesolowski (1994) to calculate the signature of liquid water in soil. Standard curves were calculated twice daily using the same procedure performed on liquidwater standards provided by Los Gatos.
All samples were expressed in per mil (‰) delta notation relative to Vienna Standard Mean Ocean Water (VSMOW) via where R sample and R VSMOW are the ratios of heavy to light species in the sample and standard VSMOW water, respectively.

| Statistical analysis and data processing
We sought to identify whether different rainfall intensities could cause preferential flow to exhibit a threshold behaviour ( These differences were far more pronounced for 18 O compared to 2 H. While the median 0.81‰ δ 2 H separation between saturation and soil matrix water was less than one unit of lowest instrument precision (i.e., 1.18‰ δ 2 H), the highest δ 18 O separation of 3.4‰ was an order of magnitude higher than our lowest precision (0.33‰ δ 18 O).
These results suggest that the soil medium itself altered the isotopic signature of saturation water, and that this effect differed between the two isotopes. For example, the median matrix water isotopic value after saturation was 0.81‰ more enriched in 2 H than the labelled input source, while drainage water was 0.61‰ more depleted in 2 H than saturation water. This enrichment of matrix water relative to drainage is consistent with kinetic fractionation of water molecules with greater contact to the mineral surface (Newberry, Nelson, et al., 2017;Newberry, Prechsl, et al., 2017). Preferential inner-sphere binding of heavy water isotopes (Gat, 1996;Newberry, Prechsl, et al., 2017) may have been amplified by saturating initially dry soil (Oerter et al., 2014), enriching the measured soil water in 2 H relative to drainage. In contrast, the matrix and drainage water were both enriched in 18 O relative to initial saturation water. We note here that -before labeling-these air-dried soils may have contained highly fractionated matrix-bound water that isotopically enriched the added water. However, even assuming that the bound water resembled extremely fractionated desert soil water (e.g., 14‰ δ 18 O; [Sprenger et al., 2016;Wu, Zhou, et al., 2014]), the mean air-dried water content of 2% (Table S1) (Gaj et al., 2017;Newberry, Prechsl, et al., 2017), we speculate that the enrichment process that altered the soil water isotope values during saturation may have modified the more mobile drainage through Event 1. We minimized kinetic fractionation from evaporative losses and temperature fluctuations by sealing the column tops using ice packs to prevent overheating, and we discarded samples of low volume, thus it is likely that the soil itself altered the isotopic composition of water in flux. We further note that this strong soil fractionation effect occurred without any artificial alteration to soil surface chemistry such as drying at >100 C (Gaj et al., 2019;Newberry, Nelson, et al., 2017). It is possible that 18 O equilibrium fractionation factors between free liquid and vapour used to calculate the soil matrix signature (Horita & Wesolowski, 1994) differ from the actual fractionation factors between bound pore water and vapour (Lin & Horita, 2016)   F I G U R E 3 Depth weighted isotope δ values from the soil matrix and volume-weighted drainage water following initial saturation of columns. The dashed line represents the delta values (‰) of water that was used to saturate the soil-pore space prior to the irrigation experiments and shaded areas represent lowest analytical precision for conservative representation

| Preferential flow increases with more natural soil structure and intense rainfall
Here we assume that, for each simulated storm, 38 mm of event water infiltrated as a uniform wetting front that filled all available pore space. Prior to each event, the soil columns each had approximately 0.10 cm 3 cm À3 of available pore space, based on the difference between mean soil porosity of 0.41 (Table 1) and mean pre-event water content of 0.31. Thus, each simulated storm would have reached 38 cm in depth. This calculation suggest that a uniform wetting front did not reach the bottom of our 60 cm soil columns, and that sampled drainage represented a combination of pre-event matrix storage and preferential flow. Therefore, our assumption that event water detection in drainage represented preferential flow appears to be valid.
We calculated the fraction preferential flow, f PF , based on drainage samples per column and event, which allowed us to estimate the proportion of preferential flow to total drainage throughout each event ( Figure S3). Intact structure and no-exchange treatments displayed similar drainage curves, often producing twice as much drainage as the reduced structure treatment in the first 24 h following simulated rainfall ( Figure 6). Increasing rainfall intensities had only a minor influence on the overall drainage pattern, regardless of soil structure. These drainage patterns also appear to have little bearing on the amount of preferential flow transiting these soils. For example, the reduced structure drainage curves were nearly identical for all three rainfall events, yet the amount of preferential flow increased by >1 order of magnitude between the 26 and 110 mm h À1 storms.
Drainage patterns were most similar at the highest simulated rainfall intensity (Figure 6 and S4); however, the fraction of preferential flow F I G U R E 4 Raw soil water and drainage samples in dual isotope space by event (rows) and treatment (columns). The signature of each rain event is denoted by the corresponding number on the plot where '0' is the saturation label, and each row considers measurements from a separate event whose signature is identified in red differed noticeably between all 3 treatments. The cumulative total drainage (30-40 mm) was similar between all events and treatments.
Preferential flow often persisted after the first 24 h of high drainage rates, though the highest rainfall intensity appears to have restricted most of the preferential flow in no-exchange columns to the first day of drainage. More cumulative preferential flow generally corresponded to greater f PF values detected in the first 24 h of drainage of high drainage rates ( Figure S4), yet preferential flow often persisted through the inter-event drainage period. This result suggests that the preferential flow detected throughout each event had an inherent travel time distribution that was selective to rapid transport (e.g., more event water appearing in drainage in first 24 h) but nonetheless incorporated a relevant proportion of water that required multiple days to drain (Pluer et al., 2020).
The fraction of preferential flow generated throughout each event differed widely, with intact structure having a mean f PF of 0.004 ± 0.01 at 26 mm h À1 versus 0.37 ± 0.18 at 110 mm h À1 (Figure 7). Preferential flow accounted for just 0-1% of drainage under the lowest rainfall intensity, indicating that flow was primarily sourced from pre-event matrix storage. The 60 mm h À1 storm produced the highest experimental f PF values; however, volumeweighted preferential flow was statistically similar between all treatments for this event (one-way ANOVA, p = 0.19) as mixing calculation uncertainty was also highest (Table S2). The columns with intact structure produced a significantly larger fraction of preferential flow compared to both the no-exchange and reduced structure treatments during the 110 mm h À1 storm (Tukey, p < 0.037).
F I G U R E 5 Depth-weighted soil water isotope values in single (top) and dual (bottom) isotope space. In all panels red, blue, and green colours correspond to low, intact B horizon structure, and no-exchange macropore treatments, whereas shapes are also used to distinguish rainfall events in the bottom panel. Dashed lines illustrate a scenario where soil matrix retained isotope composition measured after saturation for the duration of experiment (0% or no mixing). Solid lines depict a running, 100% mixing calculation of incoming rainfall with pre-event storage based on initial matrix δ values and per-event storage volume derived from the water content of destructive samples. Numbers inside of symbols denote samples taken following that respective event number (bottom). Arrows illustrate the temporal evolution of mixing calculations from saturation (event 0) to event 3, as depicted in the black coloured diagram in the bottom panel. Error bars represent standard error (SE) of the mean (n = 3-4) per treatment. Points without error bars depict samples where only one full replicate existed (four separate cores per column aggregated to one)

| Domain separation may require extreme rainfall
The mixing scenarios in this study provided mechanistic insight necessary to assess the plausibility of ecohydrological separation via domain separation in soils. For example, most of the incoming water mixed with pre-event water already present in the soil matrix during the first event (Event 1, 26 mm h À1 , 2 y storm; Figures 4 and 5). Further, similar preferential flow estimates between all treatments during Event 1 indicate that this mixing happened independent of soil structure (Figure 7). At a simulated rainfall intensity of 60 mm h À1 (Event 2, 10 y storm), preferential flow ranged between 8% and 40% for all treatments, yet the soil water isotope data (though noisy) suggested that mixing was prevalent regardless of their soil structure. The final, highest intensity rainfall event (Event 3, 110 mm h À1 , 50 y storm) then produced significantly higher preferential flow in the dualdomain system (intact structure) than the single domain dominant system (reduced structure) (Figure 7). As a result, the intact structure columns had less evidence of mixing within the pore water ( Figure 6).
The final rainfall event appears to have triggered a two water worlds scenario for two reasons: (1) incomplete mixing between matrix and rainwater was strongly apparent and (2) there was a significant difference in preferential flow between the intact and the reduced soil structure treatments. Here, preferential flow was five times greater in the intact structure system, in which macropores conducted water at a rate high enough to significantly reduce mixing, compared to the noexchange system, in which macropores were physically isolated from the less conductive matrix. Thus, the latter two events depict a transition from uniform single domain flow to one that favoured a disjunction between storage in the soil matrix and more rapid flow in macropores.
Though we were able to generate domain separation in this study, the conditions required to produce this disjunction elicit closer evaluation. Under a 2-y storm, pore water velocities appear to have been slow enough to facilitate exchange from coarse to fine pores per conventional mechanistic model assumptions (Horton & Hawkins, 1965;Van Genuchten, 1980). Because domain separation required a 50-y storm, there is a 98% (or greater) chance that a given rainstorm will not trigger a two water worlds scenario in this field soil. Thus, it F I G U R E 6 Cumulative drainage following each rain event. The propagation of preferential flow through time per event is represented as a proportion of total lower boundary flux. We (1) calculated the fraction of preferential as a function of cumulative drainage volume using isotope delta values and volumes from drainage samples collected from each column ( Figure S3), then (2) computed a time-integrated fraction of preferential flow using drainage water fluxes recorded the recorded lower boundary flux data ( Figure S4) to as a running average per treatment from 3 to 6 columns over intervals of 6 mm appears more likely that matrix and macropores can thoroughly mix when exposed to extreme rainfall that occurs just once in 10 years (e.g., Event 2, Figure 4), than it is for two water worlds to exist under a 2-y storm (Event 1, Figure 4). Here we note that higher uncertainty in the mixing calculation for Event 2 may partially obscure this assertion. However, even if domain separation was clearly apparent during this event, these rare hydrological conditions appear insufficient to explain widespread interpretations of ecohydrological separation. Worthington, 2019), but some amount of exchange between domains is necessary to initiate and sustain macropore flow (Klaus et al., 2013;Weiler & Naef, 2003). The concept of matrix-derived preferential flow further highlights the importance of domain exchange and disputes the concept of two distinct 'water worlds' in soils.

| Implications and transferability
The two water worlds hypothesis maintains that a mobile domain of soil water transits the vadose zone while avoiding exchange with less mobile storage (Brooks et al., 2010;McDonnell, 2014). This disconnection, driven by contrasts in pore water velocity (Berkowitz & Zehe, 2020) and tension (Berry et al., 2017), is argued to produce two isotopically distinct and ecohydrologically separated pools of water: storage that sustains transpiration versus water that rapidly recharges groundwater and streams (Sprenger & Allen, 2020b). Our results suggest that this physical separation, as described, should be rarely observed in nature. Rather, our study shows that the isotopic composition of soil matrix water is often far more sensitive to changes in rainfall signature than drainage. It also reveals that mixing can prevail in macroporous soils, and that exchange between domains may trigger higher amounts of preferential flow under heavy rainfall (e.g., intact vs. no-exchange, Figure 7). These data further justify the inclusion of storage exchange elements in many dual-domain transport models (Gerke & Van Genuchten, 1993;Germann & Beven, 1985;Jarvis, 1994;Weiler, 2005) despite difficulties observing this phenomenon directly (Hoogmoed & Bouma, 1980;Klaus et al., 2013).
Here we note that this controlled study does not completely capture the full range of field conditions, yet the findings produced from this conservative experimental design should warrant a re-evaluation of ecohydrological separation. For example, our experimental conditions paralleled the wet winter season in Brooks et al. (2010) that produced the highest degree of apparent domain separation in soils. While we did not quantify preferential flow under dry conditions (e.g., lower than field capacity), dry soils would likely amplify any domain exchange (Coppola et al., 2015), as mobile event water is pulled into smaller pores with lower matric potentials (Shao et al., 2017). We also did not include plants in our study, which, depending on evaporative constraints, can increase preferential flow (Luo et al., 2019;Radolinski et al., 2018) or decrease gravitational drainage (Radolinski et al., 2019). Nonetheless, our results are in agreement with studies reporting high degrees of mixing during controlled (Vargas et al., 2017) and field experiments (Dubbert et al., 2019) that incorporated growing plants.
The extreme fractionation of 18 O by soil water at saturation also raises concerns for the use of stable isotopes of water as a conservative tracer for subsurface source partitioning. We found up to 3‰ difference between saturation label and post-saturation soil water signatures, suggesting that soil organo-mineral surfaces may strongly select for heavier isotopologues, independent of evaporative enrichment (Aragu as- Aragu as et al., 1995;Oerter et al., 2014). The strong preferential fractionation in 18 O compared to 2 H complicates the use of precipitation offset metrics such as lc or d-excess to distinguish soil water from ground and stream water. For example, the soil water and drainage water lc-excess values calculated at saturation for our reduced structure columns indicate complete separation of the two pools (by factor of $2), whereas using only 2 H suggests uniform mixing (Figure 8). This phenomenon may artificially inflate apparent separation between plant and soil water versus stream and groundwater observed via analyses of large isotope datasets (e.g., (Evaristo et al., 2015)). Though we used one soil type in our study, soil-derived fractionation of 18 O was observed for all soil depths and treatments across a range of physiochemical (Table 1)  F I G U R E 7 Fraction of preferential flow per treatment separated by simulated rainfall intensity. The gradient depicts a transition from a more homogenous ('one water world') single-domain dominant system versus a more heterogenous 'two water worlds' scenario where soil structure is associated with significantly different preferential flow volumes and incomplete mixing between domains is apparent (e.g., Figures 4 and 5). Error bars represent SE (n = 3-6) and different letters denote statistically significant differences (α = 0.05) hydrological domain separation, at least until mechanisms of flow and surface chemistry can be better resolved.
If a two water worlds scenario does indeed exist in natural soils then the conceptual framework should first be confirmed or rejected at a scale relevant to soil physics and augmented to explain larger hydrological systems. Horton and Hawkins (1965) used pore scale approaches to explain the field observations of celerity versus transport, which were then adopted by Hewlett and Hibbert (1967) to propose translatory flow and variable source area theory, still cornerstones of modern hydrological models (Beven & Germann, 2013). We did not find a definitive mechanism for persistent hydrological domain separation in macroporous soils in this study, and we also were unable to identify a satisfactory mechanism in current literature. Consequently, we call for more controlled hypothesis testing using alternative tracers and soils with varying reactivity to confirm or deny the prevalence of mixing in the vadose zone.