Hydrobiogeochemical Function of Soil Based Onsite Wastewater Treatment Systems: Insights from High-Resolution O 2 Imaging

: Innovative/alternative onsite wastewater treatment systems (I/A OWTS) are increasingly installed to protect groundwater and surface water from excess nitrogen. Improved understanding of the relationships between hydraulic and biogeochemical processes can help to make more informed decisions on OWTS design to maximize their performance and resilience. In this study, for the first time, oxygen (O 2 ) imaging was used to investigate redox dynamics within a soil-based I/A OWTS during the dosed application of artificial wastewater. Mesocosms consisted of an unsaturated sand matrix (nitrifying layer) above a water-saturated sand-sawdust matrix (denitrifying layer). Mesocosms consistently removed > 85% of total inorganic nitrogen, with estimated nitrification and denitrification rates of 8.9 – 14.3 and 13.8 – 14.3 nmol N h − 1 cm − 3 , respectively. Dosed application of artificial wastewater resulted in intermittent oxygenation of the upper horizon of the sand-sawdust layer due to the pulsed release of oxic water from the overlying unsaturated sand layer. Hydraulic disturbance, simulated by effluent port closure, inhibited this pulsed water release and resulted in the formation of water-saturated pockets in the sand layer, which rapidly turned anoxic. Water ponding at the surface and the loss of interconnected air-filled pore spaces inhibited diffusive O 2 supply from the atmosphere, which caused O 2 partial pressure to decline at rates between 4 . 2 and 6 . 7 matm h − 1 throughout the sand layer, most likely as a result of nitrification. When the hydraulic function was restored, significant volumes of the sand layer remained water-saturated, and anoxic regions intensified, indicative of capillary hysteresis associated with wetting and drainage at low suction pressures during the hydraulic disturbance and sustained increased water content following the disturbance. Consequently, the nitrifying performance of the system declined, with effluent NH þ 4 concentrations linearly increasing with decreasing areas of the oxic unsaturated sand. Overall, this study highlights the critical role of unsaturated matrix composition and compaction for the system performance and resilience, and illustrates the potential of chemical imaging to coregister hydrological and biogeochemical processes in soil-based I/A OWTS.


Introduction
Nitrogen pollution of groundwater and surface water is a global problem, posing serious human health risks and environmental concerns. Consuming nitrate-contaminated water has been linked to methemoglobinemia in infants (Fewtrell 2004) and cancer in adults (Ward et al. 1996;Weyer et al. 2001). Environmental consequences of surface water nitrogen pollution include eutrophication and increased algae biomass in fresh and coastal water bodies (Bowen and Valiela 2001), tendency toward harmful algae blooms (Gobler et al. 2012), coastal hypoxia (Bergondo et al. 2005;Diaz and Rosenberg 2008), and marine organism mortality (Gobler et al. 2008). In densely populated areas with infrastructure that does not provide tertiary treatment, domestic wastewater is often a dominant source of nitrogen (Andreoli et al. 1979;Bowen and Valiela 2001;Lloyd 2014). Over 60 million US residents rely on conventional onsite septic systems that provide primary and secondary treatment but are not designed to remove nitrogen and thus are a major cause of impaired groundwater quality (EPA 2002;Harman et al. 1996;Withers et al. 2014). Practical, cost-effective, and reliable onsite wastewater treatment solutions are needed to reduce the nitrogen input to groundwater and surface water.
Nitrogen removing biofilters (NRBs) are examples of innovative/ alternative onsite wastewater treatment systems (I/A OWTS) that can help reduce the residential contribution of ammonium (NH þ 4 ) that leads to groundwater nitrate-nitrite (NO x ) contamination (Anderson and Hirst 2015;Waugh et al. 2017). Following primary treatment in a septic tank, septic tank effluent (STE) is distributed evenly by dosed application through pressurized pipes across the NRB surface. Water flows by gravity downward through an unsaturated sand layer where oxic conditions promote microbially mediated nitrification (Langlois et al. 2020). Water subsequently flows through a sand-woodchip or sand-sawdust layer where oxygen (O 2 ) is rapidly consumed through mineralization of the organic substrate (e.g., sawdust or woodchips) by aerobic microorganisms. The sand-sawdust layer is designed to maintain water saturation by either using a plastic liner or by mixing a fine-textured soil with the sand-sawdust to retain moisture. The resultant anoxic conditions then allow the mineralization to proceed anaerobically via denitrification, leading to the removal of fixed nitrogen from the system by converting the dissolved NO x species into dinitrogen gas (N 2 ) (Blowes et al. 1995;Robertson and Cherry 1995). 1 NRBs and other soil-based I/A OWTS are increasingly being implemented across the US and elsewhere (Anderson and Hirst 2015;Heufelder et al. 2008;Lopez-Ponnada et al. 2017;Rask et al. 2013;Waugh et al. 2017), but the interplay between hydrological and biogeochemical processes during dosed STE application, or in response to hydrological and chemical disturbances, is complex and not yet fully understood (Beal et al. 2005). I/A OWTS performance is typically evaluated by monitoring influent and effluent water chemistry; however, hydrobiogeochemical processes that take place within the nitrifying and denitrifying matrices during normal operation and system failure are not well understood. Such insights can help optimize the design and operation of I/A OWTS (Lancellotti et al. 2017;Lopez-Ponnada et al. 2017). For example, STE is typically applied to soil-based OWTS through pressurized pipes in doses (Costa et al. 2002;Otis et al. 1980), but it is unclear how the dosing pattern affects preferential water flow, water hydraulic retention time (HRT) in the unsaturated sand layer, geochemical dynamics at the redox transition, and, ultimately, nitrogen removal. Fine silt and clay particles are sometimes added to the unsaturated sand layer to increase the HRT (MacQuarrie et al. 2001;Otis et al. 1980;Robertson and Cherry 1995) and to serve as an alkalinity source to counteract alkalinity consumption by nitrification (Wilhelm et al. 1994). However, these fine particles also increase capillary forcing and thus can result in the formation of poorly vented anoxic regions or in ponding above the infiltration surface (Siegrist et al. 2000). These hydraulic conditions can reduce the volume of interconnected air-filled pore spaces and limit diffusive O 2 supply, decreasing the nitrification efficiency of the sand (Siegrist et al. 2000;Strong et al. 1999). As previously emphasized by McKinley and Siegrist (2011), insights into the three-dimensional (3D) hydrobiogeochemical complexity of soil-based I/A OWTS are critically needed for making more informed decisions on system geometries, matrix composition and compaction, hydraulic loading, and dosing patterns to maximize system performance and resilience and also to minimize installation costs and risks of system failure.
In this study, high-resolution planar optode imaging was used to investigate O 2 dynamics in a bench-scale mimic of a soil-based I/A OWTS during dosed artificial wastewater application and during a hydraulic disturbance. Experimental mesocosms were filled with aged NRB material consisting of an unsaturated loamy sand layer on top of a water-saturated sand-sawdust layer to simulate the sequential nitrification and denitrification function of the in-ground lined NRB. The goal was to coregister hydraulic and biogeochemical processes in functional and hydraulically disturbed mesocosms and to relate O 2 imaging data to the nitrogen removal performance of the system.

Matrix Material
Matrix for a laboratory experiment was obtained from an in-ground (7.6 × 4.5 m wide and 0.92 m deep) nonproprietary NRB that had been operated at the Massachusetts Alternative Septic System Test Center (MASSTC) for approximately 2 years until excavation in late November 2016. The in-ground system was composed of approximately 46 cm of loamy sand in the unsaturated layer and 46 cm of a 1∶1 mixture (v/v) of sand and sawdust in the watersaturated layer. Water-saturated conditions were maintained using a plastic liner, which encased the sand-sawdust layer. The NRB received STE in 8 doses per day at a time-integrated loading rate of 2.4 cm day −1 (225 gallons day −1 over 375 ft 2 surface area) with 40-60 mg N L −1 (Langlois et al. 2020), which is a common application rate in soil-based OWTS and a typical range of total nitrogen (TN) concentration in domestic STE (Beal et al. 2005;Cogger et al. 1988;Lowe et al. 2009). The NRB was excavated because of its declining nitrogen removal performance. Increased NH þ 4 concentrations in the effluent suggested that reduced performance was largely due to decreased nitrification efficiency. The excavation provided a unique opportunity to obtain an aged matrix containing disturbed but potentially active nitrifying and denitrifying microbial communities.

Experimental Optocosms
Six rectangular mesocosms, hereafter referred to as optocosms, were made of 12 mm thick polycarbonate sheets (Lexan) with inner dimensions of 19 × 12.5 × 5.2 cm (Fig. 1). At the bottom of each optocosm, a permeable support structure was included [ Fig. 1(b)] to allow even water percolation. Optocosms were first filled with aged sand-sawdust (denitrifying matrix) and saturated with deionized water (DiW). Effluent ports were connected to Tygon tubing (3.1 mm ID) and T-connectors, which were elevated to mimic water-saturated conditions in the sand-sawdust layer of the lined in-ground NRB. Subsequently, optocosms were filled with moist aged loamy sand (nitrifying matrix) and compacted by placing a standardized mass equivalent to 0.43 N=cm 2 (4.3 kPa) on the matrix surface. The sand and sand-sawdust layer in the optocosms were both 7 cm thick compared to 46 cm in the in-ground system. The sand surface area (65 cm 2 ) was exposed to the atmosphere to allow gas exchange across the matrix surface. An absorbent mat (PIG MAT231) placed on the sand surface ensured even trickling of the water through the entire surface [ Fig. 1(a)]. Water was supplied through a capped PVC pipe (15 mm ID) with two drilled holes (3 mm ID).

Artificial Wastewater Solution
Simplified artificial wastewater solution, hereafter referred to as influent, was prepared by dissolving ammonium nitrate in DiW resulting in final concentrations of 707 μM NH þ 4 (9.9 mg N-NH þ 4 L −1 ) and 742.5 μM NO x (10.4 mg N-NO x L −1 ) for a total nitrogen (TN) concentration of 1,449 μM (20.3 mg N L −1 ). No organic carbon was added to the influent water. The rationale for this simplified influent composition was, on the one hand, to provide the essential substrates for the nitrifying (NH þ 4 ) and denitrifying (NO x ) microbes; while on the other hand, it was to avoid aerobic organic carbon mineralization, which would have complicated data interpretation. Furthermore, the intent was to test whether the organic matter provided by the sawdust was sufficient to promote denitrification in the sand-sawdust layer. The chosen NH þ 4 concentration was 4-to 6-fold lower than what the in-ground system typically received during operation to account for the 6-fold lower thickness of the unsaturated loamy sand layer and to avoid high NH þ 4 supply to the water-saturated sand-sawdust layer.

Experimental Protocol
The experiment was conducted in a temperature-controlled room at 15°C. For the first 2 days, the operation of all optocosms resembled the operation of a lined in-ground system. Specifically, each optocosm received influent in doses of 13 mL delivered over 15 min every 3 h using a computerized multihead peristaltic pump . This corresponded to a loading rate of 1.6 cm day −1 , which is slightly less than what the in-ground system received (see "Matrix Material"). Effluent from each optocosm was collected after a dose application on day 2 to confirm that effluent volume was equal to dose volume (see also Effluent flow rate measurements). After 2 days, a hydraulic disturbance was initiated in three of the six optocosms (Opt 1-3) over the subsequent 8 dosing cycles (∼24 h). This was done by closing the effluent ports while the dosed application continued, thus simulating clogging of the effluent drainage pipe. Optocosms receiving this treatment are hereafter referred to as hydraulically disturbed (HD) optocosms, whereas the optocosms whose effluent ports remained open (Opt 4-6) are referred to as control optocosms (CT). Dosed influent application was continued in all optocosms for the remainder of the experiment. One week after reopening effluent ports, effluent samples were collected during four dosing events over 2 days to quantify nitrogen removal in each optocosm.

Effluent Flow Rate Measurements
On day 2 of the experiment, time-lapse images were taken at 5 s intervals with a DSLR camera (Canon Rebel T5) while collecting the effluent in graded cylinders during and after a dose application for a total of 45 min. These images were used to determine how long after a dosing event the optocosms produced effluent and to determine the cumulative effluent production curves.

Time-Lapse O 2 Imaging
Each optocosm was equipped with a planar O 2 optode ( Fig. 1), which is a fluorescence-based foil sensitive to O 2 partial pressure (pO 2 ), enabling rapid measurements of O 2 concentrations in air and in water (Wang and Wolfbeis 2014). The optodes were prepared, as described in Precht et al. (2004). Optical quantification of pO 2 is based on the detection of luminescence emitted by a porphyrin dye whose lifetime (and intensity) depends on O 2 (Glud et al. 1996). Imaging of pO 2 was done using the luminescence lifetime imaging system originally developed by Holst and Grunwald (2001) and modified by Matsui et al. (2011). The luminescence was excited with blue light-emitting diodes (λmax ¼ 455 nm, LXHL-LR5C, Philips Lumileds) and detected with a CCD camera (Cook, PCO1600). The spatial resolution of the system, determined by the imaged area and the resolution of the camera (1,600×1,200pixels), was 0.37 mm per pixel. Integration time during image acquisition depended on LED intensity and distance but was typically around 1 s to optimize signal-to-noise ratio while ensuring that signals in all pixels, including those in anoxic regions (which have the brightest signal) were below image saturation.
Images of pO 2 at the optocosm wall were recorded in 30 s intervals over the entire experiment (12 days). For each optocosm, the optode was calibrated using a two-point calibration based on the fluorescence lifetime measurements in the atmosphere above the sand matrix and in the anoxic water within the sand-sawdust layer. Then, using the pO 2 derived from the measured fluorescence lifetime, concentrations of O 2 in air and in water were calculated using the ideal gas law (pO 2 ¼ 210 matm corresponds to 8.9 mmol O 2 L −1 of air) and Henry's law (pO 2 ¼ 210 matm corresponds to 0.314 mmol O 2 L −1 of water), respectively, at the experimental temperature of 15°C.
Image processing and analysis were done in MATLAB version R2015a. Based on O 2 images from 21 h of continuous measurement (7 consecutive doses), images of pO 2 and oxic-anoxic oscillation frequencies were estimated as previously described (Volkenborn et al. 2012). Oxic conditions were defined as pO 2 > 10 matm, and the duration of a redox oscillation was defined as a time interval needed to change from anoxic (<10 matm) to oxic (>10 matm) and back to anoxic conditions. The pO 2 threshold of 10 matm was used to ensure that the probability of detecting oxygen under anoxic conditions, such as due to analytical noise in individual pixels, would be negligible (<0.1%). Previous studies showed that the standard deviation of intensity was 5%-8% of the intensity signal (Glud et al. 1996) and that choosing a 5% threshold was a robust approach to discriminate between oxic and anoxic conditions using planar O 2 data (Volkenborn et al. 2012).   The influent was pumped into the distribution pipe, and water percolated through the unsaturated loamy sand layer (7 cm thick) and into the saturated sand-sawdust layer (7 cm thick). Effluent was collected from the tube whose end determined the top level where the sand-sawdust layer was water-saturated. The rectangle represents the oxygen optode which was attached to the inner surface of one side of the optocosm. A support structure seen in panel (b) elevated the sand-sawdust matrix off the bottom of the optocosm creating a water-filled area that allowed even water flow throughout the matrix.
Average pO 2 and number of oxic pixels within regions of interest were calculated for selected time intervals and used to estimate the area of the oxic and anoxic regions at the optocosm wall. The time series of oxic area was then processed in R (version 3.5.1, using the decompose function from the fpp package) to remove the oscillation imposed by dosing (frequency of 1=3 h −1 ). This was done to estimate the average areas of the continuously oxic and anoxic matrix in each optocosm. Subsequently, this area was multiplied by the width of the optocosms in the direction perpendicular to the optode wall (5.2 cm) to estimate the volume of the oxic and anoxic matrix within the optocosms. This calculation assumes that the oxic and anoxic area measured at the wall were representative of such areas across any vertical cross-section of the optocosm. Finally, in each optode pixel, volumetric rates of O 2 consumption were calculated as previously described (Polerecky et al. 2005) using the slope of the linear decline in O 2 concentrations between dosing events.

Matrix Properties
In parallel to the optocosm setup, triplicate cores with either sand or sand-sawdust matrix were prepared and used to measure saturated hydraulic conductivity, porosity (3.8 cm ID and 16 cm matrix depth), and water content (3.8 cm ID and 8 cm sand depth). Cores were prepared by adding a matrix to water-filled cores, with the saturation level maintained at ≤1 cm above the soil-water interface to minimize air entrapment that would have affected the porosity and hydraulic conductivity measurements. Saturated hydraulic conductivity was measured using constant head tests (Klute and Dirksen 1986). Porosity was measured using a gravimetric bulk volume method. Water content was measured gravimetrically using the standard test method for laboratory determination of moisture [ASTM D2216 (ASTM 2010)] after applying 18 mL of water to moist sand. These measurements were done shortly after the water stopped draining from the columns, and thus, the values represent maximal water content after primary drainage.

Analytical Methods
Ammonium and NO x samples were acidified with concentrated sulfuric acid and stored in the refrigerator (≤4°C) until analyzed with a Lachat Quikchem autoanalyzer using methods 10-107-06-2-L (equivalent to EPA 350.1) and 10-107-04-1-A (equivalent to EPA 353.2), respectively. Dinitrogen gas samples were collected in glass exetainers prespiked with 20 μL saturated HgCl 2 , and vials were carefully filled so that no headspace remained. Samples were stored at 15°C and analyzed within 2 weeks of collection using a Membrane Inlet Mass Spectrometer (MIMS) according to Kana et al. (1994). The MIMS signal was calibrated with air-equilibrated DiW with equilibrium concentrations calculated according to Hamme and Emerson (2004). Denitrified nitrate, accounted for as N 2 , was calculated by taking the difference between the effluent and air-equilibrated influent N 2 concentrations.

Statistical Analysis
To correlate O 2 with effluent nitrogen, the authors first processed 21 h of O 2 images taken on Day 8 in Matlab to calculate the number of oxic and anoxic pixels (i.e., oxic and anoxic area) in each optocosm, as previously described. Images from Day 8 were used because they reflect the O 2 conditions in the matrix approximately 1 HRT before effluent nitrogen samples were taken. Pearson correlation was used to test for significant correlations (p < 0.05) between the continuously unsaturated oxic area in the sand layer (ideal conditions for nitrification) and effluent NH þ 4 concentrations and between the continuously anoxic area in the water-saturated sand-sawdust layer (ideal conditions for denitrification) and effluent NO x concentrations. Effluent nitrogen data had normally distributed residuals, according to the Shapiro-Wilks normality test performed in R. Statistica (TIBCO) was used to perform repeated measures analysis of variance (RM-ANOVA) for effluent concentrations of NO x and NH þ 4 and for % N removal using general linear models with type III sum of squares. Corresponding F and p values are reported within the text. Bartlett's test was used to confirm homoscedasticity. Effluent NO x data were square-root transformed to fulfill the homoscedasticity criteria. For effluent N 2 data, a twosample t-test was performed, and the corresponding t and p values are reported within the text.

Hydraulic Properties
The theoretical hydraulic retention time (HRT) of the optocosms was 4.8 days, as calculated by dividing the total water volume in both the loamy sand and sand-sawdust layers by the hydraulic loading rate (104 mL day −1 ), with a 1.2-day contribution of the unsaturated sand layer and a 3.6-day contribution of the watersaturated sand-sawdust layer. The respective HRT contributions are representative of the conditions at the start of the experiment, prior to HD of optocosms 1-3. Each influent dose was applied over 15 min and all replicates started producing effluent within 5 min of influent application (Fig. 2). More than 90% of the dose-volume was recovered as effluent within 25-45 min after dose application.

Oxygen Dynamics during Dosed Artificial Wastewater Application
Oxygen patterns in the optocosms were reproducible over the first 2 days of the experiment and similar in all replicate optocosms (Video S1). The unsaturated sand and water-saturated sand-sawdust layer remained primarily oxic and anoxic, respectively, except for the uppermost layer of the water-saturated sand-sawdust matrix [ Fig. 3(a)]. In this region, the availability of O 2 oscillated once every 3 h [ Fig. 3(b)] with O 2 concentrations reaching up to 105 matm (157 μM O 2 ) shortly after the dosed influent application [ Fig. 4(a)]. Oxygen penetrated over the entire width of the optocosms with each dose, with some indication of preferential flow toward the center (Video S1). The average thickness of the oscillatory zone over the first 2 days was 1.2 cm (AE0.2 cm; n ¼ 6) [ Fig. 3(b)]. The average rate of O 2 consumption in water within this layer was 110 matm pO 2 h −1 (mean of 925 pixels within 9 reactive regions in CT optocosms), consistent with the observation that the regions intermittently supplied with O 2 went completely anoxic within 1 h between subsequent dosing events [ Fig. 4(b)].

Oxygen Dynamics Associated with Hydraulic Disturbance (HD) of the Sand Layer
The location and spatial extent of the oscillatory zone remained relatively constant for Opt 4 and 5 throughout the experiment but changed over time in HD optocosms and Opt 6 due to intentional and unintentional changes in hydraulic function, respectively [ Fig. 3(b)]. After closing the effluent ports of HD optocosms, the pronounced O 2 dynamics in the uppermost 1.2 cm of the sand-sawdust layer associated with dosed wastewater application disappeared immediately [Figs. 4(b and d); Video S1]. The oxicanoxic boundary started to move upwards, and anoxic pockets evolved in the surface layer of the loamy sand matrix [Figs. 4(a and c)]. After 3-5 doses, pO 2 started to decline throughout the unsaturated sand layer of HD optocosms, and in all cases, the decline started immediately after a dosing event [Fig. 4(d)]. Oxygen declined homogeneously and at a constant rate throughout the sand matrix, at rates varying among optocosms between 4.2 and 6.7 matm h −1 [Figs. 4(c and d)]. Reopening of the effluent ports after a 24-h closure resulted initially in a pronounced reoxygenation of the sand layer [ Fig. 4(b)]. However, most of the reoxygenated areas reestablished anoxia within hours and remained anoxic until the end of the experiment [ Fig. 4(b); Video S1]. Four days (i.e., approximately 1 HRT) before nitrogen samples were taken, the oxic unsaturated sand area at the optocosm wall in HD optocosms was 38 AE 5 cm 2 , while it was 62 AE 1 cm 2 in CT optocosms [ Fig. 3(d)].

Discussion
In this study, high temporal resolution O 2 imaging of bench-scale mesocosms resembling soil-based I/A OWTS was used to study hydrobiogeochemical processes within intact saturated and unsaturated matrices during normal operation and following an intentional hydraulic disturbance. The data show that the dosed application of influent led to the establishment of a redox oscillatory boundary layer between the major geochemical zones, i.e., the oxic unsaturated sand layer and the anoxic water-saturated sand-sawdust layer. Additionally, the optode imaging provided new insights into the hydrobiogeochemical processes that are relevant for the function of soil-based I/A OWTS, as discussed in the following section.

Hydrobiogeochemical Function of Soil-Based I/A OWTS Inferred from Oxygen Imaging
In NRBs and other soil-based I/A OWTS, the supply of O 2 to the sand layer is critical, and the lack thereof is one of the main causes of system failure (Beal et al. 2005). Due to the stoichiometry of complete nitrification, simple mass balance calculations show that O 2 supplied with influent water is insufficient to allow for complete nitrification for typical NH þ 4 concentrations in STE, even if influent is 100% equilibrated with air. Furthermore, in-ground systems receive STE with high carbonaceous biological oxygen demand (CBOD) and must also be able to oxidize carbon to function properly (Wilhelm et al. 1994). At 15°C fully air-equilibrated water has an O 2 concentration of 314 μM, while NH þ 4 concentrations in STE are variable but typically much higher (714-7,140 μM N, equivalent to 10-100 mg N-NH þ 4 L −1 ) (Lowe et al. 2009;McCray et al. 2005;Oakley et al. 2010). Based on the ideal gas law, the pO 2 at 1 atm pressure in ambient air at 15°C relates to O 2 concentration in air of 8.9 mmol L −1 . Consequently, despite its rather small volume fraction, the air-filled pore space in NRBs and other soil-based systems, such as leaching fields and SAS, constitutes the major O 2 reservoir.
From a mass balance perspective, the O 2 ∶NH þ 4 ratio in this study was sufficient to allow complete nitrification within the unsaturated sand layer as the amount of O 2 contained in 1 cm 3 of unsaturated sand in the optocosms was 0.88 μmol, and the starting amount of NH þ 4 contained in 1 cm 3 of moist sand matrix in our experiment was about 0. the water content dependency of nitrification rates and predicted maximal nitrification rates at volumetric water contents between 0.15 and 0.26. Higher water contents (>0.26) can promote efficient nitrification if interconnected air-filled pore spaces exist that allow for rapid O 2 diffusion from the adjacent vadose zone (Bradshaw et al. 2013). Measurements of water content variance associated with dosing and drainage were not performed in the present study, but the estimate of water content (0.28) measured shortly after primary drainage is well within the range of good nitrification performance, and optode images indicated that O 2 was efficiently supplied and distributed throughout the loamy sand layer, despite optocosm walls that prevented horizontal O 2 diffusion. Complete nitrification in the matrix with similar volumes of air-filled pore spaces as in the optocosms has been documented for moist, unsaturated sand in field systems and laboratory column experiments (Beach et al. 2005;Walker et al. 1973;Yamaguchi et al. 1996). To maintain optimal nitrifying capability, soil-based I/A OWTS efficient venting is critical to balance O 2 consuming processes such as nitrification and aerobic organic matter mineralization. Oxygen imaging data illustrate that both advective transport (by water moving vertically and drawing air into the system) and diffusive transport within the sand matrix (through interconnected air-filled pore spaces) are important to sustain O 2 supply to the system. Effluent port closure initiated changes in the hydraulic performance of optocosms as inferred from O 2 imaging. Clogging of the effluent pipe in in-ground systems may arise due to microbial growth and fouling, airlocks in the system, or abundant mineral oxide precipitation once anoxic effluent is reoxygenated. While such The areas correspond to averages over the same 21-h interval which panels (a) and (b) were calculated over. The effluent samples were collected on Days 11-12. Thus, the areas represent O 2 conditions within the optocosm matrix approximately 1 HRT before effluent nitrogen samples were taken. Optocosms 1-3 experienced intended hydraulic disturbance (HD), and optocosms 4-6 were controls. Each experimental point of panels (c) and (d) corresponds to an individual optocosm. Opt 1-3 and Opt 4-6 are represented by white and black squares, respectively. scenarios may not represent frequent problems, effluent port closure established conditions that resemble a hydraulically overloaded system, such as during a heavy precipitation event. Effluent port closure also mimicked reduced hydraulic continuity due to gas bubble growth within the matrix as a result of N 2 and CO 2 outgassing from the water-saturated sand-lignocellulose layer below. Effluent port closure was, therefore, a simple hydraulic disturbance that allowed analyzing hydrobiogeochemical processes during system malfunction to gain insights into processes that are critical to maintaining NRB function. On day 7 of the experiment, one of the controls (optocosm 6) developed partial saturation and anoxic conditions in the unsaturated sand layer without effluent port closure, suggesting that similar changes in hydraulic performance could occur without closing effluent ports.
Once effluent ports of HD optocosms were closed, the pronounced advective O 2 supply to the uppermost part of the watersaturated sand-sawdust layer stopped immediately [ Fig. 4(b), Video S1], indicating that the release of oxic water to the sand-sawdust layer depended on effluent leaving the system. Furthermore, upon effluent port closure, the advection of oxygen-rich air into the upper portion of the sand layer along with each dose ceased [ Fig. 4(d); Video S1]. In contrast, O 2 concentrations in the unsaturated sand layer of CT optocosms slightly increased after dosing and remained stable around 200 matm pO 2 throughout the 3-h interval between dosing events [ Fig. 4(d)].
Given the volume of each dose (13 mL) and of the estimated airfilled pore space (0.09 mL cm −3 ) in the optocosms, each dose could potentially have saturated a 2 cm thick sand layer during effluent port closure. Thus, over 8 doses during effluent port closure, the sand layer of HD optocosms could have become completely water-saturated. However, the air-filled pore space not simply saturated bottom-up [ Fig. 4]. Instead, the lack of water release from the sand layer after effluent port closure resulted in the formation of water-saturated pockets in the uppermost centimeter of the sand layer, which was seen as regions with a much faster pO 2 decline than in the regions that remained unsaturated [ Fig. 4(c)]. With continued dosing, a film of water accumulated above the sand while air was trapped within the matrix. Water ponding and gas entrapment have been documented in other studies (Beal et al. 2005;Heilweil et al. 2004;McKinley and Siegrist 2011;Rice 1974). In turn, pO 2 started to decline linearly within the sand matrix of HD optocosms at rates of 4.2 − 6.7 matm h −1 [Figs. 4(c and d)]. The onset of the  linear pO 2 decline throughout the unsaturated sand layer [Fig. 4(d)] most likely reflected the time-point when water ponding effectively prevented, or severely reduced, the diffusive O 2 supply from the atmosphere . However, diffusive exchange within the sand layer through interconnected air-filled pore spaces still occurred and was likely responsible for the relatively homogenously distributed rate of pO 2 decline in HD optocosms during effluent port closure [ Fig. 4(c)]. Others have identified a similar mechanism related to biomat clogging in field soil absorption systems (SAS) and laboratory sand columns (Rodgers et al. 2004;Wilhelm et al. 1994). Observations of a slow O 2 decline immediately following effluent port closure and more rapid O 2 decline once a continuous water film at the surface established suggest that both advective O 2 supply associated with dosing, and diffusive O 2 supply through interconnected air-filled pore spaces are critical to compensate for O 2 consuming processes within the system and to sustain its functionality.
Effluent port reopening caused pronounced oxygenation throughout the sand layer and the upper portion of the watersaturated sand-sawdust layer [ Fig. 4(b); Video S1]. This was most likely due to drainage of oxic water that had accumulated above the matrix and the release of oxic water from unsaturated regions of the sand matrix into the uppermost layer of the sand-sawdust layer. However, this oxygenation was ephemeral, and O 2 was rapidly depleted, specifically within water-saturated regions, suggesting that the hydraulic disturbance of the sand matrix, which was intended to be temporary, resulted in sustained water-saturated regions, indicating that gravity was insufficient to release water from those regions when effluent ports were reopened. Consequently, the oxic areas measured in HD optocosms remained smaller than the oxic areas in CT optocosms for the remaining week (Video S1). Under the assumption that the oxic areas detected at the optocosm walls were representative of the conditions throughout the optocosms (see further discussion under Method Limitations), the oxic matrix volume conducive to nitrification was reduced in HD optocosms for the remainder of the experiment.
The spatial patterns of saturation and associated anoxia observed in HD optocosms, which lasted for days after effluent ports were reopened, is consistent with wetting and drainage hysteresis in a heterogeneous matrix with variant pore structures (Morrow 1970). Water can be trapped in such regions due to either low hydraulic continuity of small pore spaces with the surrounding media so that capillary pressure was insufficient to allow drainage or small pore throats surrounding large water-saturated pores (Ahrenholz et al. 2008). Optocosms were dry-filled with sand with residual moisture content. The packing procedure may have allowed for fine particle redistribution and associated decreased hydraulic conductivity following saturation. The finer texture of the loamy sand Opt 6 did not produce effluent for 2 out of 4 doses due to unintentional clogging. Note that in panels (b)-(d), the left y-axis is in mg N L −1 , and the right y-axis is in μM.
above the coarser sand-sawdust layer may have further contributed to lasting saturation because high residual saturation occurs before water discharges from a comparably fine matrix into the coarser layer below (Gardner 1979). Slow recovery of pore spaces that were initially important for diffusive O 2 supply may have further contributed to the observed O 2 decline.
It seems plausible that similar processes led to the decline in the nitrification performance of the in-ground NRB system from which the matrices used in this study originated. The characteristics of the dry-packed loamy sand with a fine particle content of 3.5% by weight (M. Shahraki, personal communication, 2019), combined with regions characterized by strong capillary forces, seemingly constrained the hydraulic resilience of the system after a single saturation event. Although the use of loamy sand is in accordance with the present regulatory guidelines (Otis et al. 1980), the present study suggests that finer-textured sand layer with a heterogeneous pore structure above coarser sand-sawdust layer, such as inlined NRBs, may increase the risk that water-saturated and anoxic conditions establish within the bottom portion of an otherwise unsaturated and oxic sand layer leading to a decrease in its nitrifying capacity. Such anoxic regions may be more prone to clogging due to higher biological biomass compared to aerobic regions (McKinley and Siegrist 2011), hence their development can further impede proper hydraulic function and may result in a positive feedback mechanism promoting the establishment of more anoxic regions.

Nitrification Rates Inferred from Oxygen Imaging
Unlike real septic tank effluent (STE), the water supplied to the optocosms in this study did not contain any organic matter. In addition, NH þ 4 concentrations were 4-to 6-fold lower than what the field NRB received to account for the 6-fold lower sand layer thickness. Thus, the experiment was designed to support complete nitrification. Additionally, the concentration of labile particulate organic matter in the sand matrix was likely very low because it had been exposed to conditions that allowed for aerobic organic matter mineralization for several days prior to the experiment setup.
The linear pO 2 decline during the failed hydraulic performance [Figs. 4(c and d)] suggests inhibited O 2 supply to the unsaturated sand layer. Under the assumption that no process other than nitrification contributed to O 2 consumption, the liner O 2 decline be used to estimate nitrification rates.
In the unsaturated sand matrix, pO 2 declined at rates between 4.2-6.7 matm h −1 [Figs. 4(c and d)]. Considering that the volume fraction occupied by air after primary drainage was 0.09, this corresponds to volumetric oxygen consumption rates (OCR) of 17.9-28.5 nmol h −1 per cm 3 of matrix. Based on the stoichiometry of complete nitrification , this OCR relates to volumetric nitrification rates between 8.9 and 14.3 nmol N cm −3 h −1 , which are within the range of measured and modeled nitrification rates in OWTS Morales et al. 2016). Homogenous pO 2 decline throughout the unsaturated sand layer [Fig. 4(c)] suggests that nitrification occurred homogenously throughout the entire matrix following zero-order kinetics. However, NH þ 4 and NO − 2 oxidation hot spots within the matrix, which may have been present in the optocosms due to preferential water flow and microbial growth, could have been masked by fast diffusion of O 2 within the sand layer through interconnected airfilled pore spaces.
When considering the higher estimates of volumetric nitrification rates (14.3 nmol cm −3 h −1 ), complete nitrification is plausible even during hydraulic disturbance in Opt 1-3, considering the hourly NH þ 4 application rate (3.1 μmol h −1 ) and the oxic unsaturated volumes of sand matrix that remained as a result of intended hydraulic failure (38 cm 2 × 5.2 cm ¼ 198 cm 3 ). Because nitrification only happened in the water-filled pore space, the estimated rate of NH þ 4 concentration decline in the matrix water was 31.8-51.1 μmol N L −1 h −1 . These rates imply that during normal operation, water would have to stay in the unsaturated sand matrix within the optocosms for 14-22 h before all NH þ 4 would be oxidized. This is slightly shorter than the HRT of 1.2 days estimated for the sand layer in the optocosms (see section "Hydraulic Properties") and, therefore, consistent with almost complete nitrification in this experiment.
Upscaling the estimated nitrification rates of 8.9-14.3 nmol cm −3 h −1 to the in-ground NRB implies that 1 m 2 of a 46 cm thick unsaturated sand matrix can nitrify 98 − 158 mM NH þ 4 m −2 day −1 (1.4 − 2.2 g N m −2 day −1 ). At a loading rate of 20.6 L −1 m −2 day −1 that was used in the in-ground system, this means that the system would be capable of completely nitrifying STE with an NH þ 4 concentrations between 4.9 and 7.6 mmol L −1 (68-107 mg N L −1 ), which is well within the typical range of STE (Lowe et al. 2009). Based on the rate of NH þ 4 decline in the matrix water in the optocosms (0.76 − 1.23 mmol N L −1 day −1 ), water would, therefore, need to stay for 6.4 days before NH þ 4 in the STE with a concentration of 4.9, and 7.6 mmol N L −1 would be completely nitrified. This is consistent with the HRT of the in-ground system: assuming that the water content in the unsaturated matrix was similar to the estimates presented in this study (0.28), it would take 6.2 days to replace 127.8 L of water within 1 m 2 of a 46 cm thick matrix at a loading rate of 20.6 L m −2 day −1 .
Apart from relatively slow pO 2 decline throughout the unsaturated sand layer (4.2 − 6.7 matm h −1 ), pO 2 declined much faster (111 matm pO 2 h −1 ) in confined regions within the top of the layer [ Figs. 4(a and c)]. These were regions that most likely became locally water-saturated and turned anoxic due to the limited O 2 supply. Considering that the pore space in these regions (porosity 0.37) was then filled completely with water, the measured pO 2 decline (111 matm h −1 ) corresponds to a volumetric OCR of 38.9 nmol O 2 h −1 cm −3 of matrix. This is slightly higher but, given the uncertainty in the local porosity, well within the range of volumetric OCR estimated for the unsaturated regions (17.9-28.5 nmol h −1 cm −3 ). Thus, due to the much smaller reservoir of O 2 in water than in air, similar OCRs, most likely driven by nitrification, resulted in much faster O 2 decline and in the rapid formation of anoxic regions when the matrix became water-saturated.

Denitrification Rates Inferred from Oxygen Imaging
The optocosms in this study removed 87%-97% of N from the artificial influent with 1,449 μM N (20.3 mg N L −1 ) [ Fig. 5(a)]. In HD optocosms, 100% of the sand-sawdust layer was continuously anoxic [ Fig. 3(a)], and TN removal was 97% [ Fig. 5(a)]. The reduction of continuously anoxic sand-sawdust matrix volume, due to the establishment of a 1.2 AE 0.2 × 5.2 × 12.5 cm oscillatory zone, is consistent with the decreased nitrate removal in control optocosms [Figs. 3(b) and 5(a)]. Effluent NO x concentrations strongly correlated with the areas of continuously anoxic sand-sawdust matrix [ Fig. 3(c)], suggesting that the twodimensional (2D) oxygen data were representative of conditions within the 3D matrix. It also suggests that, due to low organic carbon availability, there was likely minimal denitrification occurring in the sand layer regions that became anoxic. Others have similarly reported that less than 5% of the applied nitrogen was removed by denitrification in fine sands (Beggs et al. 2011;Tindall et al. 1995).
Considering the matrix volumes that were conducive to denitrification (continuously anoxic and presence of organic carbon) in control optocosms, effluent NO x concentrations correspond to a volumetric denitrification rates of 13.8-14.3 nmol N h −1 cm −3 matrix (equivalent to 4.6-4.8 mg N day −1 L −1 of matrix or 8.6-8.9 mg N day −1 L −1 of matrix water). At this denitrification rate, the time needed to completely denitrify influent nitrogen is approximately 58 h (or 2.4 days). This is consistent with almost complete nitrogen removal in optocosms with a 3.6-day HRT in the sandsawdust layer. Denitrification rates presented in this study are consistent with published rates (Jaynes et al. 2016) and similar to rates determined by sealed incubations of sand-sawdust matrix obtained from the same in-ground system during excavation (Waugh et al. 2020). Upscaling this rate to the in-ground NRB with 46 cm sandsawdust layer, a denitrification rate of 13.8-14.3 nmol N cm −3 h −1 is equivalent to complete nitrogen removal in an NRB operating at a loading rate of 20.6 L m −2 day −1 with wastewater nitrogen concentration of 7.4-7.6 mmol N L −1 as NO x (103-107 mg N L −1 ).

Oscillatory Redox Conditions in NRBs
High-resolution oxygen imaging allowed the identification of a 1.2 cm thick horizon characterized by oscillatory redox conditions in the uppermost layer of the water-saturated sand-sawdust matrix of the control optocosms [ Fig. 3(b)]. Oxygen dynamics in this zone were clearly linked to the dosed wastewater application and most likely reflected the release of oxic water from the unsaturated sand layer to the water-saturated sand-sawdust layer after a dose was applied [ Fig. 4(b); Video S1]. In HD optocosms, the oscillatory zone existed further up in the loamy sand layer after the hydraulic disturbance, and the spatial extent was smaller than in CT optocosms, suggesting that porosity and texture differences between the sand and the sand-sawdust layer may have amplified the spatial extent of the oscillatory redox zone. It is expected that a redox oscillatory zone also exists in in-ground systems. However, due to the longer flow path and increased CBOD in STE, the magnitude of the redox fluctuations, as well as the location and thickness of the oscillatory zone, may differ from those observed in our experimental optocosms.
Nevertheless, the oscillatory redox zone deserves further attention because it constitutes a unique microenvironment whose function and environmental impact are presently not understood. Nitrification and denitrification are both multistep pathways regulated by O 2 availability (Kuypers et al. 2018;Zhu et al. 2013). Oscillations between oxic and anoxic conditions can lead to the formation of the metabolic intermediates, including nitrite (NO − 2 ), nitric oxide (NO), and nitrous oxide (N 2 O) (Caranto and Lancaster 2017;Gabarro et al. 2014;Kampschreur et al. 2009). For example, O 2 concentrations between 1% and 30% air saturation observed at the redox interface of the NRB optocosms [ Fig. 4(a)] represent conditions that can result in the production of N 2 O by nitrifiers (Codispoti 2010). Even if the volume of the redox oscillatory zone in full-scale NRB were small compared to the total volume of the system, this dynamic microenvironment might be characterized by high rates of biogeochemical cycling. Furthermore, the gaseous intermediates produced in this zone, which are environmentally relevant because of their greenhouse (N 2 O) and ozone-destruction (NO) potential, may rapidly escape into the atmosphere through the interconnected air-filled pore spaces in the unsaturated sand layer. Further research is required to quantify the fluxes of these gases into the atmosphere and thus better understand the environmental impact of the oscillatory zone in NRBs.

Method Limitations
Planar optode imaging measures the concentration of a chemical species at the mesocosm wall. In some of the calculations presented in this study, it was assumed that the measured values of pO 2 were representative for the O 2 distribution throughout the 5.2 cm wide optocosms. The fact that the oxic areas in the unsaturated sand layer and the continuously anoxic areas in the saturated sand-sawdust layer scaled linearly with effluent concentrations of the reduced and oxidized nitrogen species, respectively, is consistent with the assumption that the 2D data were representative for the 3D hydrobiogeochemical conditions. Microsensor measurements in the bulk volume of the matrix should be included in future studies to verify this assumption.
As in all laboratory incubations, the sidewalls of the optocosms constituted impermeable barriers that prevented horizontal O 2 flux through the unsaturated loamy sand. In in-ground systems, horizontal O 2 fluxes from the surrounding vadose zone into soilbased OWTS may allow for continued O 2 supply even if the water is ponding at the surface and residual water content is high (Bradshaw et al. 2013). Thus, if top-down O 2 fluxes are reduced or prevented by surface ponding, O 2 limitation in the unsaturated matrix may be less of a problem in in-ground systems. However, it is worth mentioning that in some systems, the plastic liner may extend above the water-saturated layer so that lateral oxygen diffusion along the outer edges of the NRB unsaturated sand layer may be inhibited as well.
It is important to note that the rate estimates presented in this study are very sensitive to the fraction of air-filled pore space. A two-fold larger volume fraction of air-filled pore space would relate to a two-fold higher nitrification rate estimate. Unsaturated sand water content most likely varied during wetting and draining phases associated with dosed water application. The air-filled volume fraction of 0.09 used in calculations of oxygen consumption and nitrification rates is likely a conservative estimate because it was based on water contents measured shortly after primary drainage. Similar air contents have been reported for hydraulically overloaded or overcompacted systems (Beggs et al. 2011;De and Toor 2015;Motz et al. 2012).
The relative and absolute volumes occupied by particles, water, and air, as well as the connectivity of air-filled pore spaces, also depend on the magnitude of compaction during NRB construction and natural consolidation over time. In fact, compaction of soils with high concentrations of fines has been shown to reduce the size of air-filled pore space, resulting in decreased venting of the unsaturated sand layer (Walker et al. 1973). Performance of overcompacted soils and soil-based OWTS with clogged biofilms may be restored by drying and aeration to regenerate air-filled pore space, but sometimes complete replacement of the sand is required (Potts et al. 2004;Walker et al. 1973). Precise data on matrix porosity, water content hysteresis, and air-filled pore space, as well as on the spatial heterogeneity and temporal variations of these parameters, will be critical to further constrain nitrification rates based on O 2 imaging. Influent in this study was different from typical STE in that it contained similar concentrations of NH þ 4 and NO x and no CBOD. Although high CBOD could potentially have affected the vertical O 2 distribution within the unsaturated layer, it seems plausible that oxic conditions would have been maintained even if additional O 2 consumption by organic matter mineralization occurred, as long as interconnected pore spaces were present to allow O 2 to flux in from the atmosphere. Consequently, water released from the unsaturated matrix to the water-saturated matrix after dose application would have been oxic and would have resulted in similar oxic-anoxic oscillations in the uppermost water-saturated layer. However, in HD optocosms CBOD may have resulted in a more rapid and extensive establishment of anoxic regions. Under these circumstances, and in the presence of organic-rich STE, denitrification could have been promoted in anoxic regions.
The process rates presented in this laboratory study were determined at a constant temperature and thus cannot simply be applied to in-ground systems in which temperatures and nitrification rates vary seasonally (Bradshaw et al. 2013;Hwang and Oleszkiewicz 2007;Zhu and Chen 2002). However, the agreement between the nitrification rates derived from O 2 measurements, nitrogen transformations measured in this experiment, and previously published rates suggest that O 2 imaging is a promising tool to further investigate the temperature and water saturation dependence of microbially mediated pathways involving nitrogen transformations in soil-based I/A OWTS.

Conclusions
For the first time, planar optode imaging was used to study O 2 dynamics within laboratory mimics of a soil-based I/A OWTS. Real-time visualization of O 2 dynamics confirmed the expected geochemical zonation in unsaturated and water-saturated two-layer systems such as NRBs and generated new insights into hydrobiogeochemical processes in NRBs during dosed wastewater application and under hydraulically disturbed conditions. First, O 2 imaging enabled identification and real-time monitoring of regions with optimal conditions for nitrification and denitrification. The differences in the spatial extent of these regions scaled linearly with effluent NH þ 4 and NO x concentrations, suggesting that 2D O 2 measurements were representative of the 3D conditions within the optocosms and that zero-order kinetics were appropriate to describe the nitrogen transformation pathways. Nitrogen removal rates in the optocosms were consistent with efficient nitrogen removal of in-ground NRBs. Second, O 2 imaging revealed the presence of a redox oscillatory zone at the interface between the unsaturated sand and water-saturated sand-sawdust layers. Further research is required to better understand the microbial diversity, biogeochemical function, and environmental impact of this distinct microenvironment. Third, the imaging approach helped identify the likely mechanism responsible for the declined nitrification performance of the in-ground NRB from which the material for this experiment derived. A single saturation event resulted in the establishment of lasting water-saturated regions, in which capillary forces were too high to release water by gravity. Partial water saturation limited the diffusive O 2 supply to the sand layer (due to the loss of interconnected air-filled pore spaces) and, ultimately, led to a decrease in its nitrification efficiency. Although the use of loamy sand in soil-based I/A OWTS is in accordance with the present regulatory guidelines, this study suggests that loamy sand may not be an ideal nitrifying matrix when situated above a coarser textured denitrifying matrix in a lined NRB. Overall, this study highlights the critical role of unsaturated matrix composition and compaction for the system performance and resilience and illustrates the potential of high-resolution chemical imaging to coregister hydrological and biogeochemical processes in soil-based I/A OWTS. Future studies that make use of planar optodes sensitive to other chemical species such as pH or pCO 2 have the potential to further improve current understanding of I/A OWTS performance as a function of matrix composition and compaction, hydraulic loading, and dosing patterns. Such data are critically needed to improve the parameterization of reaction-transport models and to assist in making more informed decisions on soil-based I/A OWTS design and operation.