Modeling Onsite Wastewater Treatment System Contaminants in Current and Climate Changing Conditions

The use of onsite wastewater treatment system (OWTS) is a common practice in the U.S., especially in rural areas where the access to centralized wastewater treatment systems is limited. Onsite wastewater treatment systems include a soil treatment area or drainfield where contaminants are removed or attenuated. Ineffective OWTS are a source of microbial pathogens (bacteria and viruses), biological oxygen demand (BOD) and nutrients, which are among the major causes of contamination and water quality impairments in surface water in the U.S. The main objective of this research was to model the different chemical, physical processes, and removal mechanisms that influence the fate and transport of OWTS-derived contaminants using the HYDRUS 2D/3D software. In the first part of this study, segmented mesocosms (n=3) packed with sand, sandy loam or clay loam soil were used to determine the effect of soil texture and depth on transport of two septic tank effluent (STE)-borne microbial pathogen surrogates – green fluorescent protein-labeled E. coli (GFPE) and MS-2 coliphage – in soil treatment areas. In all soils, removal rates were >99.99% at 25 cm. The transport simulation compared (1) optimization, and (2) trial-and-error modeling approaches. Only slight differences between the transport parameters were observed between these approaches. Independent of the fitting procedure, attachment rates computed by the model were higher in sandy and sandy loam soils than clay loam, which was attributed to unsaturated flow conditions at lower water content in the coarser-textured soils. In the second part of this research, bacteria removal efficiencies in a conventional soil-based wastewater treatment system (OWTS) were modeled to elucidate the fate and transport of E. coli under environmental and operational conditions that might be expected under changing climatic conditions. The impact of changing precipitation patterns, initial bacteria concentrations, hydraulic loading rates (HLR), and higher subsurface temperatures at different depths and soil textures on bacteria removal was evaluated. Modeled effects of initial bacteria concentration shows that greater depth of treatment was required in coarser soils than in fine textured ones to remove E. coli. The initial removal percentage was higher when HLR was lower, but it was greater when HLR was higher. When a biomat layer was included in the transport model, the performance of the system improved by up to 12.0%. Lower bacteria removal (up to 5%) was observed at all depths under the influence of precipitation rates ranging from 5 cm to 35 cm, and 35 cm rainfall combined with a 70% increase in HLR. C Increased subsurface temperature due to climate change (23 C) increased bacteria removal relative to a lower temperature range (5 C to 20C). It appears that the performance of OWTS may be impacted by changing climate. In the third part of this research, we also simulated the fate and transport of N in three different types of OWTS drainfield, or soil treatment areas (STA) using 2D/3D HYDRUS software to develop a N transport and fate model. Experimental data from a laboratory mesocosm study, including soil moisture content, and NH4 and NO3 concentration, was used to calibrate the model and a water content-dependent function was used to compute nitrification and denitrification rates. Three types of drainfields were simulated: (1) pipe-and-stone (P&S), (2) pressurized shallow narrow drainfield (SND) and (3) Geomat (GEO), a variation of SND. The results showed that the model was calibrated with acceptable goodness-of-fit between the observed and measured (average root mean square errors (RMSE) ranged from 0.18 to 9.65 for NH4 and NO3). The model predicted the N losses from nitrification and denitrification in all STAs. The modeled N losses occurred mostly as NO3 in water Outputs, accounting for more than 82% of N inputs in all drainfields. The highest N losses by denitrification were computed for the P&S drainfield and accounted for 17.60% of the influent total N. Our results showed that HYDRUS is a useful tool to predict the fate and transport of nutrients and microbial contaminants and help to provide practitioners with guidelines to estimate pathogens and nutrients removal efficiencies for OWTS under the effect of different operational and environmental factors. In addition, the modeling approach presented in this study, will be useful to predict the extent of contamination and spatial distribution for identifying non-point sources, and establish total minimum daily loads (TMDLs).

climatic conditions. The impact of changing precipitation patterns, initial bacteria concentrations, hydraulic loading rates (HLR), and higher subsurface temperatures at different depths and soil textures on bacteria removal was evaluated. Modeled effects of initial bacteria concentration shows that greater depth of treatment was required in coarser soils than in fine textured ones to remove E. coli. The initial removal percentage was higher when HLR was lower, but it was greater when HLR was higher. When a biomat layer was included in the transport model, the performance of the system improved by up to 12.0%. Lower bacteria removal (up to 5%) was observed at all depths under the influence of precipitation rates ranging from 5 cm to 35 cm, and 35 cm rainfall combined with a 70% increase in HLR. C Increased subsurface temperature due to climate change (23 o C) increased bacteria removal relative to a lower temperature range (5 o C to 20 o C). It appears that the performance of OWTS may be impacted by changing climate. In the third part of this research, we also simulated the fate and transport of N in three different types of OWTS drainfield, or soil treatment areas (STA) using 2D/3D HYDRUS software to develop a N transport and fate model. Experimental data from a laboratory mesocosm study, including soil moisture content, and NH 4 and NO 3 concentration, was used to calibrate the model and a water content-dependent function was used to compute nitrification and denitrification rates. Three types of drainfields were simulated: (1) pipe-and-stone (P&S),   . Those Table 3. 4 Modeled water-filled pore space or relative water content for all STA types.131 Table 3. Onsite wastewater treatment systems treat residential wastewater -consisting of both black water (urine and fecal matter) and grey water (shower, laundry, kitchen) (EPA, 2002) -in a series of steps that begin with primary treatment, or initial settling of bulk solids, in the septic tank. In a conventional pipe and stone (P&S) system, septic tank effluent (STE) is dispersed directly to a soil treatment area (STA), or drainfield, without further treatment. In advanced technologies, prior to dispersal onto the STA, STE is treated further to achieve substantial removal of contaminants, reducing the contaminant load to the STA and lessening reliance on the soil for wastewater renovation.
Ineffective OWTS are a source of microbial pathogens (bacteria and viruses), which are one of the major cause of contamination and water quality impairments in surface water in U.S. (US EPA, 2010). Many pathogenic microorganisms require relatively small numbers to cause infection and induce illness in humans. In order to avoid microbial contamination, US EPA recommends a separation distance of 45 cm  between the infiltrative of surface of the STA and the water table, regardless of soil   chemical and physical properties (US EPA, 2002).
Another wastewater-derived contaminant that can produce water impairments is nitrogen (N), particularly nitrate (NO 3 -). The presence of nitrate in drinking water is the principal cause of methemoglobinemia, which affects the ability of red blood cells to bind oxygen (Shuval & Gruener, 2013). Non-Hodgkin's lymphoma is assumed to be caused by NO 3 when present in concentrations above 4 mg N/L in drinking water (Ward et al., 1996). In addition, high nitrate loads discharged into surface water or marine environments can affect water quality by stimulating eutrophication (Brandes et al., 1974;Weiskel and Howes, 1992  . These models consider that the movement of solute through the soil is the result of the physical process of convection, or mass flow, of water and the chemical process of diffusion in response to a concentration gradient (Addiscott and Wagenet, 1985). The solute transport is characterized in a porous media by an ensemble of pore velocities that exist due to microscopic and macroscopic variations in pore size, tortuosity of flow path and the distribution of both water and solutes within partially water-filled pores. A solute introduced into such a flow system will thereby spread, or disperse, as it is convectively and diffusively transported through the soil.  VIRALT, is a modular semi-analytical and numerical code that simulates the single-source transport and fate of viruses in the saturated and unsaturated zones (Park et al., 1992). The code considered that viruses are transported by advection, dispersion and sorption. Also includes virus inactivation.
 CANVAS, is a model derived from VIRALT and in addition, the code models transport by colloidal matter and the simulation of multiple contaminants sources (Park et al., 1993).
 VIRTUS, which stands for "virus transport in unsaturated soils", the model predicts the virus fate and transport in unsaturated soils and allows the virus inactivation as a function of changes in soil temperature (Yates and Ouyang, 1992 were removed completely by mechanical filtration ) and adsorption to soil particles , respectively. The model determined water-content dependent attachment-detachment rates for microbial pathogens in order to calculate reduction values, which was higher in fine-texture soils than in granular soils at a given water content.
Projections of climate conditions in parts of US, including the Northeast, indicate that sea level, rainfall rates and temperatures have been on the rise and will continue to do so during the next 100 years (Kirtman et al., 2013). The effect of climate of change and sea level rise may affect the performance of the OWTS in coastal areas or in areas with shallow water tables. Changes of ambient temperature will influence the availability and consumption of oxygen by soil microorganisms, which have an effect on microbial processes. Sea level rise, as well as increased precipitation and infiltration may reduce the vertical separation between the infiltrative surface and the groundwater. As a result, less unsaturated soil will be available and the ability of soil to remove contaminants (BOD, N, P and pathogens) through chemical and biological processes may be diminished.
In  (2) trial-and-error modeling approaches. Only slight differences between the transport parameters were observed between these approaches.
Treating both the die-off rates and attachment/detachment rates as variables resulted in an overall better model fit, particularly for the tailing phase of the experiments. Independent of the fitting procedure, attachment rates computed by the model were higher in sandy and sandy loam soils than clay, which was attributed to unsaturated flow conditions at lower water content in the coarser-textured soils. Early breakthrough of the bacteria and virus indicated the presence of preferential flow in the system in the structured clay soil (clay soil, GA), resulting in faster movement of water and microbes throughout the soil relative to a conservative tracer (bromide). are known affect STU performance, which may lead to differences in removal of viruses and bacteria [1][2][3].

INTRODUCTION
A number of studies have investigated the removal efficiency of bacteria in STAs and the processes involved. Crites [4] suggested that bacterial removal or inactivation in STAs is associated with predation by bactrivorous organisms and exposure to sunlight.
Mechanical filtration and adsorption, and flow rate also have a significant effect on removal of pathogenic bacteria [2,5,6]. All of these processes are influenced by soil texture and structure. Fine textured and poorly structured soils are expected to remove bacteria mainly through mechanical filtration because of the smaller pore sizes and lower hydraulic conductivity of those soils. Together with a greater surface area, this results in higher rates of bacteria adsorption [2]. In contrast, coarse and well-structured soils have larger pores and lower porosity values, which allows for better aeration that promotes microbial predation and attenuation [7].
Viruses are thought to be removed in STAs through adsorption to soil particles rather than by mechanical filtration [8,9]. Viruses have a smaller diameter compared to soil pores, which prevents them from being trapped in the pore space. Adsorption of viruses is a function of the physical and chemical properties of the soil, particularly pH, organic matter content and water content [10][11][12][13].
Mature OWTS systems develop a biological growth layer of low permeability at the infiltrative surface of the STA, known as a biomat. Typically, the biomat extends up to 2 cm below the water-soil interface [14,15]. It may enhance the inactivation of microbes through mechanical filtration because partial clogging of smaller soil pores results in reduced infiltration rates and the development of unsaturated flow conditions in the underlying soil profile [14][15]. Unsaturated flow conditions result in longer contact times between microbes and soil particles, which improves the pathogen removal efficiency of the soil treatment zone [15][16].
The retention of microorganisms in soil can be affected by preferential flow, which may be associated with pathways created by plant roots and earthworms, the presence of interaggregate spaces [17,18], and differences in hydraulic conductivity within the soil strata [19]. Preferential flow increases the travel velocity of the aqueous phase, allowing for faster and deeper movement of microbes into the soil profile [20][21][22].
The complex nature of pathogen removal and inactivation in the STA presents a difficult problem with respect to predicting OWTS effectiveness. Contaminant transport models can be used to predict the microbial transport in soils and to help elucidate the factors that control microbial fate as STE moves through the soil profile.
Several models have been developed to simulate virus and bacteria transport in soil. The commercially-available HYDRUS software package is widely used to simulate microbial transport and fate processes, including the transport of viruses, bacteria, and colloids based on either attachment/detachment theory or filtration theory in variably saturated porous media [22][23][24][25][26]. The model supports an interactive graphics-based user interface, and the computational program numerically solves the Richards equation for variably saturated water flow, and the advection-dispersion equations for both heat and solute transport. There are HYDRUS versions available with one-, two-and three-dimensional transport modeling capabilities.
The use and calibration of sophisticated transport models, like HYDRUS 2D/3D, permits investigation of the role of microbial inactivation, removal, and transport processes in homogeneous/heterogeneous soil media by quantifying parameters, such as die-off rates in water and soil or attachment/detachment rates [23]. The calibrated transport parameters can be used to calculate microbial removal as a function of distance between the infiltrative surface and the water table, thus permitting comparison among different soils.
HYDRUS is a valuable and accepted tool for drinking water protection and water resources management purposes. Because of its many capabilities and multi-dimensional functionality, HYDRUS 2D/3D was chosen for modeling our test data.
The objectives of our research project were to: (1) determine the extent to which removal of two microbial pathogen surrogates -a coliphage virus and a tracer bacterium -is affected by soil texture and depth, (2) measure the survival of the coliphage virus and tracer bacterium in sterile and non-sterile unsaturated soil and in STE, and (3) model microbial transport and estimate transport parameters. The results were intended to define and evaluate the potential risk of microbial contamination of groundwater resulting from soil-based treatment of STE. In this paper, we focus on the modeling of microbial transport and how different approaches to modeling -numerical optimization versus visual assessment -best describe experimental data. .  STE was analyzed for dissolved oxygen immediately after collection using the azide modification of the Winkler titration method [27]. The pH was determined using a combination pH electrode and a Model UB-10 pH meter (Denver Instruments, Denver, CO). STE was analyzed for fecal coliform bacteria using the membrane filtration method [27], and, for bacteriophage capable of growing on E. coli (K12), using the plaqueforming assay of Adams [28] [1]. Five-day biochemical oxygen demand (BOD 5 ) was determined following standard procedures [27]. Total P and total N were measured in STE using the persulfate digestion method [27], followed by colorimetric analysis [29][30].

Bromide tracer
Bromide (Br-) is a conservative tracer that permits measuring the breakthrough time of the aqueous solution and relates it to the (retarded) transport of either the bacterial or viral tracers. Tracer tests were conducted by spiking the STE influent with KBr (~20 mg Br -/l). Bromide concentrations were measured using the method of Lepore and Barak [2].
The bromide tracer test data were analyzed with the public domain model CXTFIT to determine the dispersivity (λ) value of each test material [32]. The data were then used for calibration of the transport model.  [27].

Virus tracer
The bacteriophage MS-2 was used as a tracer. MS-2 is a single-stranded RNA coliphage with a 25-nm diameter and an isoelectric point of 3.9 [13]. E. coli strain K12 was used as the host for the bacteriophage. MS-2 bacteriophages were obtained from the Colorado School of Mines (Golden, CO). For each virus addition experiment, MS-2 was diluted in PBS to ~5 × 10 6 pfu/ml and added as described above for the E. coli tracer experiment.
The bacteriophage in the collected samples were enumerated using the plaque-forming assay of Adams [28] on LB agar plates, which were incubated for ~4 h at 37°C, followed by incubation at room temperature overnight before counting plaques in the host lawn.

Survival in soil and STE
Experiments were conducted to determine the survival of the microbial pathogen surrogates in soil and STE. For soil, 2 g (air-dry weight) of soil from each of the three soil types were placed in plastic scintillation vials, in triplicate. Prior to use, the soil was either air-dried or sterilized (121°C for 60 min on 5 consecutive days

Soil Properties
The three soils (sand, sandy loam, and clay loam soil, respectively) were analyzed prior to the start of the experiment and after STE dosing for 27, 31 and 44, weeks. After 27 weeks, all mesocosms had developed a biomat layer that extended over the entire thickness of the gravel layer (4 cm) at the infiltration surface. The total carbon and nitrogen content of the soil was determined using a Carlo Erba EA1108 CHN analyzer (Lakewood, NJ). The soil pH was determined using a 1:5 soil/water ratio with a combination pH electrode and a Model UB-10 pH meter (Denver Instruments). Particle size analysis was conducted using the pipette method [35]. The water content was determined gravimetrically.

Bacteria and Virus Transport Modeling
HYDRUS 2D/3D was used to simulate the transport of microbes in the segmented mesocosms at different depths. The model simulates virus and bacteria transport and fate processes based on a modified form of the convection-dispersion equation [23] (Eq.1): were k a is the first-order attachment coefficient [T -1 ] and k d the first-order detachment coefficient [T -1 ]. According to Simunek [23], the attachment and detachment coefficients are strongly dependent upon the water content, with attachment significantly increasing as the water content decreases. Linear adsorption kinetics were assumed. The chemical non-equilibrium model was used with 50% of all sorption sites assumed to sorb instantaneously and the other 50% are governed by kinetic sorption. properties of the porous materials were obtained from the HYDRUS soil catalog [36].
Based on literature data, the diameter was set at 1.1 μm for E. coli and 0.025 μm for the MS-2 coliphage [37].

RESULTS
The bromide tracer test data and the code CXTFIT 2.1 was used to determine the column system dispersivity (λ) for all three soils. Model fits were good with R 2 values ranging from 0.97 to 0.99. The dispersivity value calculated by CXTFIT 2.1 was approximately 0.289 cm, which is typical for these types of experiments, and is consistent with the range of values (0.06 to 0.816 cm) reported by others [38][39][40][41][42][43]. Next, the hydraulics of the HYDRUS model domain was calibrated using the conservative tracer breakthrough curves (BTC). The tracer test results were fitted for each of the three soil column depth intervals (0-4 cm, 4-10.5 cm, and 10.5-25 cm). The data obtained at the 31 cm sample port was not fitted because E.coli and MS-2 phage concentrations were always below detection limit at that depth. The model results were plotted against the observed data (Figure 1. 3).
The experimental bacterial transport data were fitted to HYDRUS utilizing the model's attachment/detachment module. The data were fitted in two steps: (1) inverse solution, keeping constant the STE and soil die-off rates values (Table 1. 3), to determine the optimized attachment/detachment rates, and (2) a trial-and-error process in which dieoff and attachment/ detachments rates were modified simultaneously until an acceptable graphical fit was achieved. During the trial-and-error process, the emphasis was on achieving the best fit of the tailing end of experimental data. The best-fit simulations of the bacteria and virus test data from mesocosms are shown in Figure 1. 4 and 1.5 .
The experimental data and the model results were plotted both as normal-normal and lognormal graphs to emphasize the two principal phases of these experiments, i.e. the early, high concentration breakthrough and subsequent tailing phase characterized by low microbe concentrations. GFP E.coli concentrations were generally underestimated by the optimization simulation, while a fairly good fit was achieved by the trial-and-error procedure, particularly for the tailing phase. The normal-normal and log-normal plots of the modeled E.coli concentrations captured the oscillations caused by periodic dosing of the column system with STE. The measured bacteria data do not show these "oscillations" because the effluent sampling frequency was not sufficiently high to capture these changes. Initial and peak concentrations simulated at each sampling port and soil type tested are shown in Table 1 Table 1. 5. The goodness-of-fit (R 2 ) of the model was 0.83 or greater for the bacteria simulations, and R 2 =0.76 for the virus data. The nature of the graphical best-fit procedure precluded calculation of R 2 values for the trial-and-error simulations. The liquid (Sink L ) and solid (Sink S ) phase GFP E. coli die-off rates in the optimization and trial-and-error simulations were generally within a factor of three of each other, except for Sink S for the clay loam soil, which varied by about an order of magnitude. Overall, the trial-and-error die-off rates tended to be lower than the measured values used in the optimization procedure. Lower trial-and-error die-off rates appeared to have been compensated for by attachment rates that were approximately 2 to 3 times greater than those obtained by optimization. In the case of the sandy soil attachment rate, the results from both estimation methods resulted in identical outcomes. By contrast, detachment rates were 37 to 74 times lower than attachment rates for trial-and-error and Only the sand soil experiment produced sufficient breakthrough data to attempt a simulation of the virus data. The liquid phase die-off rate was about half of the solid phase die-off rate in case of the optimization procedure (R 2 =0.76), but more than an order of magnitude greater for the trial-and-error simulation. On the other hand, the trial-anderror virus attachment and detachment rates were very different to each other, attachment rates were 2 orders of magnitude higher than detachment rates. Overall, the results indicate that virus attachment rates were more than an order of magnitude higher than those for bacteria, while bacteria and virus detachment rates were similar.

DISCUSSION
The bacteria die-off rates measured for the three soil types were different (Table 1. 5), which provides evidence for the effect of local environmental soil conditions on bacteria die-off rates [44]. Chao and Feng [45] studied the survival of E. coli HB101 strains added to a silt loam soil at 30 0 C, resulting in die-off rates ranging from 0.04 d -1 to 0.20 d -1 (0.0017 hr -1 to 0.0083 hr -1 ). Powelson and Mills [46] reported E. coli die-off rates of 0.0259 hr -1 and 0.0693 hr -1 in sand columns under saturated and unsaturated conditions, respectively. E. coli isolated from STE collected from an OWTS near Lake Okareka, New Zealand, were investigated to elucidate microbial attenuation and transport through pumice sand aquifers [47]. The results of that study showed soil-attached E.coli die-off rates ranged from 2.59 hr -1 to 4.47 hr -1 . These few studies suggest that solid phase bacteria die-off rates have to be determined under environmental conditions representative of the location where the construction a new OWTS system will be built.
The measured die-off rates reported here for all three soil types may be different from in situ rates where the soils were collected (Colorado, Georgia, Rhode Island). It is also likely that the liquid phase die-off rates differ among locations because of differences in the chemical, physical and biological properties of wastewater. For the trial-and-error simulations, both the solid and liquid phase die-off rates were treated as variables, whereas they were fixed to the measured values during the optimization procedure (Table   1. 5). The attachment/detachment rates were fitting variables in both procedures. Based on the assumption of potentially location-specific solid and liquid phase die-off rates, treating these rates as variables may be considered for the optimization procedure. A better fit could be obtained by treating the die-off rates as variables, particularly during the tailing phase of each experiment, further research is needed to confirm this approach.
Average attachment rates, derived from either optimization or trial-and-error procedures, were highest for the sandy soil (0.163 h -1 for E. coli, 0.91 h -1 for MS-2). This result was unexpected because higher attachment rates are typically reported for fine-grained clay materials, rather than sandy soils. In general, the intrinsic lower surface area of coarser soils should result in less adsorption of microbes compared to finer textured soils [10,48,49]. In addition, the smaller pores that are prevalent in fine-grained soils are more effective for mechanical filtration (straining) of microbes than those in coarser porous soils. Conversely, unsaturated soils tend to retain more microbes than saturated soils.
That is, with decreasing water content, higher retention of bacteria and viruses in the soil has been observed [38,[50][51][52]. Because the air-water interface increases at decreasing water content, the removal and retention of microbes in fine-grained, such as the clay loam, should be, at a given water content, greater than in granular soils [38,52,53]. In our study, the water content of the sand and sandy loam soils at the end of the experiment was lower (0.15 and 0.23 g/g, respectively), compared to the clay loam soil (0.32 g/g; Table 2). Therefore, the higher air-water interface in the coarser soils could explain the higher attachment rates, since more water-free surface area is available to interact with the microorganisms. Measurements of the air-water interface area at different saturations in various soil materials would be necessary to confirm this proposition.
The effects of soil texture on microbial removal are expected to be different for bacteria and viruses. In our experiments, MS-2 phages were removed much more effectively than bacteria. Sandy loam and clay loam soils removed phages more extensively than sandy soil did. Two main mechanisms have been considered for pathogen removal in soil: (i) mechanical filtration and (ii) adsorption. For instance, Powelson et al. [13] investigated the fate and transport of a Salmonella phage in structured soils and found a reduction in virus concentration of about 60% to 90% in clay, clay loam and silt loam soils. In a review of the literature, Amador et al. [53] concluded that, although coarser textured soils tend to remove fewer bacterial pathogens than finer textured soils, the depth of treatment is important in order to obtain acceptable removal rates (close to 100%). The authors suggest that, because preferential pathways are more common in large-grained, textured soils, these pathways facilitate the transport of microbes to deeper depths relative to fine textured soils. In addition, they suggest that the soil texture and depth of soil treatment are not well-correlated variables in virus removal, which is consistent with the hypothesis that virus removal occurs by adsorption processes rather than mechanical filtration. Virus removal by adsorption processes is in agreement with our results, where the model computed higher attachment rates for viruses than bacteria in sandy soil.
On average, the detachment rate for both bacteria and viruses in all soils was 1.6% of the attachment rate (Table 1. 5). The lowest detachment rate values were observed in the structured clay soil, which suggest that bacteria and virus attachment in those soils is practically irreversible. Under those conditions, detachment can be considered negligible. This is consistent with previous studies, which concluded that the attachment of microbes to soil particles is an irreversible process [24,[54][55][56].

CONCLUSION
Modeling results showed only small differences between attenuation parameters (microbial attachment and detachment rates) obtained by optimization and trial-and-error simulation processes, i.e. results were generally within a factor of three of each other.
The microbe detachment rates were about two orders of magnitude lower than the corresponding attachment rates. Low or negligible detachment rates suggest quasiirreversible adsorption of microbes to soil. GFP E.coli concentrations were generally underestimated by the optimization simulation, whereas a better fit was achieved by the trial-and-error procedure, particularly for the tailing phase of each experiment. In case of the liquid and solid phase GFP E.coli die-off rates, the results of the optimization and trial-and-error simulations were generally within a factor of three of each other. Overall, the combination of lower die-off rates and higher attachment rates resulted in a better description of the tailing phase when using the trial-and-error procedure.
In general, the fit obtained in the optimization process should improve when concentration of bacteria or virus is measured more frequently. In addition, the results of the E.coli and MS-2 phage die-off rate experiments support the findings by Foppen and Schijven [44] that these measurements should be ideally collected under in situ conditions of the sample location rather than under standard laboratory conditions. This change in procedure would contribute to a better understanding of the effects of the local conditions on the soils and the resulting degradation/attenuation of those microbes.
The experimental data for the structured clay loam suggests that early breakthrough of the bacteria occurred. Although the presence of preferential flow pathways in the mesocosms likely influenced the results, it is not possible to simulate those conditions 27 with existing models. To better simulate the preferential flow effect on transport and fate of pathogenic contaminants in the soil it is necessary to evaluate the in situ spatial distribution of soil hydraulic properties. In the interim, a dual permeability model may be used to diversify the different flow patterns that might occur in the soil profile [43,57,58]. Numerical modeling limitations were also evident when simulating the transport of microbes because the model neglects processes that intervene in the attenuation of microorganism in the field (i.e., straining, size exclusion).

This study was funded by a grant from the Water Environment Research Foundation
where K r is the relative hydraulic conductivity and K s the saturated hydraulic where C and S are the (virus, bacteria) solution concentration [NcL where µ r is the values of a particular coefficient (rate constant) at the reference water content, θ r , µ is the value at the actual water content θ, and B is a solute dependent parameter (usually 0.7). The reference water content, θ r , which may be different for different soil layers, is calculated from the reference pressure head, h r , which is considered to be constant for a particular compound.
The die-off rates and transport parameters were first determined from mesocosm experiments and fitted using the inverse solution algorithm included in the HYDRUS model (Morales et al., 2014). These parameters were then imported into a model that simulated a conventional OWTS trench with intermittent dosing (Fig. 1A).    Tables S1 and S2).

Temperature dependence
HYDRUS accounts for temperature dependence of transport and reaction rates by

Calibration and Validation
The model was calibrated utilizing the bacteria transport data obtained in previous laboratory mesocosm experiments (

Soil Texture
We simulated how a trench system responds to sand, sandy loam and clay loam as native soil in the drainfield. The initial E. coli concentration in STE was 10 5 cfu mL -1 .
The results showed that the bacteria concentration was significantly reduced (99.99% reduction) in the first 30 cm of soil in the sandy and clay loam (Figure 2. 2a and b).
Similar results were observed for the sandy loam, except that a greater soil depth was required to reduce bacteria below 1 cfu/100 ml. Deeper bacteria movement occurs in the sandy loam due to lower solid phase die-off and attachment rates than those in the sandy and clay loam soils (Supplemental table S2). For example, the sandy soil has die-off and attachment rates that are 2.08 and 6.27 times higher, respectively than the sandy loam.
Hence, in sandy loam, fewer bacteria are removed and attached to the soil grains, which allows them to travel deeper through the soil profile.  Table S2), indicating that the soil particles have a greater level of physical interaction with microbes (due to a high specific surface area), and thus more bacteria are retained and removed on the particle surface.

Hydraulic Loading Rate
A range of 50% lower to 170% higher HLR was modeled (initial HL of 0.424 cm h -1 ) and the effluent E. coli concentrations were recorded at observation nodes located along two vertical profile cross-sections (Figure 2. 3a and b). At shallow depth (10 cm), and directly beneath the trench, the STA removed 10% more the HLR was low (0.212 cm h -1 ) relative to a higher HLR (0.72 cm h -1 ). At a greater depth (17 cm), more than 90% of bacteria were removed in all three soil types, independent of HLR. These results indicate that, although the initial removal percentage was lower when the HLR was higher, the rate of removal increased with depth. The relatively lower reduction in E. coli design HLR method that depends on the type of wastewater treatment system and soil textural class. For a conventional OWTS trench, our optimal design HLR value is slightly higher than those suggested by Siegrist (2007) for sand (0.167 cm h -1 ), sandy loam (0.083 cm h -1 ) and clay loam (0.021 cm h -1 ).

Presence of a Biomat
The For example, at 28-cm depth, more than 99% of E. coli influent concentration was removed when the hydraulic conductivity (K s ) was reduced 10 or 100 times. Removal increased by 9.5%, 12.0% and 2.6% in sandy, sandy loam and clay loam soils, respectively, relative to removal when the model was run with K s initial values or those generated by Rosetta lite ) and assigned to all three soils (Figure 2. 4). These results show that the presence of a biomat layer improved the performance of the STA in terms of E. coli removal. However, the increase in bacteria removal due to the biomat layer is relatively modest, and its benefits must be weighed against the potential consequences of excess clogging and hydraulic failure. Siegrist (1987) suggests that the biomat layer helps to reduce bacteria concentration by increasing the biogeochemical activity, straining, and promoting unsaturated conditions below the infiltration surface. Gerba (1975) showed that the highest bacteria removal rates occur between 2 cm and 6 cm below the infiltrative surface of the STA. These results are consistent with our data, which showed that E. coli influent concentration was reduced by >99% between 23.3 cm and 28 cm beneath the pipe and 5 cm away from the trench wall when the simulated biomat layer's hydraulic conductivity was reduced by 1 or 10 orders of magnitude relative to the initial K s . All three soils had greater removal rates at higher hydraulic conductivity values at observation profile points located 5 cm lateral distance from the trench sidewall (Figure 2. 4b). This is because, at a higher biomat hydraulic conductivity, STE no longer ponds on the biomat layer or the trench wall (Finch et al., 2008). As illustrated in a -c, flow around the biomat layer results in more treatment because the water flow is forced to pass over the sidewall trench and bacteria are transported through a longer path, which also results in more interaction with the soil matrix. However, when most of the STE infiltrates through the biomat, and the conductivity of the biomat is higher ( Figure   2.5 d) , any water flowing sideways from the trench must pass through the biomat wall layer, resulting in a higher bacteria concentration.

Precipitation events
The precipitation scenarios that were modeled to evaluate the influence of infiltrated precipitation on E. coli removal in the three soils are summarized in Table 2. 3.
For the sand and sandy loam soils, no surface runoff was observed during the simulation of any precipitation events, indicating that applied rainwater was infiltrated completely.
However, results could not be obtained for the clay loam soil, because when the Decreased bacteria removal rates in response to increasing amounts of rainfall may be due to the development of near-saturated or saturated flow conditions, which occur temporarily during rain events (Table 3). This is because bacteria survival is greater in moist soil than in dry soil (Campbell and Beiderbeck, 1976;Kibbey et al., 1978). In the model, this phenomenon is caused by soil water content variations and, as a result, the die-off rates are affected (i.e. water-content dependence of die-off and attachment rates, Eq. 3 and 4). Furthermore, soils exposed to prolonged dry periods -and consequently  Gerba (1975) found that low temperatures support the survival of enteric pathogenic bacteria for months or even years. Some researchers (Shaw, 1970;Fletcher, 1977)

MODEL IMPLICATIONS ON BACTERIA REMOVAL RATES
The removal of bacteria is influenced by the variable environmental and operational conditions assumed for each of the simulated scenarios for the conventional OWTS. The changes in removal rates are most evident, specifically, for the first 23 cm below the distribution pipe. Higher removal rates were computed by the model (at shallower depths) when the hydraulic loading rate was lowered 50% and the wastewater infiltrated in a clay loam soil. These results explain the importance of soil texture and flow rates for system design. A simulated biomat also improved the bacteria removal percentages due to a lower hydraulic conductivity or clogged soil pores on the surface (modeled biomat growth). This is consistent with studies that showed a higher removal efficiency of bacteria in clogged soil treatment areas or sand filters compared to unclogged systems . The modeled precipitation event scenarios did not cause significant changes in the model outputs or removal rates, and no variation was observed in the OWTS performance. A higher rainfall intensity needs to be applied over the soil surface to ensure that more bacteria are detected on the effluent concentration because of water saturation.
All of the modeled scenarios and conditions may be considered as OWTS performance evaluation. Our results can help to define system design (i.e., size and type of system) by incorporating data on wastewater, soil physical/chemical properties and site properties. Inappropriately designed or failed OWTS are sources of surface and groundwater contamination, which present a serious public health risk (US EPA, 2002).
Bacteria are of great concern because they can be transported for long distances in water bodies, causing illness through body contact or ingestion of contaminated water.

CONCLUSIONS
Although perhaps considered a shortcoming, we do not consider the lack of model suggests that the effectiveness of the STA will increase as the average soil temperature rises. Our findings also identified a role for soil texture in E. coli reduction, with finer textured soils removing more bacteria than coarser textured soils. The simulation of variable stress conditions suggests that environmental and operational factors influence the performance of soil-based wastewater treatment, and that this treatment will likely respond to changing temperature and precipitation patterns predicted by climate change models. Which of these factors becomes more influential, and how these factors correlate with other environmental or operational factors not considered in this study, remains to be evaluated.  10 cm 1 cm h -1 for five hours on day 1 followed by2.5 cm h -1 for two hours on day 12.

TABLES
15 cm 1 cm h -1 for five hours on day 1 followed 0.5 cm h -1 for ten hours on days 4 followed by 2.5 cm h -1 for two hours on day 12.
25 cm 1 cm/hr for five hours on day 1 followed 0.5 cm h -1 for ten hours on days 4 and 8 followed by 2.5 cm h -1 for two hours on day 12.

cm
1 cm h -1 for five hours on day 1 followed 0.5 cm h -1 for ten hours on days 2, 4, 6, 8, and 10 followed by 2.5 cm h -1 for two hours on day 12. 35 cm / HLR As scenario "35 cm" but with HLR increased 1.7 times (from 0.424 cm h -1 )        total suspended solids (TSS), pathogens and nutrients (i.e. N, P). However, these systems are not designed for removal of nitrogen (N) [5,6] or emerging organic contaminants, such as personal care products and pharmaceuticals [7,8]. Furthermore, their use is limited in areas where a shallow water table lies beneath the STA, as well as in many coastal areas. Advanced OWTS are used in areas that are at risk of water use impairments (i.e., pathogen and nutrient contamination) because of a shallow-placed infiltrative surface.

E. coli
A conventional OWTS consists of septic tank, distribution box and a gravity-dosed STA, which treats septic tank effluent (STE) as it infiltrates and percolates through the soil.
The STA has a pipe-and-stone (P&S) configuration: a horizontal drain constructed from perforated pipes located in an excavated trench backfilled with gravel or crushed stone.
Advanced OWTS integrate engineered treatment units (i.e., sand filters) that provide additional treatment. The STE can then be pressure-dosed to a type of STA, known as a pressurized shallow narrow drainfields (PSND). In advanced and conventional OWTS, the STA is dosed with STE or advanced-treated effluent (ATE), and is usually installed 15 -30 cm and ~ 60 cm below the ground surface, respectively [9]. The shallow depth in the STA of advanced OWTS increases the vertical separation distance, or unsaturated zone, and enhances the potential for treatment before the effluent reaches the water table [10][11][12]. A thicker unsaturated zone increases the opportunity for O 2 diffusion and attenuation of contaminants [13][14][15][16]. There are other advantages of PSND relative to conventional STAs. For example, pressurized systems disperse the effluent more uniformly over the STA, which avoids overloading (ponding) and supports complete infiltration [17]. A shallow drainfield also enhances the transformation of nutrients by microorganisms and their uptake by plants because effluent distribution takes place closer to the soil surface, within the root zone, where microbial activity is highest [11].
OWTS can be sources of surface and groundwater contamination and they are one of the top 10 probable sources of impairments in rivers, lakes, and coastal shoreline in U.S. [18]. Pathogens and nutrients are frequently cited causes of impairments in water bodies.
Nitrogen is of particular concern because its presence in high concentrations may stress the functioning of surface and coastal water ecosystems. Approximately 32% of stream length have been reported to be stressed or affected by N in U.S. [19][20][21]. Excess N in coastal areas and some freshwater ecosystems can result in eutrophication, decreased dissolved oxygen levels and habitat degradation [19][20][21]. N in wastewater is found as organic nitrogen, ammonium (NH 4 + ), nitrate (NO 3 -) and nitrite (NO 2 -) [22]. The nitrogen speciation in OWTS effluent is dependent on the type of treatment processes. In conventional systems, the STE is typically composed of 10-30% organic nitrogen and 70-90% NH 4 + [9,23]. The STA of advanced systems receives effluent from an advanced treatment system (ATE) such as a single-pass sand filter, where the concentration of NH 4 + is reduced and converted to NO 3 -. Therefore, N speciation in ATE is 18% organic N, 26% NH 4 + and 56% NO 3 - [24].
As STE and ATE are loaded to the drainfield, N species can be transformed or removed in the soil below the infiltrative surface. Nitrogen transformations in conventional and advanced STAs have been studied to some extent [24,25]. Nitrification and denitrification are thought to be the main processes that contribute to N speciation in the drainfield [26]. In nitrification, NH 4 + is oxidized by autotrophic bacteria to NO 3 in the STA under aerobic conditions. Nitrate can be subsequently reduced by heterotrophic denitrifying bacteria to nitrogen gas (N 2 ) or nitrous oxide (N 2 O), which results in net removal of N from wastewater.
The fate and transport of N in OWTS drainfield is a complex process controlled by many factors, including pH, temperature, moisture content, carbon availability, and oxygen diffusion. Computer-aided numerical models have been developed to understand N dynamics in in the STA. A broad variety of models have been used that include OWTS as a N source, but most of these only simulate NO 3 transport groundwater, but hydrodynamic processes (advection-dispersion) are not included [25][26][27].
Other researchers have used HYDRUS 1D, 2D and 3D models to predict the fate and transport of N in OWTS [28][29][30][31]. HYDRUS is a commercially-available computer program used to simulate water flow, solute and microbial transport [32], heat transport, and colloid transport in variably-saturated porous media [33,34]. For instance, Hassan [28] used HYDRUS 2D to simulate an onsite wastewater subsurface drip irrigation system (SDIS) dosed with pre-treated wastewater in a sequential batch reactor (SBR) .
The wastewater was collected from a restaurant and contained oil and grease with high organic matter content. Together with a grease trap and aeration unit, the SBR was used as a pre-treatment unit, where NH 4 + was nitrified and entered the SDIS as NO 3 --N. The model included NO 3 transport, plant uptake, and denitrification in order to estimate an N mass balance for the SDIS-SBR system. In addition, soil water pressure head data was collected and modeled. Based on this model, it was estimated that 48% of NO 3 was stored in the soil profile, 27% was taken up by plants, 22% removed by denitrification, and 0.4% NO 3 left with the drainage water.
Heatwole and McCray [29] used HYDRUS 1D to model fate and transport of N in a conventional STA. The model was developed to evaluate the concentration of NO 3 reaching groundwater using site-specific data and input transport parameters estimated from statistical distributions. The results showed that no NH 4 + was detected at 30-cm depth below the infiltrative surface or deeper in the model domain. Also, NO 3 concentrations were predicted to be below maximum contaminant level (MCL = 10 mg N/L) when the median value for denitrification rate was applied.
HYDRUS 2D/3D was used to fit experimental soil pressure head and N and chloride (Cl -) data collected from a conventional OWTS with a drainfield installed in a clay soil [30]. Geomat. We determined nitrification and denitrification rate coefficients for the three drainfield types and used this to estimate N losses from simulations and compared to actual experimental data. The information obtained from these models is expected to aide designers of OWTS and regulators to make informed decisions about the most effective treatment practicse for removal of N species in the STA.

Experimental setup
Replicated mesocosms (n = 3) were engineered to mimic the soil treatment area and   [24].

Modeling approach
HYDRUS 2D/3D version 2.0 was used to simulate water flow and solute transport in soils under variably-saturated conditions. The HYDRUS program numerically solves the Richards equation for saturated-unsaturated water flow (Eq. 1): where K r is the relative hydraulic conductivity and K s the saturated hydraulic conductivity [LT -1 ].
HYDRUS allows the user to select among several analytical models to describe the soil water retention and unsaturated hydraulic conductivity functions. In our model, the van Genuchten [39] equation was applied to compute the soil hydraulic properties (Eq. 3-5): where α (L -1 ), m (dimensionless), and n (dimensionless) are fitted parameters, θ(h) is the volumetric water content (L 3 L -3 ), θ s is the saturated volumetric water content (L 3 L -3 ), and θ r is the residual volumetric water content (L 3 L -3 ). The unsaturated hydraulic conductivity function K(h) (LT -1 ) is written as follows: where m = 1-1/n and l is the pore connectivity parameter, which it is assumed to be about 0.5 [40]. The model permits the application of the convection -dispersion equation in the liquid phase to simulate solute transport and fate. Chemical equilibrium and linear adsorption is described by the following mass balance equation: where c is dissolved solution concentration [ML −3 ], t is time (T), K d is the adsorption coefficient (L 3 M -1 ), μ represents the solute transformation or degradation rate in the liquid phase, x is the solute travel distance (L) and z is depth (L). D w ij is the dispersion

Model domain and boundary conditions
The model domain was developed to resemble the engineered mesocosm columns, not only physically but also in terms of operational conditions. The geometry of the domain properties reproduced the two shallow and trench drainfields described previously [24]. The native soil, used for the mesocosms, was described as a Bridgehampton silt loam (coarse-silty, mixed, active, mesic Typic Dystrudept) (S1 Table). The infiltrative surface was placed 20 -25 cm below the ground surface for PSND and GEO (A horizon), and 84 cm (C horizon) for P&S. Based on field observations, two layers were used to simulate B (gravelly loamy sand) and C (gravelly coarse sand, 40 -45% gravel) horizons. For the purpose of this study and because of their similarities in the particle size distribution, sublayers B w and 2B w were assumed to be B horizon and modeled as one single layer.
Finite element mesh with a maximum element size of 3.90 cm was generated automatically with 478, 537 and 614 nodes for P&S, GEO and PSND, respectively. A denser grid was defined around the simulated distributed pipes and the PVC cover.
Elements size in that area was 0.45 cm. Observation nodes were located along the soil profile to compare the observed against modeled data. Two observation nodes were placed 15 cm and 30 cm below the infiltrative surface and one was located at the bottom of the model domain and one at the column outlet.
Atmospheric boundary condition was assigned to the top of the columns (Figure 3. 2).
The sides and bottom of the column were treated as no-flux boundaries. As wastewater infiltrates, it accumulates on the bottom and flows out when the soil is saturated or a hanging water table is formed. In order to account for this condition, a seepage face boundary was selected for one of the nodes at the bottom right of each soil column ( Figure 3. 2). In the HYDRUS model, this assumption is that the water is removed by overland flow when saturated conditions prevail [33].

N transformation modeling
Nitrogen losses in STA are attributed to NH 4 + conversion to NO 3 or nitrification followed by reduction of NO 3 to N 2 O or N 2 through denitrification. Therefore, we developed a decay model to simulate the N species fate and transport in conventional and advanced STA in which N was assumed to be transformed as follows [26]: where μ is described as the zero-order reaction rate for nitrification. N species were modeled using sequential decay reactions built into HYDRUS [41]. In this approach, the program provides nonlinear non-equilibrium reactions (adsorption-desorption) between the solid and liquid phases (soil-water interface) based on the two-site sorption concept [42,43]. It is considered that the sorption sites are composed of two fractions, sorption in one of the fractions is assumed to be instantaneous, while on the remaining site is timedependent. Also it is assumed that the solute transport takes place by convection and dispersion. The measured total N (TN) was modeled as an input concentration to include all N infiltrated in the drainfield. Thus, the influent organic N was considered to be transformed to NH 4 + through ammonification.
Several researchers have reported the water content dependency of nitrification and denitrification [44,45]. Nitrification is an aerobic process that occurs at low soil water content because high soil water content increases tortuosity and, as a result limits oxygen diffusion and the activity of nitrifying bacteria [46]. On the other hand, denitrification takes place under soil-saturated conditions which promotes anoxic conditions.. Thus, HYDRUS was modified to account for the effect of soil water content and aeration conditions on N transformation on OWTS. A water content dependency function was built in HYDRUS that allows computing of nitrification and denitrification rates at low water saturation or unsaturated conditions. The program incorporates the water content dependency function implemented in DRAINMOD-N2 [47], an agricultural computer program used to model N transformation and the impact of water content. DRAINMOD-N2 simulates nitrification and denitrification using Michaelis-Menten kinetics [48]. For nitrification the model uses a stepwise function to model the influence of nitrification inhibitors on decay rates. Denitrification is modeled as a function of the organic content decrease with depth [49]. The following expression describes the nitrification rate: , , where μ nit is the calculated nitrification rate, μ nit,max is the maximum nitrification rate, C NH4 is the ammonium-nitrogen concentration, and K m,NH4 is the half-saturation constant, which is the ammonium-nitrogen concentration at which the nitrification rate is half its maximum value. The value of f sw is soil-water content dependency functions (Eqn. 10): where fsw varies between 0 and 1. The term, fs is the value of fsw at full saturation, fwp is the value of fsw at the wilting point, S is the water-filled pore space (or relative saturated water content), Sh is the upper saturation boundary for optimal nitrification, S l is the lower saturation boundary for optimal nitrification, swp is the saturation level at the wilting point, and e 2 and e 3 are fitting exponents. The denitrification rate equation included in the modified HYDRUS version is written as follows: where μ denit is the denitrification rate, μ denit,max is the maximum denitrification rate, C NO3 is the nitrate-nitrogen concentration, and K m,NO3 is the half-saturation constant, which is the nitrate-nitrogen concentration at which the denitrification rate is half its maximum value.
The terms f t , and f z are temperature-dependency, and carbon dependency functions, respectively.
, 0 1 (14) f t varies between 0 and 1, T is the temperature, T opt is the optimum temperature for nitrification, β and a are fitting parameters, and z is depth below the infiltrative surface.
The term f sw,dn is the water content-dependency function, s dn is a threshold saturation value for denitrification, s is the actual soil saturation, and f is a fitting exponent.

Calibration and parameter sensitivity
Model calibration was carried out to determine input parameter values for obtaining the best fit between the predicted and measured soil data. The model was calibrated by coupling HYDRUS with UCODE, a computer program used to estimate parameters through inverse modeling by nonlinear regression [50]. The nonlinear regression problem is solved by minimizing a weighted least squares objective function with respect to the parameter values using a modified Gauss-Newton method.
A sensitivity analysis in UCODE was performed to identify which of the parameters influenced the model output results and their uniqueness. Composite scaled sensitivities (CSSs) were calculated to identify the influence of the observed data on the estimation of a parameter. CSS is the measure of the total amount of information provided by the observations to estimate one parameter. Larger CSS values indicate that those parameters are likely to be estimated more precisely with the proposed model and observations. The ratio of the CSS of a parameter to the maximum CSS was used to compare relative sensitivity among estimated parameters. Parameters with CSS ratio less than 0.01 are not sensitive and denote that a regression will not converge. Therefore, in some cases, parameters with CSS ratio < 0.01 were excluded from the inverse modeling process. Fitting parameter values were assigned to the entangled plastic filaments (GEO) and crushed stone (P&S) systems. Both materials were considered highly-conductive (K s = 3,000 cm day -1 ) with low porosity and residual water content that was similar to a coarse gravel soil. Initial parameter values for native soils were estimated using Rosetta [51] and fitted with UCODE, whereas values for the plastic filaments and gravel layers were kept fixed. Initial N transformation rates were selected from McCray [26] and initial NH 4 + and NO 3 soil concentration were set to zero. Water dependency function parameters were selected from McCray et al. [52]. Finally, the model was run for 3-months (90 days). The predicted N species concentrations were computed to estimate a N balance produced by each of the three OWTS.
The best fit between the predicted and observed data were evaluated the root mean squared error (RMSE) (eq. 15). (15) where ŷ i is the predicted value, y i is the observed value, and n is the number of observations. A RMSE value closer to zero indicates the best of fit to observed data.

Water content
The model was calibrated using soil moisture data to simulate the unsaturated soil profile beneath the infiltrative surface, and to account for moisture changes associated with N transformation processes. Given that variations in water content around the measured moisture data were minimal, the mesocoms simulations were under steady state conditions. The soil hydraulic parameters (θ r, θ s, α, K s, n and l) were determined for each of the soil layer (silt loam and gravelly-coarse sand); only the pore connectivity parameter value was not calibrated or changed (l was equal to 0.5, as recommended [33]).
Ten parameters were calibrated for the advanced OWTS technologies and five for the conventional one. In advanced STAs, the measured water content (cm -3 cm -3 ) ranged from 0.11 to 0.13 and 0.02 to 0.05 at 15 cm and 30 cm below the infiltrative surface, respectively. Even though the intact soil cores were collected in close vicinity to each other, water content variations were expected at greater because of the increasing influence of variable physical properties on soil moisture and water flow with depth.
Also, the amount of water retained in the upper soil layer was expected to affect the hydraulic properties of the deeper soil layers. This more heterogeneous behavior of the soil system is illustrated for one of the three PSND mesocosms that showed higher water content (0.23 cm -3 cm -3 ) at the 15-cm depth compared to the other (0.11 to 0.13 cm -3 cm -3 ).
These variations are indicative of soils with low residual and high saturated water content characteristics.
Overall, the model results showed a good fit between the observed and simulated water content data for PSND, GEO and P&S (Figure 3.  was computed for the P&S drainfield mesocosms (Table 3. 1). It is most likely that a biomat developed over time above the infiltrative surface, which provides unsaturated conditions and a reduced saturated hydraulic conductivity (Figure 3. 4).

Sensitivity analysis
For the sensitivity analysis of the PSND and GEO, five soil hydraulic parameters (θ r, θ s, α, K s and n) for each of the two horizons (silt loam and gravelly-coarse sand) were calibrated simultaneously (10 parameters total). For P&S, the moisture data for the gravelly-coarse sand was calibrated with the 5 parameters mentioned above. The sensitivity of soil moisture to soil hydraulic properties for each of the mesocosms and soil horizons is shown in Table 3. 2. Most of the selected parameters were significantly sensitive (CSS ≥ 0.01) to the water content data. In most advanced STAs, the silt loam soil properties were sensitive to the simulated soil moisture. Not unexpectedly, the soil properties were found most important for the calibration of the hydraulic parameters along the soil profile. Generally, the most sensitive parameters were θ ss and n s .
Conversely, K ss , θ rs , θ rg and K sg were not significant or least sensitive parameters to the moisture data. For P&S, the saturated and residual water content (θ rg and θ sg ) were very important parameters determining the soil moisture distribution along the profile. Also, the hydraulic conductivity (K sg ) (range: 908.88 to 942.48 cm day -1 ) was more sensitive compared to PSND and GEO (CSS = 0.21 to 0.25).
In one of the PSND columns (Table 3. 2, column #3) the K sg was not a sensitive parameter to the fitted water content data (CSS < 0.01). In this mesocosm, the water content of the silt loam was almost two times higher (0.23 cm -3 cm -3 ) than those values observed for the other two PSND columns (0.11 cm -3 cm -3 to 0.13 cm -3 cm -3 ). These variations are likely linked to soil heterogeneities and affected the sensitivity of K sg as reflected in the model output data.

Nitrogen transport and fate
Nitrification and denitrification were modeled using a water content-dependent function to account for changes in oxygen diffusion and availability in the mesocosms.
The function uses water-filled pore space or relative saturation to mimic soil aeration during water infiltration [49]. Based on this approach, NO 3 production is achieved with a water-filled pore space (WFPS) of 0.20 and the maximum nitrification rate is reached when WFPS is more than 0.35. Denitrification takes place when WFPS is more than 0.60 and the highest N 2 gas production is observed at saturation (WFPS = 1.00) [55,56]. Linn and Doran [56] reported that organic carbon decomposition associated with N mineralization and immobilization occurs when WFPS ranges from 0.5 to 0.6 and near saturation as well. Therefore, WFPS variation may affect the denitrification rates in the soil drainfield. However, it must be emphasized that the aqueous solution used in those experiments [55,56] had a higher dissolved oxygen concentration compared to the STE and ATE used in this study. These observations show that the relationship between WFPS and relative rate of microbial nitrification and denitrification may be affected during N transformation, and nitrification and denitrification may occur at lower WFPS.
The nitrification and denitrification rate coefficients were computed using Eq. 9 through 14, and parameter values were selected from literature data [52]. The fitted parameter values for the water-content dependent transformation rates are shown in Table 3. 3.
Initially, the model was adjusted until the best fit between the observed and predicted data was achieved. As a result, the parameters for nitrification and denitrification dependency functions are median values that best reproduced the observed data [52].
The fitted water content was important to elucidate the N transformation and decay in the mesocosms and the application of the water content dependent functions.
The results showed that the WFPS was higher than 0.27 (P&S gravelly-coarse sand) in all drainfields types ( Nitrate tended to increase with depth along the soil profile in all mesocosms, with the highest concentration detected in the outlet (seepage boundary). For ATE, the model output included NO 3 already in the influent water as well as NO 3 produced in situ from NH 4 + conversion. In PSND and GEO, influent total N included NO 3 and NH 4 + . Some of the nitrate resulted from NH 4 + being nitrified in the sand filter that preceded the treatment system from which the ATE was collected from. The model suggests that the remaining NH 4 + will be transformed to NO 3 in the drainfield.
The predicted NO 3 concentrations showed an acceptable goodness-of-fit with the observed data, with RMSEs that ranges from 4.45 mg L -1 to 9.65 mg L -1 in all STA types ( Figure 3.6). Lower RMSE values were observed for predicted NO 3 data for PSND and GEO compared to P&S. The ATE was assumed to be more uniformly distributed over the infiltrative surface in the PSND and GEO in the absence of an overlying layer (i.e., crushed stone) that influences the water flow and solute transport.

Nitrification and denitrification rates
The processes involved in N transformation and removal are mainly nitrification and denitrification. In addition, NH 4 + sorption to soil can affect the fate and transport of N in some OWTS drainfields. Because of the low sorption capacity of the soils used in drainfield mesocosms (Supplemental Table S1), NH 4 + sorption was not simulated in this model. Therefore, all NH 4 + moves with soil water and can be readily nitrified. Average simulated nitrification and denitrification zero-order reaction rates were computed to analyze the N dynamics and conversion in all drainfield types (Table 3.  Relative to each other, denitrification rates were higher in P&S than GEO and PSND. This finding was consistent with the experimental results presented in [24], where denitrification was higher in P&S compared to the other STAs. Besides anaerobic conditions, denitrification requires organic carbon to proceed [52]. Because ATE has a low organic carbon content, it may have limited the extent of denitrification in the advanced drainfield mesocosms. This is consistent with [24].

N losses and comparison between simulated and real systems
Average modeled N losses were calculated and compared with the experimental data from all of the advanced and conventional drainfield mesocoms. The calculations were based on the 90-day simulation period and accounted for all N species produced. An N mass balance was calculated from the modeled N species for influent and effluent water. In P&S, the modeled effluent N was comprised of dissolved NO 3 -(82.72%) and were underestimated by the model. It is likely that not all organic N was converted to NH 4 + and as a result, less NH 4 + was nitrified (for our modeling approach, it was assumed that organic N has been completely transformed to NH 4 + before entering the treatment system). Organic N was found to account for 14% to 16% [24] of the total N in the effluent water in P&S, which is a significant amount for N loss . Also, a fraction of the influent organic N is likely non-biodegradable or recalcitrant (not amenable to ammonification), which means it might not be transformed in the treatment system, passing through the drainfield unchanged. For GEO and PSND, the modeled N losses occurred mostly as NO 3 -(90.75% and 88.45%, respectively). No significant amount of NH 4 + was observed during the 90-days simulation period (ranging from 0.23 to 1.41% for all drainfield types). Nitrogen losses as N 2 were more evident in P&S compared to the advanced technologies.

CONCLUSIONS
A model was developed to predict the N fate, transport and transformation in a conventional P&S drainfield and in two types of shallow narrow drainfields (PSND and GEO). The model was calibrated using water content, NH 4 + and NO 3 concentration data.
From these inputs, water flow and solute transport parameters were determined.
Nitrification and denitrification rates were computed as function of the soil water content and the WFPS. The model was capable to determine nitrification and denitrification zeroorder rates with acceptable goodness-of-fit between the observed and simulated data.
These results allowed quantification the N losses in all OWTS drainfield types and an estimation of the N species fluxes. This information is useful to better understand the N transport and transformation mechanisms and to identify potential contamination sources of groundwater.  Table 3. 2 Composite scale sensitivity ratios to the measured soil moisture data for the silt loam and gravelly-coarse sand soils for PSND, GEO and P&S.   Table 3. 5 Average zero-order nitrification and denitrification rates for the selected soils and materials in advanced and conventional drainfield mesocosms. Values are means ± SD (n = 3).

CONCLUSIONS
This study demonstrates that HYDRUS 2D/3D is a useful tool to predict the fate and transport of microbial and nutrient contaminants under different operational and environmental conditions. The model was able to estimate microbial (bacteria and virus) transport in an onsite wastewater treatment system (OWTS) soil treatment area and attachment-detachment rate coefficients were determined to better understand the transport parameters that control the pathogen concentration in the porous media. The effects of variable environmental conditions on OWTS performance was evaluated in a simulated OWTS trench. For instance, warmer soil temperature and light-to-heavy rainfall events affected the transport of bacteria, which indicates that climate change may influence the OWTS performance, particularly in the soil treatment area. In addition, our model predicted nitrification and denitrification rates and N losses in a conventional and two advanced OWTS, using measured NH 4 + , NO 3 and water content data of the soil matrix as calibration parameters. Nevertheless, our model proved that OWTSs are not able to remove nitrate in the STAs. Therefore, the need for additional design modifications in all OWTS types may be necessary to comply with the water quality standards established for nitrate (10 mg L -1 ). Finally, the modeling approach presented in this study will help to predict the extent of contamination and their spatial distribution which aides in identifying non-point sources, and to establish total minimum daily loads (TMDLs).