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, 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,blueberry pot 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. The clusters were thus automatically grouped by decreasing levels of denitrification and microbial activity. 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 HSreduction 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,nursery pots with silt loams having higher water holding capacity and lower pore gas velocities compared to sandy loams.
Consequently, lower NO3 – concentration and lower NO3 – :Clratio are predicted in the silty loam vadose zone as compared to the sandy loam column . It is interesting to note that while greater NO3 – loss and denitrification are predicted for the silty loam vadose zone, carbon concentration associated with the shallow vadose zone are comparatively lower than for the sandy loam column. Moreover, the calculated pH is lower and iron concentrations are higher in the silt loam profile below the top meter when compared to the same depths within the sandy loam column . This suggests that chemolithoautotrophic reactions could be more important for these finer textured sediments. While both heterotrophic and chemolithoautotrophic reactions would be expected to result in a pH decrease , the greater decline in pH and concomitant increase in Fe+3 concentration suggests the importance of Fe and S redox cycling associated with the chemolithoautotrophic reactions in silty loam sediments . Evolving from these steady state conditions, scenario S1 suggests that denitrification is enhanced as floodwater infiltrates into the silt loam column. Model results indicate that saturation increases to 80% from 1 to 4 m depths and O2 decreases from 2.1 x 10-4 mol L-1 to 1.7 x 10-4 mol L -1 , resulting in 43% of the NO3 – being denitrified for this scenario . In comparison to the homogeneous profiles, the sandy loam with silt loam channel stratigraphy has higher calculated water contents and slightly lower O2 concentration within and surrounding the silt loam channel than the homogenous sandy loam column under steady state conditions . Calculated NO3 – concentrations are also similar between the homogenous sandy loam column and SaSi case, except for within and below the silt loam channel where lower NO3 – concentration was predicted . For scenario S1, water content for the SaSi case increased in a manner similar to the homogenous sandy loam, except for within the silt loam channel, which increased from 60 to 81%. Figure 4 further demonstrates that the infiltrating floodwater resulted in an increase in NO3 – concentration between 1 and 3 m within the sandy loam textured soil, but a decrease elsewhere. Within the channel itself , lower nitrate and NO3 – :Clratio are predicted, suggesting higher rates of denitrification . Overall, the model results indicate that an average of 37% of the NO3 – concentration is denitrified in the SaSi case 60 days after flooding, with 35% denitrification occurring in the sandy loam matrix and 40% occurring within the silt loam channel. This suggests that the silt loam channel acts as a denitrification hotspot. Furthermore, the silt loam channel has lower carbon and higher Fe+3 concentrations similar to the homogenous silt loam column again suggesting the importance of both heterotrophic and chemolithoautotrophic denitrification in these finer textured sediments. In comparison to the SaSi case, calculated water saturation and O2 profiles were markedly different between the homogenous silt loam column and the silt loam with sandy loam channel under steady state conditions . In particular, the sandy loam channel has lower calculated water content than the homogenous silt loam column . Further, greater gas flux within the channel resulted in 11-19% higher O2 concentration that penetrated deeper into the vadose zone as compared to the homogeneously textured column. NO3 – concentration are also estimated to penetrate deeper into the vadose zone in the SiSa case due to the high permeability of the sandy loam channel . While carbon concentration also penetrated deeper in the vadose zone in the SiSa case, higher calculated O2 concentration did not allow for comparable rates of denitrification below 1 m in this case as observed in the homogenous silt loam profile. This is further confirmed by the lower NO3 – :Clratio, which indicates that transport processes dominate biogeochemical fluxes within this column . With scenario S1, the calculated water content increased to 48% saturation while the O2 concentration remained the same within the channel. The high permeability channel allowed for NO3 – to move faster and deeper into the vadose zone. Overall, calculated denitrification was lower in the SiSa case as compared to the homogeneous textured column. In the simplified ERT stratigraphy, similar patterns were observed such that high permeability channels transported water, O2, and NO3 – faster and deeper into the subsurface than low permeability regions . As a result, concentration profiles showed significant variability across the modeled domain even under steady state conditions. For example, the calculated O2 and NO3 – concentrations are an order of magnitude lower in the shallow vadose zone below the limiting layer than within the preferential flow channel. Higher NO3 – :Clratio within the channel further confirms that preferential flow paths transport higher quantities of dissolved aqueous species without their being impacted by other processes such as denitrification . Dissolved carbon in particular is predicted to have a lower concentration in the preferential flow channel and the matrix surrounding the channel than below the limiting layer. In contrast, the Fe+2 concentration is estimated to be higher in the matrix surrounding the preferential flow channel and below the limiting layer . For scenario S1, model results indicate that NO3 – moved through the preferential flow path faster and deeper into the profile, while the limiting layer acts as a denitrification barrier as evidenced by the decrease in NO3 – :Clratio. The highest denitrification was estimated to occur in the matrix adjacent to the preferential flow channel , followed by intermediate nitrate reduction below the limiting layer and far away from the channel , while the lowest denitrification was estimated to occur within the channel itself .