Water isotopic composition traces source and dynamics of water supply in a semi‐arid agricultural landscape

Changes in seasonality and form of precipitation alter the structure and function of grassland and steppe ecosystems and pose challenges for land management and crop production in regions like the Northern Great Plains, North America. This research uses isotopic composition of water (δ18O and δ2H) to explore the sources and fate of soil water in lower‐elevation agricultural areas of the Judith River watershed, in the headwaters of the Missouri River, USA. Extensive non‐irrigated cereal crop production in this area occurs on well‐drained soils and depends on careful water management. Our observations indicate that colder precipitation contributes isotopically distinct water to cultivated terrace soils relative to downgradient groundwaters and streams. Riparian waters also exhibit a higher fraction of contributions from colder precipitation relative to terrace groundwaters and streams. Apparent contributions from colder precipitation in terrace and riparian soil waters suggest that snowmelt is a key component of the water supply to these systems. Riparian waters also show evidence of evaporation suggesting that water spends sufficient time in some ponds and open channels in the riparian corridor to reflect fractionation by evaporation. The evolution of water isotopic composition from soils to shallow aquifers to stream corridors indicates source water partitioning as precipitation moves through this semi‐arid agricultural landscape. The apparent mixing processes evident in this evolution reveal source water dynamics that are necessary to understand plant transpiration, solute processing, and contaminant leaching processes.


| INTRODUCTION
In the semi-arid landscapes of the Northern Great Plains (NGP), trajectories of the seasonality and form of precipitation due to climate change present existential challenges for sustainable management of non-irrigated crop production and grassland ecosystems (Hoell et al., 2020;Lee et al., 2008;Qin et al., 2020).Climate analyses in the NGP and similar ecoregions suggest that the area is experiencing warmer winters, wetter springs, drier summers, and shifts in the composition of annual precipitation from snow to rain (Cross et al., 2017;Xu et al., 2013).These changes in the NGP climate appear to propagate shifts in grassland plant communities and regional greening (Berdugo et al., 2020;Brookshire et al., 2020;Currey et al., 2022).The NGP supports one of the largest areas of intact temperate grasslands (ca.55% of NGP land cover), which are highly vulnerable to economic pressures for conversion to cropland (currently ca.28% of NGP land cover; Auch et al., 2011;Currey et al., 2022;Drummond, 2007;Lark et al., 2015).Land use changes are shifting patterns in the partitioning of precipitation to evapotranspiration, hydrologic storage, or runoff across the NGP landscape (Liu et al., 2019;Miller et al., 1981;Sigler et al., 2020).Understanding the consequences of the trajectories of land use and climatic changes in the NGP will require a detailed understanding of how soil water storage, groundwater recharge, and stream flow generation respond to spatial and temporal dynamics of the precipitation regime.
Snow accumulation and melt drive seasonal and spatial patterns in hydrologic and thermal connectivity that determine the behavior of water and solutes across northern landscapes (Costa et al., 2020;Jasechko et al., 2014;Jencso & McGlynn, 2011;Lundberg et al., 2016;Martin et al., 2018;Pellerin et al., 2012).However, the partitioning of melt from ephemeral snow cover to evapotranspiration, hydrologic storage, or runoff has not been studied as extensively at lower elevations as the flushing effect of melt from higherelevation snow accumulations (Baulch et al., 2019;Boyer et al., 1997;Jencso et al., 2010;Pacific et al., 2010;Pomeroy & Li, 2000).A lack of understanding of snowmelt partitioning limits our ability to predict the trajectories of solute processing and transport at lower elevations.
Lower elevations receive higher proportions of precipitation as rainfall compared to higher elevations that receive higher proportions of precipitation as snow.Snowmelt or rainfall that occurs during colder temperatures when evapotranspiration is limited (Wever et al., 2002) contributes to soil water storage and is particularly important in nonirrigated and annually cultivated agricultural systems.The timing of this increase in water availability stimulates crop productivity and drives solute production and leaching (Kort, 1988;Sigler et al., 2022;Simão et al., 2021;Van Meter et al., 2020).Identifying the contribution of cold versus warm precipitation to agricultural soils at low elevations in the NGP can reveal precipitation partitioning relevant to seasonal drivers of water availability and biogeochemical processes.
Isotopic analysis of environmental water provides a useful tool for assessing the history of water found in soils, terrace aquifers, riparian corridors, and stream channels, which are conceptualized as separate process domains (Montgomery, 1999) that vary in hydrologic and biogeochemical function.Precipitation reflects seasonal variation in the isotopic composition of meteoric water with isotopic values typically lower in colder precipitation than in warmer precipitation relative to the local climate (Clark & Fritz, 1997;Criss, 1999;Dansgaard, 1964;Rozanski et al., 1993;Sharp, 2017).The relationship of δ 2 H with δ 18 O values in precipitation is typically linear with a regionally representative slope, known as the local meteoric water line (Bowen & Wilkinson, 2002;Clark & Fritz, 1997;Dansgaard, 1964;Rozanski et al., 1993;Sharp, 2017).When compared to precipitation, patterns in the isotopic compositions of soil water, groundwater, and stream corridor waters illuminate patterns of infiltration and mixing over seasonal or longer timescales (Liu et al., 2021(Liu et al., , 2022)).Additionally, evaporation fractionates water by preferentially transferring the lighter isotopes of O and H to the atmosphere (Craig, 1961).Because the effect of variation in O isotopes on the molecular mass of water is generally stronger than the effect of variation in H isotopes, the isotopic composition of residual water typically falls below the local meteoric water lines (Gammons et al., 2006;Sprenger et al., 2016).
Considering these potential effects, comparisons of water isotopic compositions across process domains improve our understanding of critical transport and storage dynamics among soils, groundwaters, and stream corridors.
The physical properties of soil shape the response of soil water and solute dynamics to changing precipitation regimes (Rabot et al., 2018;Vogel et al., 2006).Soil weathering creates zones of finer texture and structural units within horizons that alter hydrologic and solute residence times as a function of soil development (Ewing et al., 2006;Maher & Chamberlain, 2014;Schaetzl & Thompson, 2015;Vidon & Hill, 2004).
Soil development is shaped by landform geomorphic and stratigraphic controls (e.g., bedrock, topography) that in turn dictates overall soil water holding capacity (Schaetzl & Thompson, 2015;Sigler et al., 2020;Vidon & Hill, 2004).Soil texture and vertical structure influence water movement and storage, therefore affecting the isotopic composition of 'mobile water' and 'bulk water' in soils (Adams et al., 2020;Liu & Guo, 2022;Sprenger & Allen, 2020).In addition, physical transformation of surface horizons by cultivation and replacement of native perennial vegetation with annual crops transforms soil hydrological processes, which influence groundwater recharge and evapotranspiration in irrigated and non-irrigated systems (Scanlon et al., 2005;Sullivan et al., 2022).Cultivation of soils with inherently limited water holding capacity can lead to increased water loss and nutrient mobilization that impairs ground and surface water quality (Ranalli & Macalady, 2010;Sigler et al., 2020).The inherent development of soils in conjunction with transformations due to cultivation establish soil physical properties that dictate soil hydrology.
Patterns of water isotopic composition across landscape process domains reflect soil hydrological processes and allow for an improved understanding of the feedback between soil development and soil hydrology.
The goal of this research is to understand the contributions of seasonally variable forms of precipitation to soils, groundwaters, and streams across a strath terrace landform managed for agricultural production.We assess water content in cultivated soils and patterns in δ 18 O and δ 2 H values in water samples from precipitation, soils, upland groundwaters, riparian groundwaters, and streams to infer potential causes of variable partitioning of colder precipitation, warmer precipitation, and evaporated waters among these process domains (Figure 1).We define colder precipitation as occurring during lower temperature weather conditions that dictate lower water isotopic values and can occur October through May in northern latitudes and can arrive as either snow or rain.We define warmer precipitation as occurring during warmer temperature weather conditions that can occur from April through October and arrives as rain.We hypothesize that water supply, storage, and hydrologic connections among process domains vary across a landscape to dictate spatial and temporal patterns of hydrologic inputs, mixing, and loss to evaporation.We address this hypothesis through three research objectives: Objective 1: Determine how soil conditions regulate seasonal inputs of precipitation.
Objective 2: Determine how the terrace aquifer mediates contributions of water from soils and delivery of water to riparian corridors.
Objective 3: Determine how the riparian corridor transforms water received from the terrace aquifer.
We investigate these objectives through the assessment of precipitation sources and delivery across landscape process domains to understand soil hydrological processes that drive nutrient transport and biogeochemical processing.et al., 2017).Producers incorporate the practice of fallow in nonirrigated cropping systems by using herbicide to suppress vegetation in a field through a growing season with the goal to increase soil water storage for the following cash crop rotation.

| Study area
The climate of the JRW is semi-arid, with lower elevations receiving mean annual precipitation of 39 cm and a mean annual temperature of 6.6 C (based on the 30-year normal, 1991-2021, PRISM Climate Group).Precipitation occurring in late spring and early summer accounts for nearly half of the annual total (Figures 3,S1,and S2).
Rain is the dominant form of precipitation in this region with snow comprising approximately 7%-25% of the annual precipitation over the period of study (Figures 3 and S2).Severe winds and limited vegetation cover create heterogeneous patterns of snow cover, contributing to uncertainty in determining the total snow-derived water and resulting in substantial spatial and episodic variation in infiltration of snowmelt to soil water storage (Pollock et al., 2018).
Crop production in the JRW occurs on alluvial fans and terraces that host well-drained calcareous soils.Soils are classified as Calciustolls and Argiustolls in the US Taxonomy (U.S.Department of Agriculture, 2022).Soils exhibit generally clay loam textures and relatively shallow gravel zones (40-100 cm) overlying subsurface horizons with >40% coarse gravel content of variable thicknesses extending to between 1.5 and 10 m depths (Sigler et al., 2018).This soil and subsurface architecture results in substantial nitrate leaching from cultivated soils, with observations of pervasive high nitrate concentrations The terrace-riparian-stream continuum on the Moccasin terrace.Conceptual diagram of anticipated precipitation infiltration, mixing, and evapotranspiration processes and subsequent solute processing and transport along a hydrologic gradient such as that of the Moccasin terrace.In this landform, nitrate concentrations are generally 20 mg N L À1 in the terrace groundwater, near 0 mg N L À1 in the riparian groundwater, and 10 mg N L À1 in the stream.
(ca. 20 mg L À1 nitrate-N) in shallow terrace aquifers (Schmidt & Mulder, 2010;Sigler et al., 2018).Initial sampling of precipitation and landscape waters occurred more broadly across the JRW starting in 2012, with subsequent sampling focusing on the Moccasin terrace landform through 2022 (Figure 2 and Tables 1 and S1).The Moccasin terrace is an extensively cultivated strath terrace with a shallow aquifer that is hydrologically isolated from streams or other water sources flowing from the Little Belt mountains (Figure 2; Sigler et al., 2018).The shallow aquifer is perched on relatively impermeable Cretaceous-age marine shales (Vuke et al., 2002).The shale units limit connectivity with deeper groundwater, so water stored in the shallow aquifer is sourced exclusively from local precipitation infiltrating through predominantly cultivated soils and ultimately exported through stream corridors.As a result of this landform structure, water in this shallow aquifer is relatively young, exhibiting apparent water ages of 1-2 years from tritium-helium ( 3 H-3 He) relationships in wells and springs, and average residence times on the order of 20 years estimated from total aquifer volume divided by volumetric discharge (see Supporting Information; Sigler et al., 2018Sigler et al., , 2020)).Terrace groundwaters flowing into stream corridors may contribute directly to the stream channel via channelized surface flow or indirectly via subsurface flow paths through riparian aquifers.Thus, the Moccasin terrace provides a simplified context for studying the movement of water and solutes through terrace soil, into groundwater and ultimately stream corridors (Figure 1).

| Sample collection and handling
Water was sampled at varying frequencies and analyzed for δ 18 O and δ 2 H from 2012 through 2022.Rain and snow cover were sampled near Stanford, Moccasin, and Moore, MT (Tables 1 and S1 and   Figure 2).Precipitation was also sampled at CARC in collaboration F I G U R E 2 Study area and site locations.The Moccasin terrace (center) is located within the Judith River Watershed (inset, top left) and the Northern Great Plains (inset, top center).Detailed sample site locations within each stream corriodor are shown in the far right panels.Sites include terrace groundwater from wells and springs (green), riparian groundwater (purple), and stream water (red).Specific terrace soil water sample locations are omitted for producer privacy.The terrace is geographically isolated from mountain water recharge.Areas surrounding the towns of Moore and Stanford host similar soils and shallow aquifers to the Moccasin terrace but are more connected to mountain water recharge.See specific site locations in Table S1.GIS sources: ESRI, USDA, USGS, FSA, NGP boundary obtained from World Wildlife Fund.
with the Montana Precipitation Isotope Network beginning in 2021 (Tables 1 and S1; Carstarphen, 2021).Water was sampled from terrace soils using porous cup tension lysimeters (PTFE/silica; Prenart Equipment; Frederiksberg, Denmark) at seven locations across a cultivated field near Moore, MT (Tables 1 and S1 and Figure 2).At each sampling location, a single lysimeter was installed in the fine-textured soil horizon just above the contact with the underlying gravel horizon.
Depths of each sampling location range from 60 to 119 cm depending on the thickness of fine-textured soils (see Supporting Information and Sigler et al., 2018).All groundwater and surface water were sampled using a peristaltic pump (Geotech Environmental Equipment; Denver, CO) and pumped through a 0.45-μm capsule filter (Geotech  Environmental Equipment; Denver, CO) into 20 mL glass scintillation vials.Terrace groundwaters were collected from wells and hillslope spring orifices except at two springs, which required sampling from shallow PVC wells 1-2 m up-gradient from the orifice due to stagnant water at the orifice (Tables 1 and S1 and Figure 2; see detailed sampling methods in Sigler et al., 2018).Riparian waters were collected from screened groundwater wells, surface water ponds, and spring channels throughout the riparian corridor sites distributed across the Moccasin terrace (Tables 1 and S1 and Figure 2).Casings of riparian sampling wells were 2.54-cm internal diameter PVC directly inserted with a driving point to ca. 1 m depth in riparian soils.Slots were cut in the PVC to provide 30-50 cm of screening.Terrace and riparian groundwater samples were collected after purging three well volumes.
Stream channel and riparian spring channel waters were collected from turbulent locations presumed to be a well-mixed representation of stream flow.
Samples were either filtered on site using a 0.45-μm capsule filter (Geotech Environmental Equipment; Denver, CO) or filtered in the lab using 0.2 or 0.45-μm syringe disc filters.Samples were collected in duplicate as water volume allowed.Samples were filtered into 20 mL glass or HDPE scintillation vials with conical caps to eliminate headspace, kept on ice during transport, and stored at room temperature in the lab prior to analysis (detailed methods in Supporting Information).

| Isotope analysis
Water samples were analyzed for oxygen and hydrogen isotopic composition using integrated cavity output spectroscopy (Los Gatos Research; Mountain View, CA) at the MSU Environmental Analytical Laboratory (EAL).The relative isotopic compositions of the hydrogen and oxygen atoms in water were first measured as atomic mass ratios of 2 H: 1 H for deuterium and 18 O: 16 O for oxygen-18 (Sharp, 2017).Isotopic ratios were then normalized to their corresponding δ value as the fractional deviation (as per mil or ‰) from the global water standard (Vienna Standard Mean Ocean Water, VSMOW; Craig, 1961;Dansgaard, 1964)

| Local meteoric water line assessment and uncertainty
We calculated a local meteoric water line (LMWL) from the isotopic composition of precipitation at the study site (Equation 1; Dansgaard, 1964): where values for m and b are inferred by linear regression of local pre-  Rozanski et al.,1993) provides a specific regional perspective on predictable differences in the hydrogen and oxygen isotopic compositions in meteoric waters (Sharp, 2017).
To allow assessment of meaningful deviation of environmental waters from the LMWL, the analytical uncertainty estimates of 0.6‰  S3).The Water isotopes database uses the precipitation δ 18 O and δ 2 H dataset from the Global Network of Isotopes in Precipitation (GNIP; IAEA/WMO, 2015) along with latitude and altitude to spatially interpolate precipitation δ 18 O and δ 2 H values at a resolution of 5 0 Â 5 0 across the continents (Bowen et al., 2005;Bowen & Revenaugh, 2003;Bowen & Wilkinson, 2002).Latitude and altitude are geographic controls on temperature (Bowen & Wilkinson, 2002) and vapor pressure which drives relative humidity and determines precipitation isotopic composition (Criss, 1999;Dansgaard, 1964).
We used the model fit to our measurements to develop a pair of volume-weighted annual mean precipitation δ 18 O and δ 2 H values for our study period.The 95% confidence interval of the sinusoidal model fit to our precipitation isotope data encompasses the model fit Water isotope database values (Figure S3).This comparison suggested that the model based on our measurements was reasonable to fill data gaps; therefore, we used it to interpolate isotopic values of precipitation for each day of the year.We multiplied the daily isotopic values by the volume of precipitation for that day to calculate the corresponding precipitation-weighted isotopic value (U.S. Bureau of Reclamation, 2023).The sum of the precipitation-weighted values divided by the total precipitation volume over a given period yields the precipitation-weighted mean isotopic value for that period, which we calculated for each year of the study period.We then averaged these values to determine a single pair of volume-weighted mean δ 18 O and δ 2 H values for precipitation over the entire study period.
Water δ 18 O and δ 2 H values are influenced by evaporation plot to the right of the LMWL because evaporation increases δ 18 O values more than δ 2 H values in the residual water (Craig, 1961).We determined the offset of these reservoirs from the meteoric conditions at the study site by calculating the line-conditioned excess (lcexcess), defined as the difference between the measured δ 2 H value and that predicted by the LMWL regression model for a corresponding measurement of δ 18 O (Equation 2; Landwehr & Coplen, 2006).
More negative lc-excess values suggest the greater influence of evaporative fractionation in a given water volume, relative to the isotopic composition expected in local meteoric waters (Gammons et al., 2006;Sprenger et al., 2016).We calculated lc-excess values for the 2.5% and 97.5% quantiles of the Monte Carlo ensemble for the LMWL (Equation 2) to assess meaningful deviations from zero.

| Statistical analysis across process domains
We used a generalized linear mixed model to assess whether patterns in isotopic data reflect the effect of process domain (i.e., soil water, terrace groundwater, riparian water and stream water) on measured δ 18 O and  Schwarz, 1978;Wagenmakers & Farrell, 2004).Uncertainty in the model may be underestimated on the extreme ends of our data because the normality of the residual assumption was not met for process domains with multimodal results (e.g., terrace groundwater and evaporated riparian water).We performed a 95% family-wise controlled Tukey HSD post hoc test to determine which process domains are associated with a significantly different measured δ 18 O and δ 2 H and calculated lc-excess values when compared to the other process domains.Water content in terrace soils began to decline in early June through the growing season (Figure 3).Overall, the growing season soil water content values exhibited few meaningful responses of volumetric water content to rainfall, presumably due to transpiration dominating the influence of precipitation (e.g., July 2021 across all soil depths; Figure 3).Water began to accumulate deeper in the soil from September through October.Water contents of deeper soils either stabilized or continued to increase through the winter months depending on soil temperatures and contributions from snowmelt (Figure 3).An increase in soil volumetric water content was generally observed when soil temperatures exceeded 0 C in shallow depths and snow was present or precipitation was falling (Figures 3   and S1).Soil temperatures below 0 C interfered with the detection of liquid water additions, thus confounding soil volumetric water content measurements at subfreezing temperatures (Ochsner & Baker, 2008).Therefore, increases in water content with warming above freezing were evaluated relative to conditions prior to the onset of freezing.Snowpack insulates the soil and thus helps maintain temperatures above 0 C, facilitating water movement.Periods of the greatest observed snowpack (>20 mm SWE) coincided with an increase in surface soil temperature to above 0 C and a subsequent increase in soil volumetric water content (e.g., early January 2022; Figure 3).

| Isotopic composition of atmospheric waters
Variation in precipitation isotopic composition in the study area produced an LMWL with a slope of 7.14 and intercept of À10.4‰ (Equation 1; Figure S4).Isotopic values of precipitation oscillated seasonally, where values for snow and rain under colder conditions during the autumnal through vernal periods were generally lower than values for rain under warmer conditions in the vernal through summer periods (Table 2 and Figure 4).The Monte Carlo sensitivity analysis of the LMWL is expressed as a 95% confidence interval around the LMWL (Figure S4) and encompasses the global meteoric water line.
Precipitation values and the modelled volume-weighted annual mean) provide a reference against which we evaluated water sample results from soils, aquifers, and streams (Figure 4).

| Isotopic composition of landscape waters
Combining the data across all years of the study and all sample locations, the variability in terrace soil water δ 18 O and δ 2 H values generally decreased from winter through summer (Figure 4a).Lower isotope values for terrace soil waters were observed in January, March and April (there were no samples collected in February; Figure 4a).These terrace soil water values were lower than those of all other landscape process domains, except for eight terrace groundwater samples (Table 2; see below for terrace groundwater details).
Based on the mixed-model multi-comparison analysis, the population of isotopic values for δ 18 O and δ 2 H in terrace soils was statistically different from that of each of the other landscape process domains (p-value = 0.001 for δ 18 O and p-value <0.0014 for δ 2 H; denoted by the lower-case letters in Figure 5a).  2 and Figure 5a).Riparian surface water samples resulted in some of the highest values observed within this domain and across other process domains (Figure S4).The mean isotope values in stream waters were higher than all other process domains, yet confidence was low in statistical differences with terrace groundwater or riparian water (Table 2 and Figure 5a).Stream water isotope values exhibited limited variation for δ 18 O and δ 2 H values (Table 2 and Figure 5a and 6d).
We interpreted a meaningful influence of evaporation where paired δ 18 O and δ 2 H values of a sample fall below the 95% confidence interval of the LMWL (Figures 5a and 6) and where a negative lc-excess value falls below the 95% confidence interval propagated from uncertainty in the LMWL (Equation 2; Figure 5b).While no process domains had a mean lc-excess value that fell below the 95% confidence interval, terrace soil water lc-excess values exhibited the least negative mean lc-excess value, with a population of lc-excess values that were statistically different from all other landscape process domains ( p-value = 0.023; Table 3 and Figure 5b).Riparian waters T A B L E 2 Summary statistics for δ 18 O and δ 2 H results in precipitation source waters and each process domain.resulted in more negative values outside the 95% LMWL confidence intervals than other process domains, but the mean lc-excess value of this process domain showed little statistical difference from the terrace groundwater and stream water (Figures 5b and 6 and Table 3).

| DISCUSSION
The isotopic composition of terrace soil and riparian waters reflects more distinct contributions from colder precipitation compared to the terrace aquifers and stream channels.Therefore, infiltration of colder precipitation to soil storage may be important to plant productivity and soil biogeochemical processes such as nitrate production and consumption.While the isotopic composition of terrace soil water does not appear to be immediately responsive to individual precipitation events, it is highly variable over space and time and exhibits the most influence from colder precipitation.Terrace aquifers integrate heterogeneities evident in the overlying soils, while stream channels appear to integrate terrace groundwater and riparian water sources.Areas of the riparian corridor exhibit the influence of evaporation, suggesting retention of water for sufficient time and exposure to the atmosphere for this local processing to take place.This analysis illuminates the relative importance of warmer and colder precipitation source waters and their distribution as hydrologic drivers of soil water availability and landscape processes in low-elevation, semi-arid agricultural systems.

| Soil conditions control seasonal variation in infiltration of precipitation
Annual cycles in the volumetric water content, temperature, and isotopic composition of terrace soil water suggest influence from multiple potential water delivery and storage pathways.Soil waters sampled at the bottom of the fine-textured soil horizons just above the gravel layer (60-119 cm depths) generally exhibited a mixture of colder and warmer precipitation waters with limited evaporative influence (Figures 4 and 5 and Tables 2 and 3).Soil conditions conducive to water infiltration were likely variable throughout the winter when soil temperatures in the top 20 cm oscillated between ca.À6 C and +3 C (2019-2022; Figures 3 and S1; see Supporting Information).
Snowmelt infiltrated into soils during periods with the largest observed snowpack (>20 mm SWE e.g., early January 2022; Figure 3).While infiltration into frozen soils is possible and water content measurements of frozen soils are dubious (Gray et al., 2001;Lundberg et al., 2016), increases in apparent soil water content took place after soil temperatures increased above freezing at our monitoring sites.At soil depths below 20 cm, soil temperatures generally remained just above freezing throughout the year (Figure 3), which likely allows for continuous soil water mixing in deeper soil horizons.Rain that falls on snow can deliver a mixture of rain-and snowmelt-derived water into terrace soils.Relative to shallower samples, the more limited variation in the isotopic composition of waters extracted from the deeper soils (>115 cm, Figure 7) suggests integrative mixing effects of longer-term storage (Carrer et al., 2016;Sprenger et al., 2017) or limited potential for warmer precipitation to reach greater soil depths when evapotranspiration is higher.Evidence of both colder and warmer precipitation sources in these soil water samples, along with field-level spatial and annual temporal variation across sampling locations, indicates mixing processes that result from distinct patterns of water transport through soils.
Soil water contents and isotopic compositions exhibited limited sensitivity to individual precipitation or snowmelt events, indicating effective source water mixing over timescales exceeding the time between weather events and across spatial scales captured by tension lysimeters (Figures 4 and 7).Waters influenced by colder precipitation are stored in soils and can persist throughout the growing season (Figure 4a).Preservation of the colder precipitation signal in soils suggests that warmer precipitation, delivered when soil water contents are reaching saturation (Figure 3e,f), maybe more readily available for root uptake or preferentially exported to groundwater.Contributions of colder precipitation to the shallower soil (i.e., 66 cm) are evident in the steady decrease in isotope values in response to a period of consecutive low-intensity storms (e.g., April-May 2014 Figure 7b).In the subsequent season, the shallower soil also exhibited incorporation of colder precipitation, yet at higher volumes in response to higher intensity, shorter duration storms that occurred from the end of May into June (April-June Figure 7c).The initially high isotopic values of the shallower soil at the onset of each growing season suggest the influence of antecedent precipitation stored from the previous fall following plant senescence and harvest (Figure 7b,c).Limited soil water holding capacity at shallow depths and possible freezing in surface soils can restrict the additional accumulation of precipitation received in the intervening colder months, allowing a warmer precipitation isotopic signal to persist over winter.Our data illustrate how antecedent soil moisture and weather conditions are important determinants of soil water isotopic composition at the onset of a growing season, which is particularly important to managing water content in soils with limited storage capacity.Our findings also indicate that contributions of colder precipitation are important to longer-term soil water storage in these non-irrigated, cultivated systems.
Several lines of evidence speak to the heterogeneity of infiltration processes in these soils.Isotopic composition and volume sampled from individual lysimeters differed across sites and across years (Figure 7).Generally, isotopic compositions of soil water exhibit inconsistent trends across sites in the same cropped field and for an individual site across years (Figure 7).Even soils sampled at similar depths exhibited different patterns in both water content and isotopic composition (e.g., depths 117 and 119 cm in Figure 7b,c).One of the deepest soil sampling locations (119 cm) exhibited a substantial accumulation of warmer precipitation from 2013 to 2014 and then remained more consistent from 2014 to 2015 (Figure 7a-c).In contrast, the soil water isotopic composition observed at another site with a similar depth (117 cm) remained consistent from 2013 to 2014 yet exhibited a colder precipitation signal early in 2015, then showed T A B L E 3 Summary statistics for lineconditioned excess results in each process domain.active accumulation of warmer precipitation (Figure 7a-c).Assuming an equal capacity to buffer isotopic composition in these deeper soils, potential drivers of observed year-to-year changes are differences in duration and intensity of storms, changes in crop rotation, variations in snow accumulation, or atypical drying due to meteorological drought (Kleine et al., 2020;Oerter & Bowen, 2019;Sprenger et al., 2018).The variation in isotopic composition observed at similar soil depths across multiple sites within the same field may be explained by spatial variation in snow accumulation and melt or soil physical properties at the profile scale.The variation observed in soil water isotopic compositions during the growing season (June and July) suggests that the isotopic composition of water taken up by plants varies year to year according to precipitation delivery and storage patterns in the soil.An important consideration for interpreting these data is that a given point of observation captured by a lysimeter within a soil profile reflects the effects of physical soil properties that vary within the soil profile (i.e., soil structure) and dictate how mixing processes and potential preferential flow patterns manifest.Recent work demonstrates the occurrence of isotopically distinct preferential flow from soil water stored in smaller pore spaces through direct observation or modelled inference (Allen et al., 2019;Oerter & Bowen, 2017;Sprenger et al., 2015;Stumpp & Maloszewski, 2010).
Here, our observations of different isotopic and volumetric patterns at similar depths may be an expression of soil profile heterogeneities (Wanzek et al., 2018).Additionally, our data capture spatial and temporal variability within and among different soil depths that reflect a range of mixing conditions due to heterogeneities in soil at the field scale.As a result, the range of isotopic values observed in soil water exhibits increased variation and more influence of colder precipitation than the range observed in the terrace groundwater (Figure 5a).

| Mixing processes dominate in an unconfined terrace aquifer
Terrace groundwater exhibits a mixture of recharge waters from spatially variable soils that appear less influenced by colder precipitation and more influenced by evaporation than observed pore waters of overlying soils.Soil preferential flow paths facilitate the transport of isotopically different source water mixtures from surface soils to the terrace aquifer (Brooks et al., 2010;Sprenger et al., 2016).
Decreased spatial and temporal variability in terrace groundwater relative to terrace soil water indicates the integration of water received from heterogeneous terrace soils by a highly conductive aquifer (Figures 4 and 5 and Table 2; Bowen et al., 2019;Sprenger et al., 2016, Sprenger, Llorens, et al., 2019).The integrated terrace groundwater signal incorporates spatially variable processes that occur within soil profiles and across soils on the landform (Figures 4   and 7).Consistency in terrace groundwater isotopic composition over time exhibited minimal apparent influence from the seasonality of precipitation, evaporation, or variation in the timing of water movement through soils (Figures 4 and 5;Jasechko et al., 2014).Groundwater recharge occurs when soils are near saturation, which is likely to correspond with periods of high precipitation delivery to soils and, in non-irrigated agricultural systems, periods following a fallow rotation when crops are not grown for an entire season to conserve water (John et al., 2017;Sigler et al., 2020) The integrated isotopic signal in terrace groundwater suggests that while precipitation moving through soils may be isotopically variable, soil water transport and groundwater mixing were effective enough to result in groundwater isotopic values similar to the volumeweighted mean precipitation value.Assuming steady-state turnover of a well-mixed volume, the estimated 20-year groundwater residence time (Sigler et al., 2018) suggests that annual recharge received from deep percolation of soil water is a relatively small fraction (ca.5%) of the total aquifer volume.While recharge to the terrace aquifer through preferential flow paths likely occurred following wetter antecedent conditions, the limited temporal variability and minimal evaporation effects in terrace groundwater isotopic composition suggest that variation in the isotopic composition of these contributions was masked by relatively rapid mixing with the aquifer volume (Figures 4 and 6;Sprenger, Llorens, et al., 2019;Sprenger & Allen, 2020;Stumpp & Maloszewski, 2010).We note that our sampling method using tension lysimeters may bias our assessment of terrace soil water to water from larger pores or preferential pathways in the soil matrix (Orlowski et al., 2016;Sprenger et al., 2015).However, sampling with tension lysimeters does not exclusively sample the most mobile water due to the episodic availability of that water in macropores.Smaller volume lysimeter samples caused by drier soils likely reflect less mobile soil water and larger volume lysimeter samples from wetter soils likely included more mobile soil water (Figure 7).At some sites, periodically increased sample volumes corresponded with seasonally higher precipitation and increased soil water contents (e.g., 66 and 117 cm depths in Figure 7).More mobile soil water likely drives the recharge of the terrace aquifer during specific seasonal weather and antecedent soil conditions (Figure 7 and Table S1; Jasechko et al., 2014).The consistent isotopic composition observed in the terrace groundwater despite the isotopically variable recharge events suggests a relatively rapid mixing of relatively small volumes of deep percolation into a relatively large aquifer volume (Figures 5a and 7).

| Local processing and inherited signals define riparian corridors
The influence of colder precipitation on the isotopic composition of riparian water relative to terrace groundwater (Figures 4a and 5a) suggests this process domain received recharge from local colder precipitation in addition to process domains up-gradient on the landscape.
Limited isotopic variation in riparian waters and the general similarity of riparian waters to terrace groundwater suggest efficient integration of source waters by the riparian corridor (Table 2; Figures 4a and 5a).
However, some riparian groundwaters reflected the increased accumulation of colder precipitation relative to riparian surface water, terrace groundwater, and stream water (Figures 5a and S5, e.g., April and May 2021 in Figure S6).Previous studies found that a lack of evapotranspiration in winter months promotes a shallower groundwater table that buffers riparian soil temperatures (Burt et al., 2002;Smerdon & Mendoza, 2010).Increased probability of non-frozen soils in addition to increased snow catch allow colder precipitation to accumulate in riparian soil and groundwaters.Although we did not evaluate riparian soil temperature and water content data during the period of this study, we anticipate similar trends to observations in terrace soils where water accumulates in winter and spring (Figure 3).These infiltration dynamics allow for increased contributions of colder precipitation to riparian groundwater.
Riparian waters undergo local processing in addition to inherited signals.Riparian water isotopic compositions that deviated from the majority and the LMWL (ca.10% of riparian samples) provide evidence of local processes that promote evaporation (Figures 5b, 6c, and S5;Gammons et al., 2006;Sharp, 2017).Increased fractionation in riparian water was especially evident in riparian surface water where reduced velocity and increased exposure to the atmosphere facilitate evaporation (Figure S5).Riparian corridors encompass flow paths with strong variation in transit and residence times, where longer flow paths promote biogeochemical processes that can decrease nitrate concentrations in water leaving the system and can contribute to N 2 O emissions (Chapin et al., 2011;Cheng et al., 2020;Harms & Grimm, 2008;Hedin et al., 1998;Hill, 2019;Hill et al., 2014;Sigler et al., 2022).Areas of more pronounced evaporation within the riparian corridor may reflect increased water residence times and hence exert control on biogeochemical function, reflect the spatial heterogeneity of these ecosystems, and indicate the influence of riparian structure on hydrologic processes.
Stream water isotopic compositions indicate a mixture of upgradient channel water, terrace groundwater, and riparian water.The stream water reflects a relatively consistent mixture of inherited components, with limited evidence of local processes (Figures 5 and 6).
Instances of increased isotopic signals in the stream water relative to riparian groundwater suggest that colder precipitation (e.g., snowmelt) delivered to the riparian may flush warmer season water into the stream channel (Table 2, Figures 5 and 6d).Local evaporative effects on stream water may also result in isotopic signals distinct from upgradient domains.The stream water exhibited the most consistent isotopic signal and mixture of source waters observed across the process domains, which supports our conceptualization of this terrace-riparian-stream continuum as a gaining system where downgradient domains integrate the effects of up-gradient processes by way of hydrologic connection (Hill, 2019;Vidon & Hill, 2004).

| CONCLUSIONS
Our data provide a novel analysis of water isotopic compositions Our study area is located within the Judith River Watershed (JRW) of central Montana (MT), USA.The Judith River is a tributary of the Missouri River on the western edge of the NGP (Figure 2).Depositional terraces and wider floodplains in the lower-elevation areas of the watershed are extensively cultivated primarily for non-irrigated cereal crop production.Narrower stream corridors dissecting the terraces are often used as pasture for cattle production to increase producer diversification (U.S.Department of Agriculture, 2023).Dominant crop types are spring wheat, winter wheat, and barley with increasing incorporation of pulse crops in cropping system rotations to replace fallow (U.S.Department of Agriculture, 2023; John Snow water equivalent was recorded every 15 minutes with a snow scale (Sommer SSG-2) installed in February 2020 at the Montana State University (MSU) Central Agricultural Research Center (CARC) located on the Moccasin terrace (Figure 2).Soil volumetric water content and temperature were measured every 15 min at depths of 10, 20, 50, and 90-100 cm by Montana Climate Office (2023) and USDA NRCS (see Supporting Information; U.S. Department of Agriculture, 2022).Precipitation volume, air temperature, and relative humidity were measured at 1-h intervals by the U.S. Bureau of Reclamation (2023).

F
I G U R E 3 Precipitation and soil climate data by 2021 and 2022 water years at the Central Agricultural Research Center, Moccasin, MT.Data include (a, b) daily precipitation, snow water equivalent (SWE), and atmospheric temperature, (c, d) soil temperature at 20, 50, and 90 cm depths and (e, f) percent soil water content at 20, 50, and 90 cm depths (U.S. Bureau of Reclamation, 2023; Montana Climate Office, 2023).Soil volumetric water content is invalid at soil temperatures below 0 C due to changes in geophysical properties caused by state change from water to ice.T A B L E 1 Sample type, period of sample collection and total count.

δ 2 H
and calculated lc-excess values while accounting for repeated measures of the same sites over time (R software version 4.2.1 package 'nlme' version 3.1-160).We assigned 'water-year' as a fixed effect (i.e., independent variable) in the statistical model because δ 18 O and δ 2 H values in environmental waters are a function of water year-specific climatic conditions.We assigned a sampling campaign within the year as a random effect to account for the increased probability that results from samples collected closer in time were more similar than samples collected further apart in time.The respective assignments of 'water-year' and 'sample campaign' allowed us to account for autocorrelation in residual errors of the data.We also included individual sites as a random effect in the model to account for the increased probability that results from samples collected at the same site were more similar than samples collected from other sites.Accounting for site-to-site variation greatly improved the Akaike and Bayesian information criteria when compared to the model not accounting for individual sites (ANOVA p-value <0.001; .1 | Terrace soil water content and temperature Temperatures and water contents were most variable in the top 20 cm of terrace soils near Moccasin, MT (Figures 3 and S1).The effects of evapotranspiration were limited at soil depths near 50 cm, and soil water contents at these depths reach a minimum of ca.0.25 m 3 m À3 , consistent with field capacities typical of the clay loam to coarse sandy loam soils present across our study area(Montana Climate Office, 2023;Sigler et al., 2020).Soil water content at 90 cm reached a minimum less than ca.0.20 m 3 m À3 , consistent with field capacities of coarser substrate textures at depth (Figure3).Soils exhibited peak water content values at all depths during May, with the greatest variation occurring at shallower depths (Figure3).During the spring season, maximum water contents reached ca.0.40 m 3 m À3 at all depths, suggesting a nearly saturated state assuming porosity consistent with clay loam to sandy clay loam textures (U.S.Department of Agriculture, 2022).

F
I G U R E 4 Water isotope and line-conditioned excess values by day of the water year.Values of δ 2 H (a) and lc excess (b) in precipitation (rain, snow, and monthly composite in collaboration with the Montana Precipitation Isotope Network, MTPIN), terrace soil water, terrace groundwater, riparian water, and stream water with multiple years (2013-2022) represented.The dashed line (a) indicates the multi-year volume-weighted mean precipitation value, with uncertainty (dotted lines) reflecting variability in annual precipitation over the AgriMet record for 2012-2022 (U.S. Bureau of Reclamation, 2023).The solid line (b) indicates the local meteoric water line, with uncertainty (dashed lines) reflecting the 95% confidence interval derived from the Monte Carlo analysis (Figure S4, see methods).F I G U R E 5 Values of δ 2 H and line-conditioned excess by process domain.Distribution of δ 2 H (a) and lc-excess (b) values in four process domains across the terrace including terrace soils (orange), terrace groundwater (green), riparian surface and groundwater (magenta), and stream water (red).Mean values and 95% frequentist confidence intervals are indicated by black bars for each process domain.Lowercase letters indicate significant differences determined by Tukey HSD post hoc test of the impact of landscape process domain on δ 2 H and lc-excess after accounting for the water year, campaign within a year, and individual site.The dashed line on (a) indicates the multi-year volume-weighted mean precipitation value, with uncertainty (dotted lines) reflecting variance in annual precipitation over the AgriMet record for 2012-2022 (U.S. Bureau of Reclamation, 2023).The cluster of more negative values in terrace groundwater reflects results from the Montana Department of Agriculture M-1 monitoring well during a rapid snowmelt runoff episode in 2014.The solid-grey, zeroline (b) expresses the local meteoric water line (Figure S4).The 95% confidence intervals (dashed grey) are derived from the Monte Carlo uncertainty analysis (Figure S4, see methods).Values below the lower bound of the confidence interval indicate an evaporative loss signal either from within that domain or via inmixing from up-gradient.

6
Water isotope values in four process domains.Values of δ 18 O and δ 2 H from (a) terrace soil water (2014-2016); (b) terrace groundwater sampled from wells (2013-2017) and springs (2013-2022); (c) riparian water sampled from spring channels, ponds, and shallow groundwater (2020-2022); and (d) stream water from two channels (2012-2022) relative to the modelled volume weighted annual mean value (black point with error bars based on variability in annual precipitation over the AgriMet record for 2012-2022 (U.S. Bureau of Reclamation, 2023).The local meteoric water line (LMWL) and 95% confidence interval are derived from precipitation results (Figure S4; see methods).Points falling outside the confidence interval indicate transformation processes and mixing of inherited transformed waters.

F
I G U R E 7 Precipitation amounts and terrace soil water δ 2 H values across years.Daily precipitation (mm) from USBR AgriMet station: MWSM is in grey bars along the top of each panel for water years 2013 (a), 2014 (b), 2015 (c) and 2016 (d).Soil water was collected via tension lysimeter across four sampling sites (denoted by color, legend above panel a) representing a range of depths of fine soil material across the study area.Point size is proportional to the volume of the sample collected which ranges from 1.5 to 181 mL.The dashed line indicates the multi-year volume-weighted mean precipitation value, with uncertainty (dotted lines) reflecting variance in annual precipitation over the AgriMet record for 2012-2022 (U.S. Bureau of Reclamation, 2023).Only sites with four or more samples collected in each of two or more years are included in this figure.
across a non-irrigated agricultural setting.Observed δ 18 O and δ 2 H values in environmental waters across process domains show that soils play an integral role in mixing seasonal precipitation and exhibit spatial variation of mixing and throughflow processes within and across soil profiles.Our results highlight the role of soils in processing water that enters the landscape and indicates that colder precipitation influences water stored in soils across seasonal to annual and longer timescales.The persistence of colder precipitation in soils suggests these source waters contribute to crop production and nutrient processing.We found that shallow riparian aquifers likely receive local contributions from colder precipitation and act as efficient integrators of source waters.Areas within the riparian corridor where increased water retention coupled with exposure to atmospheric forces impart localized signatures of riparian isotopic compositions.Our data show that stream channel waters reflect a mixture of riparian water and terrace groundwater, which suggests the influence of direct inflows of terrace groundwater that undergo minimal riparian processing.Streams exhibiting limited riparian influence reflect corridors that are less effective at mitigating the export of nutrient loads from agricultural systems.Additional efforts to quantify how soil water storage heterogeneities mediate variable precipitation patterns will further inform the potential consequences of warmer winters with more rain and less snow across both non-irrigated agricultural production and grassland ecosystems.ACKNOWLEDGEMENTS The authors thank the Consortium for Research on Environmental Water Systems (NSF EPSCoR Cooperative Agreement OIA-1757351), the MSU Nielsen Pedology Graduate Research Fellowship (2021funding this work.Any opinions, findings, conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation or the U.S. Department of Agriculture.The authors appreciate the support and data from the Montana Climate Office, the US Department of Agriculture Natural Resources Conservation Service, the U.S. Bureau of Reclamation, and the Montana Department of Agriculture.A special thank you to Dr. Ann Marie Reinhold for her help with building the conceptual framework.The authors thank Dr. Toby Koffman, Kevin Hyde, Sam Leuthold, Simon Fordyce, and Skye Keeshin for their help with data collection, analysis, and technical assistance.The authors sincerely thank two anonymous reviewers for their time and constructive feedback.
for δ 18 O and 0.9‰ for δ 2 H were used in a Monte Carlo propagation of uncertainty to the LMWL, following the approach ofLeuthold et al.
(Bowen, 2022)nd δ 2 H values were determined for precipitation over the study period based on annual volume-weighted mean water values for the years 2012-2022.A sinusoidal model was fit to measure δ 18 O and δ 2 H values in precipitation to allow interpolation of the isotopic composition of precipitation throughout the year.We assessed whether the model fit to our precipitation isotopic data was reasonable by comparing it to a sinusoidal model fit to monthly predictions of isotopic values obtained from the Water isotopes database(Bowen, 2022)for a central location in our study area (see Supporting Information 2.3; Figure . The similarity of the terrace groundwater isotopic composition to the volume-weighted mean pre-