It is not yet clear how this legacy nitrogen may respond to changing hydrologic regimes and variations in AgMAR practices, and more importantly, if flooding agricultural sites is enhancing nitrate transport to the groundwater or attenuating it by supporting in situ denitrification. Denitrification rates in the subsurface have been reported to vary as a function of carbon and oxygen concentrations, as well as other environmental factors . While total soil organic carbon typically declines with depth , dissolved organic carbon can be readily transported by water lost from the root zone to deeper layers and can therefore be available to act as an electron donor for denitrification . Oxygen concentration in the vadose zone is maintained by advective and diffusive transport, while oxygen consumption occurs via microbial metabolic demand and/or abiotic chemical reactions . The effects of drying and wetting cycles on oxygen concentrations in the deep subsurface are not well documented. However, in 1 meter column experiments, there is some evidence that O2 consumption proceeds rapidly as saturation increases and rebounds quickly during dry periods . These variations in oxygen concentration can influence N cycling and thus, transport to groundwater. Variability in nitrate concentration has also been linked to heterogeneous subsurface properties, rainfall events, seasonality of flow and other local geochemical conditions across a diversity of settings However, a gap currently exists in quantifying N attenuation and transport from agriculturally intensive regions with a “deep” vadose zone while accounting for the many competing N cycle reactions and transformations, as impacted by different hydrological regimes imposed under AgMAR. The application of AgMAR itself can vary in terms of the hydraulic loading and rates used, as well as the duration between flood water applications. These can in turn affect water retention times, O2 availability, consumption of electron donors and consequently, denitrification rates . For example, denitrification rates were found to increase with increased hydraulic loading and with shorter times between flood applications within the vadose zone of a rapid infiltration basin system used for disposing of treated wastewater . In shallow, sandy soils, high flow rates – above an infiltration threshold – were negatively correlated with denitrification rates,plant pots with drainage suggesting that an optimum infiltration rate exists for a given sediment stratigraphy to maximize NO3 – reduction .

Given the immense stratigraphic heterogeneity in alluvial basins, such as in California’s Central Valley, a range of optimum infiltration rates may exist with implications for managing AgMAR differently based on the geologic setting of the intended site. Therefore, the objectives of this study are to: a) understand the effects of varying stratigraphy and hydrologic regimes on denitrification rates, and b) identify AgMAR management scenarios that increase denitrification rates, such that the potential for N leaching to groundwater is decreased. Herein, we focus on an agricultural field site in Modesto, California located within the Central Valley of California, which is responsible for California’s $46 billion-dollar agricultural economy . The field site typifies the deep vadose zones prevalent in this region, which are characterized by heterogenous layered alluvial sediments intercalated with discontinuous buried clay and carbon rich paleosols . These discontinuous, layered features, especially the paleosols and areas of preferential flow, are typically associated with enhanced biogeochemical activity, higher carbon content and availability of metabolic substrates such as nitrogen . These regions respond to and change depending on environmental conditions such as water content and oxygen concentration in situ that are influenced by the hydrologic regime at the surface and may be important for NO3 – attenuation and reduction prior to reaching the water table. Therefore, this study considers varying hydrologic regimes and stratigraphic variations that are prevalent in the region. More specifically, at the Modesto field site , large amounts of legacy N already reside in the vadose zone, while N fertilizer application and irrigation occurs throughout the spring and summer months. AgMAR, if implemented, occurs during the winter months as water becomes available. Therefore, we focus here on quantifying the effects of AgMAR on N cycling in the deep vadose zone and implications for NO3 – contamination of groundwater at this characteristic agricultural field site. We also investigate the specific AgMAR application rates that would increase the effectiveness of in situ denitrification under different stratigraphic configurations through the development and testing of a reactive transport model. We believe such an analysis provides important insights for the successful application of AgMAR strategies aimed at improving groundwater storage in a changing climate.

Reactive transport models can be beneficial tools to elucidating N fate and transport in deep vadose zone environments. Herein, we develop a comprehensive reaction network incorporating the major processes impacting N transport and attenuation, such as aqueous complexation, mineral precipitation and dissolution, and microbially mediated redox reactions. While using the same reaction network, we implement several numerical scenarios to quantify the range of denitrification rates possible under different AgMAR implementation strategies and stratigraphic configurations . For the latter, we used four different stratigraphic configurations with a low permeability layer on top including i) two homogeneous textural profiles, ii) a sand stratigraphy with a discontinuous silt band, iii) a silt stratigraphy with a discontinuous sand band, and iv) a complex stratigraphy more representative of the field conditions investigated by electrical resistance tomography . The top layer served two purposes, one, it allowed the net infiltration rate to be calibrated to match measured average field infiltration rates of 0.17 cm/hr and two, it represented the expected increase in sediment uniformity expected in ploughed or tilled layers in agricultural settings. While, the impact of the top layer resulted in water being delivered more slowly to the heterogenous sediments below, varying rates of percolation occurred after reaching below the more homogenous layer allowing us to examine the effects of heterogeneity on nitrate transport and fate in the vadose zone. For each stratigraphy, we further varied the frequency and duration of water per application to investigate the impact of different AgMAR implementations that are similar to recent field trials conducted throughout the state . In addition, we tested the effect of antecedent moisture conditions on N biogeochemistry within the more complex stratigraphy by setting the model with a wetter initial moisture profile. Overall, a set of 18 simulation experiments were used to isolate and understand the contribution of different AgMAR strategies to enhance or decrease denitrification rates in deep vadose zone environments with homogeneous and banded configurations. A detailed model setup and numerical implementation is provided in Section 2.3. Although our reactive transport analysis was guided by a particular field site that is classified as a “Medium to Good” site for MAR , our aim was not to replicate site conditions in its entirety, but rather to enhance our understanding of how hereogeneity might impact nitrogen transport and fate under MAR.

The study site is an almond orchard located in California’s Central Valley, southwest of Modesto, and north of the Tuolumne River . The surface soil is classified as a Dinuba fine sandy loam . The site is characterized by a Mediterranean climate, with wet winters and hot, dry summers. Average annual temperature and total annual precipitation are 17.5° C and 335 mm, respectively. As suggested above, the vadose zone typifies the valley with contrasting layered sequences of granitic alluvial sedimentary deposits consisting of predominantly silt loams and sandy loams. We therefore use these textures to design our modeled stratigraphic configurations with and without banded layers. The groundwater table in the study area typically occurs around 15 m below ground surface. Soil properties including percent sand, silt, clay, total N, total C, and pH are shown in Table 1.To specifically characterize the textural layers and subsurface heterogeneity at our site, we used electrical resistivity tomography . ERT profiles were generated along a 150 m transect to 20 m depth prior to flooding to quantify subsurface heterogeneity while the subsurface was relatively dry . Further,plastic plants pots to validate the texture profiles generated by the ERT data, a set of six cores were taken along the transect of the ERT line down to nine meters with a Geoprobe push-drill system . The first meter of the core was sampled every 25 cm. Thereafter, cores were sampled based on stratigraphy as determined by changes in color or texture. The ERT profiles were used to develop the stratigraphic modeling scenarios and the coring guided the specification of the hydraulic parameters. Redoximorphic features were noted throughout the cores. centrifuge tubes with 40 mL of 0.5% sodium phosphate and shaken overnight . Samples were hand shaken immediately before a 2.5 mL aliquot was taken 11 seconds and 1 hour and 51 minutes , respectively after shaking and placed in a pre-weighed tin. Tins were oven dried at 105°C overnight and Statistical analysis was used to help guide the development of the geochemical reaction network. First, correlation analysis was used to inform the choices of primary geochemical species on the basis of the strength of their relationship with N2O. Second, on the basis of cluster analysis, stratigraphic configurations with different textural classes were developed. In particular, a Spearman’s rank correlation was conducted on the dataset including several physical and geochemical measurements collected on the soil cores. Specific variables included pH, N2O, NO3 – , NH4 + , DOC, Fe, Mn, S, total C, percent sand, silt, and clay, and depth. Variables were standardized using the median and mean absolute distance because most variables were found to be non-normally distributed based on the Kolmogorov-Smirnov test. To further understand how the data grouped, a cluster analysis was conducted using the partitioning around medoids method for the same set of variables. Interestingly, data were found to group according to textural classes and depth, which provides a mechanism to develop the modeling strategy around these textural profiles.Several scenarios were developed based on the soil textures identified in cores and the ERT profiles to provide insights into the effect of stratigraphic heterogeneity and AgMAR management strategies on NO3 – cycling in the deep subsurface, as described in section 2 above.

The five stratigraphies modeled in this study are shown in Figure 1. The limiting layer in the ERT scenario spans 187 to 234 cm-bgs based on field core observations. For each lithologic profile, three AgMAR management strategies were imposed at the top boundary between 20 m and 150 m of each modeled profile . For each AgMAR management strategy, the same overall amount of water was applied, but the frequency, duration between flooding events, and amount of water applied in each flooding event varied : a total of 68 cm of water was applied either all at once , in increments of 17 cm once a week for four weeks , in increments of 17 cm twice a week for two weeks , and all three scenarios with an initially wetter moisture profile . Note, that for all scenarios, the same reactions were considered, the water table was maintained at 15 m, and temperature was fixed across depths at 18°C, the mean air temperature for January to February in Modesto. For all scenarios, the modeling domain consists of a two-dimensional 20-meter deep vertical cross-section extending laterally 2,190 m and including a 190 m wide zone of interest located at its center, thus distant from lateral boundaries on each side by 1,000 m to avoid boundary effects. The zone of interest was discretized using a total of 532 grid blocks with a uniform grid spacing of 1 m along the horizontal axis, and a vertical grid spacing of 0.02 m in the unsaturated zone increasing with depth to 1 m in the saturated zone. A maximum time step of 1 day was specified for all simulated scenarios, although the actual time step was limited by specifying a Courant Number of 0.5, typically resulting in much smaller time steps during early stages of flooding. Before each flooding simulation, the model was run first to hydrologic steady state conditions including the effect of average rainfall . The water table was set at a depth of 15 m by specifying a constant pressure at the bottom model boundary , and the model side boundaries were set to no-flow conditions. Under these hydrologic conditions, the model was then run for a 100-yr time period including biogeochemical reactions and fixed atmospheric conditions of O2 and CO2 partial pressures at the top boundary, a period after which essentially steady biogeochemical conditions were achieved, including the development of progressively reducing conditions with depth representative of field conditions.