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,bucket flower 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, 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. 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. For these simulations, the concentrations of dissolved species in background precipitation and in groundwater at the bottom model boundary were fixed, with compositions described in Table 2 to yield similar vertically distributed NO3 – concentrations as were measured in the soil cores. Flooding scenarios were then started from the initially steady flow and biogeochemical conditions developed as described above and run for 60 days. For these simulations, a free surface boundary was implemented for scenario S1 where 68 cm of water was applied all at once. In contrast, a specified flux boundary condition was imposed for the scenarios S2-S3, where floodwater applications were broken up over a week. The flood water composition is discussed in Section 2.3.5.
The groundwater composition was taken from analyses reported by Landon and Belitz for a groundwater well located near our study site. For simplicity,cut flower bucket the background recharge from rainfall was assumed to have the same composition as groundwater except that it was re-equilibrated under atmospheric O2 and CO2 conditions prior to infiltration. In addition, the concentrations of N species in the background recharge were set to values determined from our own analyses of N at the top of soil cores. The composition of the flood water was set to that of the background precipitation diluted by a factor of 100 for most constituents except for Cl-1 . Ratios of NO3 – to Cl-1 were used to trace the difference between dilution and denitrification effects on NO3 – . Denitrification and N2O production were simulated as aqueous kinetic reactions coupled to the fate of pH, CO2, Fe, S, NO3 – , and NH4 + based on the Spearman correlation analyses discussed above . Apart from pH and nitrate species, Fe and S have been linked to denitrification through chemolithoautotrophic pathways in addition to heterotrophic denitrification , and are therefore included in our reaction network. Heterotrophic denitrification of NO3 – to N2 was represented via a two-step reduction process of NO3 – to nitrite and NO2 – to dinitrogen . Additionally, chemolithoautotrophic reduction of NO3 – to N2 with Fe and bisulfide as electron donors were implemented. Further, dissolved organic carbon was observed throughout the nine-meter profile at our field site, and CO2 and N2O profiles showed strong correlation . Therefore, DOC degradation was simulated using Monod kinetics, although individual DOC components were not simulated consistent with other modeling studies . In particular, we considered a single solid phase of cellulose in equilibrium with acetate as the source of DOC. Parameters for cellulose dissolution were calibrated using the total organic carbon concentrations obtained for each cluster. Biodegradation of acetate was coupled to multiple terminal electron acceptors, including NO3 – , Fe and SO4 2- which follow the hierarchical sequence of reduction potential of each constituent implemented by using inhibition terms that impede lower energy-yielding reactions when the higher energy yielding electron acceptors are present. These microbially mediated reactions and their kinetic rate parameters are shown in Table 5. Rates for denitrification were calibrated using the results from the acetylene inhibition assays as described above. Enzymes involved in denitrification include nitrate reductase, nitrite reductase and nitrous oxide reductase. To remain conservative in our estimates, we chose values typical for oxygen inhibition of nitrous oxide reductase L -1, the most sensitive to oxygen of the enzymes . Spearman rank correlation indicated that pH, DOC, S, NO3 – , and Fe exhibit significant correlation with N2O and therefore, these geochemical species were included in the reaction network. Cluster analysis was used to further detect natural groupings in the soil data based on physio-chemical characteristics, textural classes and the total dataset. Cluster analysis revealed three clusters representing distinct depth associated textural classes with varying levels of substrates and biogeochemical activity. Table 5 shows the median and range for N2O, CO2, NO3 – -N, Fe, S and total organic C for each of the clusters. The first cluster is dominated by sandy loams within the top meter with highest median values of total N2O, total CO2, NO3 – -N, Fe, and total organic C concentrations, indicative of greatest microbial activity and denitrification potential. The second cluster is dominated by silt loams below one meter and had average values of total N2O, total CO2, NO3 – -N, Fe, and total organic C concentrations when compared to the other groups. The third group is dominated by sands and sandy loams below 1 meter and had the lowest median values of total N2O, total CO2, NO3 – -N, Fe, and total organic C concentrations amongst all groups. While most concentrations followed a decreasing concentration trend from cluster 1 to 3, the highest median values of S were associated with cluster 2. Liquid saturation profiles and concentration of key aqueous species predicted at different times for the homogeneous sandy loam column are shown in Figure A1. The sandy loam vadose zone is computed to be 32% saturated with near atmospheric concentrations of O2. As a result of oxic conditions, model results demonstrate significant residual NO3 – concentration within the vadose zone . Evolving from these conditions, Figure A1d shows that with flooding scenario S1, water reaches depths of 490 cm-bgs and saturation levels reach 40% in the sandy loam column. Deeper in the column, lower saturation and only small decreases in O2 concentration are predicted . Calculated concentration profiles show that O2 introduced with the infiltrating water is persistent at shallow depths down to 100 cm-bgs, below which O2 declines slightly as floodwater moves below this zone. Model results further indicate higher NO3 – reduction in the shallow vadose zone including the root zone with 35% of NO3 – being denitrified . Overall, this scenario results in NO3 – concentration persisting at depth. While other redox reactions, such as iron reduction and HS reduction of NO3 – to N2, may be important, conditions needed to induce these reactions were not realized in the sandy loam vadose zone due to the high pore gas velocities of the homogenous sandy loam allowing for large amounts of O2 to penetrate the profile from the incoming oxygenated water. In comparison to the homogenous sandy loam column, the predicted water content is higher and O2 concentration is 53% lower in the vadose zone of the homogenous silt loam column at steady state . This result is expected because of the difference in porosity, with silt loams having higher water holding capacity and lower pore gas velocities compared to sandy loams.