Search tips
Search criteria 


Logo of plosonePLoS OneView this ArticleSubmit to PLoSGet E-mail AlertsContact UsPublic Library of Science (PLoS)
PLoS One. 2016; 11(6): e0158292.
Published online 2016 June 29. doi:  10.1371/journal.pone.0158292
PMCID: PMC4927103

Modeling Nitrogen Losses in Conventional and Advanced Soil-Based Onsite Wastewater Treatment Systems under Current and Changing Climate Conditions

Zhi Zhou, Editor


Most of the non-point source nitrogen (N) load in rural areas is attributed to onsite wastewater treatment systems (OWTS). Nitrogen compounds cause eutrophication, depleting the oxygen in marine ecosystems. OWTS rely on physical, chemical and biological soil processes to treat wastewater and these processes may be affected by climate change. We simulated the fate and transport of N in different types of OWTS drainfields, or soil treatment areas (STA) under current and changing climate scenarios, using 2D/3D HYDRUS software. Experimental data from a mesocosm-scale study, including soil moisture content, and total N, ammonium (NH4+) and nitrate (NO3-) concentrations, were used to calibrate the model. A water content-dependent function was used to compute the nitrification and denitrification rates. Three types of drainfields were simulated: (1) a pipe-and-stone (P&S), (2) advanced soil drainfields, pressurized shallow narrow drainfield (PSND) and (3) Geomat (GEO), a variation of SND. The model was calibrated with acceptable goodness-of-fit between the observed and measured values. Average root mean square error (RSME) ranged from 0.18 and 2.88 mg L-1 for NH4+ and 4.45 mg L-1 to 9.65 mg L-1 for NO3- in all drainfield types. The calibrated model was used to estimate N fluxes for both conventional and advanced STAs under current and changing climate conditions, i.e. increased soil temperature and higher water table. The model computed N losses from nitrification and denitrification differed little from measured losses 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. Losses as N2 were estimated to be 10.4% and 9.7% of total N input concentration for SND and Geo, respectively. The highest N2 losses, 17.6%, were estimated for P&S. Losses as N2 increased to 22%, 37% and 21% under changing climate conditions for Geo, PSND and P&S, respectively. These findings can provide practitioners with guidelines to estimate N removal efficiencies for traditional and advanced OWTS, and predict N loads and spatial distribution for identifying non-point sources. Our results show that N losses on OWTS can be modeled successfully using HYDRUS. Furthermore, the results suggest that climate change may increase the removal of N as N2 in the drainfield, with the magnitude of the effect depending on a drainfield type.


Decentralized wastewater treatment systems, such as onsite wastewater treatment systems (OWTS), are engineered technologies for wastewater management practices that protect public health and lower contamination risk. Onsite wastewater treatment systems integrate a septic tank, where solids removal takes place, and a soil treatment area (STA), or drainfield, where contaminants are attenuated and treated wastewater infiltrates to safely recharge groundwater. Passage of wastewater through the STA of conventional system attenuates the 5-day biochemical oxygen demand (BOD5), total suspended solids (TSS), pathogens and nutrients (i.e. N, P). Conventional systems are not designed for removal of nitrogen (N) [1,2] or emerging organic contaminants, such as personal care products and pharmaceuticals [3,4]. Furthermore, they are less effective in areas with a shallow water table, and in 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.

Most conventional STAs receive septic tank effluent (STE). These systems have a pipe-and-stone (P&S) configuration: a horizontal drain constructed from perforated pipes placed in an excavated trench backfilled with gravel or crushed stone. Advanced OWTS integrate engineered treatment units (i.e., sand filters, aeration units) that provide additional treatment. The advanced treated effluent (ATE) can then be pressure-dosed to an STA place closer to the surface than that of a conventional OWTS. Such a system is known as pressurized shallow narrow drainfield (PSND). In advanced and conventional OWTS, the STA is dosed with STE or ATE. These drainfields are usually installed 15–30 cm below the ground surface (bgs) for advanced OWTS and ~ 60 cm bgs for conventional systems [5]. The shallow depth of the STA of advanced OWTS increases the vertical separation distance, or thickness of the unsaturated zone, and enhances the potential for treatment before the effluent reaches the water table [68]. A thicker unsaturated zone also increases O2 diffusion and attenuation of contaminants [912]. Another advantage of PSND over conventional STAs is that pressurized systems disperse the effluent more uniformly over the treatment area, which avoids overloading (ponding) and promotes complete infiltration [13]. 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 [7].

Onsite wastewater treatment systems are sources of surface and groundwater contamination. They are one of the top 10 probable sources of impairments of rivers, lakes, and the coastal shoreline by pathogens and nutrients in the U.S. [14]. Nitrogen from OWTS is of particular concern because its presence in high concentrations has a negative impact on surface and coastal water ecosystems. For instance, approximately 32% of U.S. rivers and streams are considered stressed or affected by N [1517]. Excess N in coastal areas and some freshwater ecosystems can result in eutrophication, followed by decreased dissolved oxygen levels and habitat degradation [1517]. Nitrogen in wastewater is present primarily as organic nitrogen, ammonium (NH4+), and nitrate (NO3-) [18]. The nitrogen speciation in OWTS effluent depends on the type of system. In conventional systems, the STE is typically composed of 10–30% organic nitrogen and 70–90% NH4+ [5,19]. In ATE the typical N speciation is 18% organic N, 26% NH4+ and 56% NO3- [20].

As STE and ATE are dispersed to the drainfield, N species can be transformed or removed in the soil below the infiltrative surface. Nitrogen transformations in conventional and advanced STA have been studied to some extent [21]. Nitrification and denitrification are thought to be the main processes that contribute to N speciation in the drainfields [22]. Nitrification involves the oxidation of NH4+ by autotrophic bacteria to NO3- under oxic conditions. Nitrate can be subsequently reduced under anoxic conditions by heterotrophic denitrifying bacteria to nitrogen gas (N2) or nitrous oxide (N2O), both of which result in net removal of N from wastewater.

The fate and transport of N in OWTS drainfield is a complex process controlled by temperature, moisture content, carbon availability, oxygen diffusion and other factors. A suite of computer-aided numerical models have been developed to simulate N dynamics in the STA. However, many of these only modeled the fate of NO3- in groundwater and hydrodynamic processes (advection-dispersion) [2123]. Others used HYDRUS 1D, 2D and 3D models to predict the fate and transport of N in OWTS [2427]. HYDRUS is a commercially-available computer program used to simulate water flow, solute, heat and colloid and microbial transport in variably-saturated porous media [2830]. For instance, Hassan [24] used HYDRUS 2D to simulate an onsite wastewater subsurface drip irrigation system (SDIS) dosed with advanced-treated wastewater in a sequential batch reactor (SBR). The wastewater was collected from a restaurant and contained oil and grease with high organic C content. Together with a grease trap and aeration unit, the SBR was used as a pre-treatment unit, where NH4+ was nitrified and entered the SDIS as NO3-. The model included NO3- transport, plant uptake, and denitrification in order to estimate an N mass balance for the SDIS-SBR system. In addition, the soil water pressure head data was modeled. It was estimated that 48% of NO3- was stored in the soil profile, 27% was taken up by plants, 22% removed by denitrification, and 0.4% NO3- left with the drainage water.

Heatwole and McCray [25] used HYDRUS 1D to model the fate and transport of nitrogen species in a conventional STA. The model evaluated the concentration of NO3- reaching groundwater using site-specific data. The authors relied on input transport parameters estimated from statistical distributions of first-order nitrification and denitrification rate coefficients, dispersivity, effluent loading rate and nitrogen effluent concentration. The model predicted that no NH4+ should be detected at or below 30-cm from the infiltrative surface. Also, NO3- concentrations were predicted to be below maximum contaminant level (MCL = 10 mg NO3-−N/L) when a median value for denitrification rate (0.042 day-1) was applied.

HYDRUS 2D/3D was also 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 [26]. The model involved the application of an N transformation chain, that is non-equilibrium transport of N in sequential decay reactions (NH4+ → NO3-→ N2) with water content-dependent, first-order transformation rates for nitrification and denitrification. Unlike Heatwole and McCray [25], the model assumed that N decay occurs and aquifer recharge was considered. The authors computed N losses from the STA with the calibrated model. Based on a N mass balance, the model predicted that 52% of N was removed by denitrification. Furthermore, less than 5% of N loss occurred as plant uptake and change in N storage. The model [26] was then used by Radcliffe and Bradshaw [27] to evaluate OWTS hydraulic loading rates (HLR) and N transformations in 12 soil textural classes. As in the previous study [26], water flow and N and temperature dynamics were simulated in a 2-D drainfield trench for a two-year period. All HLRs values (range: 1.48 to 5.40 cm d-1) were found to be suitable for all soil types except for the sandy clay textural class, where the trench was overloaded (HLR = 1.48 cm d-1). The predictions for denitrification losses varied widely among soil types, from 1% in sand to 75% in sandy clay. Leaching losses of NO3- were more significant than denitrification, ranging from 27% in sandy clay to 97% in sand. The variations in leaching losses were attributed to limitations in water content and the effect of HLRs on N transformation rates.

Climate change may produce considerable variations to rainfall rates, sea level and temperatures in the U.S., including the Northeast [31]. These climatic conditions are expected to affect the performance of OWTS and subsequently N transformations in the STA because nitrification and denitrification are sensitive to flow and water content, as well as temperature.

A limited number of studies have investigated the N fate, transport and removal mechanisms of N in STAs receiving ATE [3235]. None of these studies has numerically modeled N transformations in STAs dosed with ATE. Little is known about nitrification and denitrification rates in advanced STAs, and no modeling approach has been developed to simulate these transformation processes. In this study, we addressed this knowledge gap with a calibrated HYDRUS 2D [29] model using soil moisture content and N speciation data collected from a conventional P&S drainfield and two types of shallow narrow drainfields (PSND and Geomat, a type of PSND) using replicated (n = 3) intact soil core mesocosms. We determined nitrification and denitrification rate coefficients for the three drainfield types and used these data to estimate N losses. The simulated results were compared to actual experimental data. The calibrated model was then used to predict the effect of climate changing conditions (increased temperature and rising water table) on nitrification and denitrification rate transformation in the drainfields. The information obtained from these models is expected to aide designers of OWTS and regulators to make informed decisions about the most effective treatment practices for removal of N species in the STA.

Material and Methods

Experimental Setup

Replicated mesocosms (n = 3) were engineered to mimic the soil treatment area and wastewater delivery system of a PSND, Geomat, and P&S. All systems were maintained at 20°C ± 0.7 with a water table separated 90 cm (PSND and Geomat) or 30 cm (P&S) from the infiltrative surface, representative of current climate conditions [20]. Mesocosms consisted of polyvinyl chloride (PVC) pipes (0.15 m ID, 1.5 m H) containing undisturbed soil that is representative of the soil profile used for an STA of an OWTS in southern New England (Morphological, physical and chemical properties of the soil are listed in S1 Table).

Mesocosms were dosed with domestic wastewater, based on accepted guidelines for frequency and volume of wastewater inputs for the State of Rhode Island. For P&S mesocosms, STE was applied at a rate of 400 mL d-1 in two doses of 200 mL over 1.5 h every 12 hours. PSND and GEO mesocosms were dosed with ATE at a rate of 2 L d-1, in 42-mL doses over 15 min every 30 min. The wastewater was dispersed 20 cm below ground surface for PSND, at 25 cm for GEO, and at 84 cm for P&S. The mesocosms were instrumented with probes to collect soil moisture and temperature data.

Effluent samples, along with wastewater inputs, were analyzed weekly for total N, ammonium and nitrate, and other water quality parameters. The physical, chemical and microbiological characteristics of STE and ATE are summarized in S2 Table. Detailed information about soil mesocosms setup and water analysis methodology are found in Cooper et al. [20].

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 θ is the volumetric water content [L3L-3], h is the pressure head [L], S is a sink term [T-1], xi and xj (i and j = 1,2) are the spatial coordinates [L], t is time [T], KijA are components of a dimensionless anisotropy tensor KA, and K is the unsaturated hydraulic conductivity function [LT-1] given by

K(hxyz) = Ks(xyzKr(hxyz)

where Kr is the relative hydraulic conductivity and Ks 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 [36] equation was applied to compute the soil hydraulic properties (Eqs 35):


where α (L-1), m (dimensionless), and n (dimensionless) are fitted parameters, θ(h) is the volumetric water content (L3 L-3), θs is the saturated volumetric water content (L3 L-3), and θr is the residual volumetric water content (L3 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 [37]. 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), Kd is the adsorption coefficient (L3M-1), μ represents the solute transformation or degradation rate in the liquid phase, x is the solute travel distance (L) and z is depth (L). Dwij is the dispersion coefficient tensor for the liquid phase [L2T−1], θ is the volumetric water content [L3L−3], ρ is the bulk density of porous medium [ML−3], and qx and qz is the specific discharge [LT−1] along the horizontal and vertical direction, respectively.

Model Domain and Boundary Conditions

The model domain was developed to resemble the mesocosms, 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 [20]. The model domain consisted of a 2D vertical plane (x-z) (rectangular, L = 15 cm, H = 137 cm high) (Fig 1).

Fig 1
Model domain and porous material distribution for pressurized shallow narrow (PSND), Geomat (GEO), and pipe and stone (P&S) drainfield mesocosms.

The infiltrative surface was placed below the top boundary that shaped the ground surface. The PSND system consists of lateral pipes that distribute the advanced-treated effluent by squirting it against a cover made of larger diameter pipe cut longwise. It was modeled by an arc that represents an impermeable half-pipe cover located above the drainfield. The GEO system is comprised of a core of entangled plastic filaments and a pressure distribution pipe covered with a protective layer of geotextile fabric. GEO was modeled by including a 1-cm filament core layer and a 2.54-cm diameter circle on the top, which simulates the distribution pipe. The P&S model integrates a 30-cm layer (crushed stone or gravel backfill) with an embedded 2.54-cm diameter circle of simulated perforated pipe located 60 cm below the soil surface.

The native soil in the mesocosms is a Bridgehampton silt loam (coarse-silty, mixed, active, mesic Typic Dysturdept) (S1 Table). As mentioned previously, 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 Bw and 2Bw were modeled as one single layer.

A 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. Element 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, a third node located at the bottom of the model domain and a fourth node at the column outlet.

Atmospheric boundary condition was assigned to the top of the columns (Fig 2). The sides and bottom of the column were treated as no-flux boundaries. As wastewater infiltrates, it accumulates at 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 side of each soil column, as indicated in Fig 2. In the HYDRUS program, this occurs when the water is removed by overland flow once saturated conditions prevail [29].

Fig 2
Boundary conditions for (a) pressurized shallow narrow (PSND), (b) Geomat (GEO) and, (c) Pipe and stone (P&S) soil drainfield mesocosms.

N Transformation Modeling

Nitrogen losses in STA are attributed to NH4+ conversion to NO3- or nitrification followed by reduction of NO3- to N2O or N2 through denitrification. Nitrification proceeds in a two-step reaction. First, NH4+ is oxidized to NO2-. Although not specifically investigated, this reaction is assumed to be catalyzed by bacteria through (Eq 7):

NH4++ 1.5O2 NO2+H2O+2H+

Secondly, NO2- is assumed to be oxidized to NO3- by bacteria as follows (Eq 8):

NO2+ 0.5O2 NO3

Denitrification is the transformation of NO3- to its gaseous form (N2) and represents of one of the main mechanism of nitrogen loss in soils. N2 is produced from the reduction of NO3- by heterotrophic-denitrifying bacteria. Denitrification takes place under anaerobic conditions and requires organic carbon (C). NO3- is reduced to N2 through electron transfer from the organic C by the reaction:

C6H12O6+4NO3 6CO2+6H2O+2N2

We developed a decay model to simulate the fate and transport of N species in conventional and advanced STAs in which N was assumed to be transformed as follows [22]:

NH4+  NO3  N2

Nitrite-oxidizing bacteria typically react much faster than NO2-producing bacteria, which restrict the accumulation of NO2- [3841]. For this reason we assumed that NH4+ is converted directly to NO3-, i.e. the production of NO2- (Eq 8) was excluded from the process. From Eq 7, the rate of change of nitrification and denitrification can be expressed as:



where μnit and μdenit are described as nitrification and denitrification reaction rates, respectively. N species were assumed to be dissolved in the STE and ATE and thus the nitrification and denitrification rate coefficients were only determined for the aqueous phase. Nitrogen species were modeled using sequential decay reactions built into HYDRUS [42]. 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 [43,44]. It assumes that the sorption sites are composed of two fractions; sorption to one fraction of sites is assumed to be instantaneous, while for the remaining sites sorption is time-dependent. Also, it assumes that the solute transport in the bulk of the pore space takes place by convection and dispersion. The measured total N (TN) was modeled as an input concentration to include all N infiltrated into the drainfield. The influent organic N was considered to be mineralized to NH4+ through microbial decomposition:


Several researchers have reported that nitrification and denitrification processes are dependent on the water content [45,46]. Nitrification is an aerobic process that occurs at low soil water content because high soil water content increases tortuosity and, limits oxygen diffusion and the activity of nitrifying bacteria [47]. Denitrification takes place under saturated soil conditions, which promotes anoxia. HYDRUS was modified to account for the effect of soil water content on N transformation. For this, a water content dependency function in the DRAINMOD-N2 [48] module allows for computing nitrification and denitrification rates at low water saturation or unsaturated conditions based on Michaelis-Menten kinetics [49]. Both processes are modeled as first-order kinetics when the substrate concentration (NH4+and NO3-) is limited. The reaction kinetics change into zero-order reaction when substrate concentration increase [50]. Because the three drainfield types receive STE and ATE containing high concentrations of NH4+ that are rapidly converted to NO3-, nitrification and denitrification were modeled as zero-order kinetics. For nitrification the model uses a stepwise function to simulate the influence of nitrification inhibitors on decay rates. The following expression describes the nitrification rate:


where μnit is the calculated nitrification rate (ML-1T-1), μnit,max is the maximum nitrification rate (MM-1T-1), CNH4 is the NH4+-N concentration (MM-3), and Km,NH4 is the NH4+-N half-saturation constant (MM-1T-1), which is the NH4+-N concentration at which the nitrification rate is half its maximum value. The parameter fsw is the water content dependency function and is written as follows (Eq 15):


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, Sh and Sl are the upper and lower saturation boundary for optimal nitrification, swp is the saturation level at the wilting point, S is the actual soil saturation or water-filled pore space, and e1 and e2 are fitting parameters.

Denitrification is modeled as a function of the NO3- concentration available and the organic content decrease with depth [41]. The denitrification rate equation included in the modified HYDRUS version is:


where μdenit is the denitrification rate (ML-1T-1), μdenit,max is the maximum denitrification rate (MM-1T-1), CNO3 is the NO3--N concentration, and Km,NO3 is the NO3--N half-saturation constant (MM-1T-1), which is the concentration at which the denitrification rate is half its maximum value. The terms ft and fz are temperature-dependency and carbon dependency functions, respectively.



fsw, dn={0S<Sdn(SSdn1Sdn)fSSdn

ft varies between 0 and 1. T is the temperature, Topt is the optimum temperature for nitrification, which is a system property and typically ranges from 20°C to 25°C, β and α are fitting parameters, and z is depth below the infiltrative surface. The term fsw,dn is the water content-dependency function at a threshold saturation value for denitrification (sdn), s is the actual soil saturation, and f is a fitting exponent.

Climate Change Simulation

Climate projections indicate that precipitation, sea level and temperatures have been increasing in parts of the United States, and this tendency is expected to continue [31]. Increased precipitation and sea level rise may lead to rising water table and a reduced treatment depth, which affect the performance of OWTS, particularly in coastal areas. Warmer soil temperatures are likely to affect N transformations in the drainfield because nitrification and denitrification rates are sensitive to soil temperature and water content, which might restrict the production of NO3- and subsequent transformation to N2 [26]. To predict the response of N transformations, and the fate and transport of N in the three drainfield types to changing climate conditions, the calibrated model was used to simulate nitrification and denitrification under three scenarios: (i) warmer soil temperature (23°C vs. 20°C), (ii) rising water table, where the thickness of the unsaturated zone inside the mesocosms was reduced by 30 cm (from 112 cm to 82 cm), and (iii) simultaneously rising soil temperature and water table elevation. Observation nodes were located at the water table level and the column bottom outlet to record the modeled N species concentrations (only for scenarios ii and iii). The changing climate conditions were set based on climatic projections of eastern U.S. that suggest the water table and temperature will increase about a foot (30 cm) and ±3°C, respectively [31].

Solute Transport Parameters

The transport characteristics were estimated for each of the soil drainfield mesocoms. The longitudinal dispersivity (λL) was set to be one-tenth of the soil profile beneath the infiltrative surface [51,52]. A similar approach was used by others [53,54] to model nitrification and denitrification of septic tank effluent using HYDRUS 2D. Soil bulk density was measured for the silt loam and gravelly-coarse sand (S1 Table). A bulk density value of 1.50 g cm-3 was assigned to the entangled plastic filaments and crushed stone. The diffusion coefficients for NH4+ and NO3- were set to 0.067 and 0.061 cm2 h-1, respectively [53,55].

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 [56]. The nonlinear regression problem was 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 a 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.

The model was calibrated by fitting water content and nitrogen species data (NH4+ and NO3- concentration). HYDRUS water flow and solute transport modules were applied to complete the calibration. First, water content data were fitted to obtain the soil hydraulic parameters and evaluate the impact of moisture content on N transformation. Secondly, NH4+ and NO3- concentration data were used to determine the nitrification and denitrification rates, and to estimate N losses.

The model was initially run under near saturation conditions to reach steady water flow in a shorter simulation time. Therefore, initial average pressure heads were set to -50 cm for the entire model soil profile. Atmospheric boundary conditions were assigned to the top of the model domain or simulated soil surfaces. The minimum permissible pressure was assumed to be -1,000 cm. No precipitation, evapotranspiration or root uptake was included in the simulated N transformation.

Hydraulic loading rates were modeled by assigning a variable flux boundary condition in each of the soil mesocosms. For PSND, it was assumed that wastewater was distributed uniformly over the entire infiltrative surface. For GEO and P&S, the variable flux boundary was located below the distribution pipe. SFE and STE deliveries were modeled as applied in the mesocosm experiments.

Initial values for soil hydraulic parameters were determined by Rosetta [57], a computer program that is part of HYDRUS. The software estimates soil water retention by implementing hierarchical pedotransfer functions (PTFs) based on soil textural classes. Fitting parameter values were assigned to the entangled plastic filaments (GEO) and crushed stone (P&S) systems. Both materials were considered highly conductive (Ks = 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 [57] and fitted with UCODE, whereas values for the plastic filaments and gravel layers were kept fixed. Initial N transformation rates were selected from McCray [22] and initial NH4+ and NO3- soil concentration were set to zero. Water dependency function parameters were selected from McCray et al. [58]. Finally, the model was run for a 3-month (90 days) simulation period. The predicted N species concentrations were computed to estimate the N balance produced by each of the three OWTS. The best fit between the predicted and observed data was evaluated based on the root mean squared error (RMSE). A RMSE value closer to zero indicates the best of fit to observed data.

Results and Discussion

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 found to be under steady-state conditions. The soil hydraulic parameters (θr, θs, α, Ks, n and l) were determined for each of 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 by [29]). PSND and GEO, five soil hydraulic parameters (θr, θs, α, Ks and n) for each of the two horizons (silt loam and gravelly-coarse sand) were calibrated simultaneously (10 parameters total). Five parameters were used for the conventional P&S (Table 1).

Table 1
Model soil hydraulic parameters calibrated in advanced and conventional drainfield mesocosms.

In PSND and GEO, the measured water content (cm-3cm-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 proximity to each other, water content variations were expected at greater depths because of the cumulative 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 heterogeneous behavior of the soil system is illustrated by the three PSND mesocosms, i.e. one mesocosm showed higher water content (0.23 cm-3cm-3) at the 15-cm depth compared to the other two (0.11 to 0.13 cm-3cm-3). In general, these variations are indicative of soils with low residual and high saturated water content characteristics (Fig 3).

Fig 3
Observed and simulated water content for (a) pressurized shallow narrow (PSND), (b) Geomat (GEO) and (c) Pipe and stone (P&S) drainfield mesocosms.

Overall, the model resulted in a good fit between the observed and simulated water content data for PSND, GEO and P&S (Fig 3). RMSE values for PSND and GEO ranged from 0.0010 and 0.0075 cm-3cm-3 for silt loam and gravelly-coarse sand, indicating good agreement between the simulated and measured data. The goodness-of-fit is illustrated in Fig 3, where the model output data were described by a straight line indicating steady conditions during the entire period of simulation at both observation nodes (15-cm and 30-cm depths). In case of the P&S mesocosms, these were dosed with wastewater every 12 h, which resulted in a comparatively drier soil profile and longer times of unsaturated flow between doses relative to the PSND and GEO dosing scheme. Thus, variations in soil moisture content were observed between dosing events, with soil moisture values varying by a factor of two. The water content peaked immediately after dosing (0.03 cm-3cm-3 to 0.05 cm-3cm-3) and dropped quickly (0.01 cm-3cm-3 to 0.02 cm-3cm-3) between doses. Under steady state conditions, the model reproduced those fluctuations with acceptable goodness-of-fit (RMSE: 0.0033 cm-3cm-3 to 0.0044 cm-3cm-3).

The water content data were modeled under the effect of a simulated hanging water table at the bottom of the mesocosms, where the seepage face boundary (Fig 4) caused this part of the model domain to remain saturated once the system was at steady state. The calibrated retention curve parameters are shown in Table 2.

Fig 4
Pressure head distribution as a result of the seepage boundary condition to simulate a hanging water table at the bottom of the mesocosms.
Table 2
Calibrated soil hydraulic parameters for the simulation of pressurized shallow narrow (PSND), Geomat (GEO) and pipe and stone (P&S) drainfield mesocosms.

The calibrated values differed among soil layers, which indicate that the properties of the soil at the infiltrative surface were different from the underlying soil, likely due to differences in soil texture and structure. Based on the soil moisture data, the silt loam was less conductive and had higher saturated water contents. The underlying soil (gravelly-coarse sand) for the PSND and GEO was simulated with Ksg ranging from 908.9 to 942.5 cm day-1, which were 21% to 44% higher than reported values for sandy soil [59,60]. Although variations in hydraulic conductivity values are expected, in our study the higher Ksg were likely caused by the presence of significant amounts of gravel, which accounted for 40% to 45% of the soil by weight. These differences in physical properties affect the soil properties directly, and influence the hydraulic properties and water flow. The gravelly-coarse sand layer for the P&S drainfield mesocosms was simulated with an average hydraulic conductivity of 4.51 cm day-1 (Table 2). As indicated in Fig 4, it is most likely that a biomat developed overtime above the infiltrative surface, which results in unsaturated conditions and a reduced saturated hydraulic conductivity.

Sensitivity Analysis

The sensitivity analysis was focused on five soil hydraulic parameters (θr, θs, α, Ks and n). The 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 are shown in Table 3. Most of the selected parameters were significantly sensitive (CSS ≥ 0.01) to the water content data. For the advanced STAs, the simulated soil moisture was sensitive to the silt loam soil properties. Not unexpectedly, the soil properties θss and ns were most important for the calibration of the hydraulic parameters along the soil profile. Generally, Kss, θrs, θrg and Kg were either not significant or less sensitive parameters. For P&S, the saturated and residual water content (θrg and θsg) were very important parameters determining the soil moisture distribution along the profile. In addition, the hydraulic conductivity (Ksg) was more sensitive compared to PSND and GEO (CSS = 0.21 to 0.25).

Table 3
Composite scale sensitivity ratios to the measured soil moisture data for the silt loam and gravelly-coarse sand soils for pressurized shallow narrow (PSND), Geomat (GEO) and pipe and stone (P&S) drainfield mesocosms.

In one of the PSND columns (Table 3, column #3), the Ksg was found to be less sensitive to fitting water content data (CSS < 0.01). In this mesocosm, the water content of the silt loam was almost two times higher (0.23 cm-3cm-3) than in the other two PSND columns (0.11 cm-3cm-3 and 0.13 cm-3cm-3). These variations are likely linked to soil heterogeneity.

Nitrogen Transport and Fate

Nitrification and denitrification were modeled using a water content-dependent function to account for changes in oxygen diffusion and its availability in the mesocosms. The function uses water-filled pore space to mimic soil aeration during water infiltration [61]. Based on this approach, NO3- 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 N2 gas production rate is observed at saturation (WFPS = 1.00) [62,63]. Linn and Doran [63] reported that organic carbon decomposition associated with N mineralization and immobilization occurs when WFPS ranges from 0.50 to 0.60 as well as saturation. 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 [62,63] had a higher dissolved oxygen concentration compared to the STE and ATE used in this study. Based on 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 than previously described.

The nitrification and denitrification rate coefficients were computed using Eqs 14 through 19, and parameter values were selected from literature data [58]. Also, the parameter values for the water-content dependent functions were fitted. Initially, the model was adjusted until the best fit was achieved between the observed and predicted data. As a result, the parameters for nitrification and denitrification dependency functions are median values that best reproduced the observed data (fwp = 0, fs = 0, swp = 0.154, sl = 0.665, sh = 0.809, e1 = 2.267, e2 = 1.104, sdn = 0 and f = 2.86) [58].

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 (Table 4).

Table 4
Modeled water-filled pore space (WFPS) for pressurized shallow narrow drainfield (PSND), Geomat (GEO) and pipe and stone (P&S) drainfield mesocosms.

This indicates that sufficient oxygen is available for nitrification to proceed. Compared to the gravelly-coarse sand, the silt loam material has the highest values for the modeled WFPS in both PSND and Geo (0.64 and 0.74, respectively). A similar value (0.76) was reported by Bradshaw et al. [26] when simulating nitrification and denitrification rates from an OWTS drainfield installed in a clay-textured soil using pressure head and NH4+and NO3- concentration data to simulate the system. Their model converted the pressure heads into water content values to calculate the actual WFPS of the drainfield. It also captured the effect of seasonal changes (dry and wet weather) on N transformation and found that the computed WFPS was adequate for nitrification to occur.

Our results are consistent with what is expected for the soil types and the hydraulic properties of our mesocosm materials. The data indicate that nitrification occurred in the first few centimeters below the infiltrative surface. Nitrate production in all drainfields and at shallow depths (top 15 cm) is likely caused by the oxidation of ammonia by ammonia-oxidizing and nitrifying bacteria [20].

The predicted and measured NH4+ concentrations for all drainfield types are shown in Fig 5. Except for the pipe and stone (P&S) drainfield mesocosm (Fig 5c), the model output show a good fit relative to the measured NH4+ concentrations in effluent water, with RMSE values ranging between 0.18 and 0.45 mg L-1. The P&S model resulted in higher RMSE values (2.18 to 2.88 mg L-1) because of initially elevated NH4+ concentrations. The cause for the elevated NH4+ concentrations is unknown and including these data caused the model to exceed the measured ammonia concentrations.

Fig 5
Predicted and measured NH4+ concentrations for (a) pressurized shallow narrow (PSND), (b) GEOMAT (GEO) and (c) pipe and stone (P&S) drainfield mesocosms.

The maximum NH4+ concentration was found near the infiltrative surface (top 15 cm) and decreased with depth along the soil profile. The model results showed that the NH4+ was almost completely transformed at the 30-cm depth, consistent with other studies [25,64,65]. Moreover, the lowest measured and modeled NH4+ concentrations were observed in the outflow, where almost no NH4+ was detected.

The NO3- concentration in ATE inputs and water exiting the mesocoms were measured and the data used to calibrate the models for the advanced and conventional STAs (Fig 6). In PSND and GEO, influent total N included NO3- and NH4+. Some of the nitrate resulted from NH4+ being nitrified in the sand filter that preceded the treatment system from which the ATE was collected. Nitrate tended to increase with depth along the soil profile in all mesocosms, with the highest concentration detected near the outlet (seepage boundary).

Fig 6
Predicted and measured NO3- concentrations for (a) pressurized shallow narrow (PSND), (b) Geomat (GEO) and (c) pipe and stone (P&S) drainifield mesocosms.

For ATE, the model output included NO3- initially present in the influent water as well as the NO3- produced in situ from NH4+ oxidation. The model suggests that the remaining NH4+ will be transformed to NO3- in the drainfield.

The predicted NO3- concentrations showed an acceptable goodness-of-fit with the observed data, with RMSEs ranging from 4.45 mg L-1 d-1 to 9.65 mg L-1 d-1 in all STA types (Fig 6). Lower RMSE values were observed for PSND and GEO compared to P&S. This observation was not unexpected because ATE likely is more uniformly distributed over the infiltrative surface in the PSND and GEO.

Nitrification and Denitrification Rates

Nitrogen transformation and removal is mainly controlled by nitrification and denitrification. In addition, NH4+ sorption to soil can affect the fate and transport of N in some OWTS drainfields. Because of the low sorption capacity of the soils studied herein (S1 Table), NH4+ sorption was not considered. Therefore, all NH4+ must move with the 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 5).

Table 5
Average zero-order nitrification and denitrification rates for the selected soils and materials in pressurized shallow narrow (PSND), Geomat (GEO) and pipe and stone (P&S) drainfield mesocosms.

The nitrification rates ranged from 0.5 mg L-1 d-1 to 574 mg L-1 d-1 and were similar to zero-order rate values reported by McCray et al. [22]. Geza et al. [66] developed a tool for predicting the fate and transport of nitrogen in STAs (STUMOD), which uses nitrification rates as an input parameter. Their default value is 56 mg N L-1 d-1, which is similar to nitrification rates modeled herein. Overall, the advanced OWTS drainfields showed higher nitrification rates compared to P&S. As summarized in Table 5, the average PSND zero-order nitrification rates for silt loam and gravelly coarse sand were 45.25 mg N L-1 d-1 and 49.19 mg N L-1 d-1, respectively. Lower values were computed for GEO (2.17 mg N L-1 d-1 and 25.88 mg N L-1 d-1 for silt loam and gravelly-coarse sand, respectively). The model results suggest that some nitrification occurred in the entangled plastic filaments (25.88 mg N L-1 d-1). Nitrate production at that interface may be attributed to high oxygen diffusion and ATE aeration during the infiltration process. The average nitrification rates were 3.83 mg N L-1 d-1 in the gravelly coarse sand for the P&S. Nitrification took place at a rate of 12.10 mg N L-1 d-1 in the crushed stone and was 0.5 times lower than that computed for the GEO plastic filaments (25.88 mg N L-1 d-1). These data indicate that the presence of a conductive layer on the top of the infiltrating native soils provides an additional treatment zone for N removal. Furthermore, the higher NH4+ transformation rates in the advanced OWTS suggest that the drainfield placement at a shallower depth is more effective for nitrification than the conventional systems, likely because of a larger volume of unsaturated soil available for treatment.

Denitrification was not very significant in any of the OWTS drainfields. Denitrification rate values were one to three orders of magnitude lower than nitrification rates (from 0.01 to 0.44 mg N L-1 d-1). Tucholke et al. [67] reported higher zero-order denitrification rates values, ranging between 0.033 and 127 mg N L-1 d-1. However, those values [67] were obtained under fully saturated conditions (WFPS = 100%). Because unsaturated conditions prevailed in all mesocosms discussed herein, denitrification may have been restricted, since it requires anaerobic conditions to proceed [53]. Anaerobic conditions are more likely under saturated flow conditions.

Denitrification rates were higher in P&S than GEO and PSND. This finding is consistent with the experimental results presented in [20], 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 low organic carbon content, it may have further limited the extent of denitrification in the advanced drainfield mesocosms. This would be consistent with [20].

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 NO3- (82.72%) and NH4+ (1.41%). In GEO and PSND, the modeled effluent N speciation consisted of 89–91% NO3- and 0.23–0.44% NH4+. The model results indicate that the total N losses as N2 were 10.44%, 9.65%, 17.60% for PSND, Geo and P&S, respectively. There were discrepancies between the computed and observed NO3- data, particularly for the N removal in P&S. That is, the computed NO3- data underestimated the measured data in some instances. It is likely that not all organic N was converted to NH4+ and as a result, less NH4+ was nitrified, whereas for our modeling approach, we assumed that organic N is completely transformed to NH4+ before entering the treatment system. Organic N accounted for 14% to 16% [20] 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 and passing through the drainfield unchanged. For GEO and PSND, the modeled N losses occurred mostly as NO3- (90.75% and 88.45%, respectively). Only minor percentages of NH4+ were observed during the 90-days simulation period (0.23 to 1.41% for all drainfield types). Nitrogen losses as N2 were more evident in P&S compared to the advanced technologies.

Nitrification and Denitrification under Changing Climate Conditions

Our model predicted the effect of warmer temperature and water table located at a shallower distance from the infiltrative surface on nitrification and denitrification in the three different drainfield type. The modeled NH4+ output concentrations were higher at the water table level in all simulated scenarios (Figs (Figs779), indicating that nitrification was affected by the simulated environmental stresses, which resulted in less NH4+ converted to NO3-.

Fig 7
Predicted output NH4+, NO3-, and N2 concentration for pressurized shallow narrow drainfield (PSND) mesocosm under changing climate conditions.
Fig 9
Predicted output NH4+, NO3-, and N2 concentration for pipe and stone (P&S) drainfield mesocosm under climate changing conditions.
Fig 8
Predicted output NH4+, NO3-, and N2 concentrations for Geomat (GEO) drainfield mesocosm under climate changing conditions.

Ammonium output concentrations were more pronounced when the water table was raised and the soil temperature was increased, simultaneously. For example, the average NH4+ concentration was increased about one order of magnitude in the GEO (from 0.21 mg L-1 to 3.71 mg L-1) compared to those simulated at current conditions (soil temperature: 20°C and no rising water table). For PSND, when the water table was located at a shallower distance from the infiltrative surface and warmer soil temperature, only small differences were observed on average output NH4+ concentration, which was increased from 0.11 mg L-1 to 0.83 mg L-1. The greatest effect of climate changing conditions on nitrification was observed in P&S mesocosms, where the average effluent NH4+ concentration was about six times higher than under current conditions. Since N transformation was assumed to be the result of a sequential decay, NO3- concentrations were also influenced by changing climate conditions. As a result, NO3- concentrations decreased due to the effect of rising water table and warmer temperature. The relatively lower nitrification rate at the water table level can be attributed to a reduction in the amount of unsaturated soil, particularly the vertical separation between the infiltrative surface and the water table.

Denitrification was also restricted when our model for all drainfield types was run under climate changing conditions. Variations in N2 output concentration were observed in response to the differences in temperature and water table level. For instance, denitrification was most affected when the model was run with the water table located at shallower distance from the infiltrative surface and warmer soil temperature. Similar results were obtained when both conditions were modeled separately. In most mesocosms, N2 production increased with increasing temperature. At 23°C and current separation distance, N losses as N2 increased 22%, 37% and 21% for GEO, PSND and P&S, respectively. As the water table was located closer to the infiltrative surface, denitrification accounted for 100% of input NH4+ concentration at the saturation zone (column bottom) and NO3- was not detected below the water table. This is attributed to the complete conversion of NO3- to N2 in the saturated zone. The model incorporates a carbon dependency function to account for the effect of the spatial variability of C on denitrification rate (Eq 18) [61]. In this context, C availability is simulated as a function of soil depth (C content is decreased with depth in the soil profile). Therefore, denitrification was not expected below the water table due to carbon limitation. We assumed that the NO3- plume would be likely transported and diffused rather than converted to N2 through the saturated soil.


We developed a model to predict the fate, transport and transformation of nitrogen species in a conventional P&S drainfield and in two types of shallow narrow drainfields (PSND and GEO) under current and climate changing conditions. The model was used to gain knowledge about the effect of soil water content on nitrification and denitrification rates in soil drainfields. Based on the calibrated model, UCODE was an efficient tool to determine water flow and transport parameters with acceptable goodness-of-fit between the observed and simulated data. However, when soil water content data are fitted, additional parameter information (i.e., soil hydraulic conductivity, Ks) should be gathered to increase the model precision and to provide more realistic results.

The water content dependency functions were suitable to compute nitrification and denitrification rates in all drainfield types. Nevertheless, more research has to be done to determine the relationship between WFPS and relative rate of microbial nitrification and denitrification for STE and ATE. Additionally, this information has to be incorporated in the model to obtain better results.

When the water table is located at shallower distance from the infiltrative surface and under warmer soil temperature conditions, the model predicted that nitrification and denitrification rates are likely to be affected in all drainfield types. Nitrification was restricted (NO3- concentrations were reduced) and subsequently N losses as N2 were increased between 21% and 37%. The simulation of the fate and transport of N in STAs suggest that climate change may positively influence the performance of OWTS, and also reveals that these treatment systems will response to changing temperature and rising water table conditions, predicted by climate assessment models. These results allow for the quantification the N losses in all OWTS drainfield types and an estimation of the N species fluxes under changing climate conditions. This information is useful to better understand the N transport and transformation mechanisms in OWTS now and in the future.

Our results are more representative for the temperate climates of the U.S. Northeast or Mid-Atlantic region. Both higher and lower initial soil temperatures and corresponding climate changes amplitudes are possible in other parts of the world. Therefore, additional modeling studies, together with field investigations, have to be conducted to determine if other temperature amplitudes than simulated herein would result in similar effects on contaminant removal.

Supporting Information

S1 Table

Select morphological, physical and chemical properties of the soil used in drainfield mesocosms.

Values for physical and chemical properties are means (n = 7) ± s.d. Measurements of pH, electrical conductance (EC) and cation exchange capacity (CEC) were made on composite samples. Cooper et al. (2015).


S2 Table

Characteristics of septic tank effluent (STE) and sand filter effluent (SFE) used in our study (n = 26–49).

Cooper et al. (2015).



We thank Dr. Mengistu Geza at Colorado School of Mines for his technical support.

Funding Statement

Funded by U.S. Department of Agriculture - Multistate Project (NE 1045).

Data Availability

Data Availability

All output data files are also available from the Figshare database (, DOI: For additional data, please contact the authors at ude.iru@gnivob.


1. Giblin AE, Gaines AG. Nitrogen inputs to a marine embayment: the importance of groundwater. Biogeochemistry. 1999. June 1;10(3):309–28.
2. Harden HS, Roeder E, Hooks M, Chanton JP. Evaluation of onsite sewage treatment and disposal systems in shallow karst terrain. Water Res. 2008. May;42(10–11):2585–97. doi: 10.1016/j.watres.2008.01.008 [PubMed]
3. Swartz CH, Reddy S, Benotti MJ, Yin H, Barber LB, Brownawell BJ, et al. Steroid Estrogens, Nonylphenol Ethoxylate Metabolites, and Other Wastewater Contaminants in Groundwater Affected by a Residential Septic System on Cape Cod, MA. Environ Sci Technol. American Chemical Society; 2006. August 18;40(16):4894–902. [PubMed]
4. Standley L., Swartz C, Rudel RA, Attfield KR, C J., Erickson M, et al. Wastewater-contaminated Groundwater as a Source Of Endogenous Hormones and Pharmaceuticals to Surface Water Ecosystems. 2008. [PubMed]
5. Tyler E.J., Laak R., McCoy E., Sandhu SS. The Soil as a treatment system. In: Proceedings of the Second National Home Sewage Treatment Symposium. St. Joseph, MI: ASAE; 1977. p. 5–77.
6. Hargett DL. Performance assessment of low pressure pipe wastewater injection systems. In: 4th National Symposium on Individual and Small Community Sewage Systems. St. Thomas, MI; 1985. p. 131–43.
7. Stewart L.W. and Reneau RB. Shallowly placed, low pressure distribution system to treat domestic wastewater in soil with fluctuating high water tables. J Environ Qual. 1988;(17):499–504.
8. Reneau RB, Hagedorn C, Degen MJ. Fate and Transport of Biological and Inorganic Contaminants from On-Site Disposal of Domestic Wastewater. J Environ Qual. 1989;135–44.
9. Anderson D.L., Otis R.J., McNeillie J.I., Apfel RA. In-situ lysimeter investigation of pollutant attenuation in the vadose zone of a fine sand. In: Proceedings of the Seventh International Symposium on Individual and Smaill Community Sewage Systems. St. Joseph, MI: ASAE; 1994. p. 209–18.
10. Powelson DK, Gerba CP. Virus removal from sewage effluents during saturated and unsaturated flow through soil columns. Water Res. 1994;28(10):2175–81.
11. Stevik TK, Aa K, Ausland G, Hanssen JF. Retention and removal of pathogenic bacteria in wastewater percolating through porous media: A review. Water Research. 2004. p. 1355–67. [PubMed]
12. Birkham TK, Hendry MJ, Wassenaar LI, Mendoza CA. A transient model of vadose zone reaction rates using oxygen isotopes and carbon dioxide. Vadose Zo J. 2007;6:67–76.
13. Carlile BL. Use of shallow, low pressure injection systems in large and small installations. In: In NI McClelland (ed) Proceedings of the 6th National Conference on Individual Onsite Wastewater Systems. Ann Harbor, MI: National Sanitation Foundation; 1980. p. 371–85.
14. Environmental Protection Agency U.S.. Watershed Assessment, Tracking & Environmental Results [Internet]. 2015 [cited 2003 May 20]. Available: Accessed 20 May 2003.
15. US Environmental Protection Agency; Wadeable Streams Assessment: A Collaborative Survey of the Nation’s Streams. Office of Water, US Environmental Protection Agency, Washington DC: 2006.
16. Loomis GW, Dow DB, Stolt MH, Green LT, Gold AJ. Evaluation of innovative onsite wastewater treatment systems in the green hill pond watershed, Rhode Island—A NODP II project update. On-Site Wastewater Treat Proc. 2001;506–15.
17. Brady NC, Weil RR. The nature and properties of soils. Macmillan Publishing Co; New York: 2002. 960 p.
18. Burks D.B. and Minnis MM. Onsite Wastewater Treatment Systems. Madison, WI: Hogart House, Ltd; 1994. 248 p.
19. Siegrist R.L., Tyler E.J., Jenssen PD. Design and Performance of Onsite Wastewater Soil Absorption Systems. National Research Needs Conference: Risk-Based Decision Making for Onsite Wastewater Treatment. Palo Alto, CA: EPRI; Report; 2001.
20. Cooper JA, Loomis GW, Kalen D V, Amador JA. Evaluation of Water Quality Functions of Conventional and Advanced Soil-Based Onsite Wastewater Treatment Systems. J Environ Qual. 2015. February. [PubMed]
21. Beal CD, Gardner EA, Menzies NW. Process, performance, and pollution potential: A review of septic tank-soil absorption systems. Australian Journal of Soil Research. 2005. p. 781–802.
22. McCray JE, Kirkland SL, Siegrist RL, Thyne GD. Model parameters for simulating fate and transport of on-site wastewater nutrients. Ground Water. 2005. p. 628–39. [PubMed]
23. Valiela I, Bowen JL, Kroeger KD. Assessment of models for estimation of land-derived nitrogen loads to shallow estuaries. Appl Geochemistry. 2002;17(7):935–53.
24. Hassan G, Reneau RB, Hagedorn C, Jantrania AR. Modeling Effluent Distribution and Nitrate Transport through an On-Site Wastewater System. 2008. September;1937–48. [PubMed]
25. Heatwole KK, McCray JE. Modeling potential vadose-zone transport of nitrogen from onsite wastewater systems at the development scale. J Contam Hydrol. 2007;91(1–2):184–201. [PubMed]
26. Bradshaw JK, Radcliffe DE, Šimůnek J, Wunsch A, McCray JE. Nitrogen Fate and Transport in a Conventional Onsite Wastewater Treatment System Installed in a Clay Soil: A Nitrogen Chain Model. Vadose Zo J. 2013. August;12(3):1–20.
27. Radcliffe DE, Bradshaw JK. Model Test of Proposed Loading Rates for Onsite Wastewater Treatment Systems. 2014;(6):97–107.
28. Morales I, Atoyan JA, Amador JA, Boving T. Transport of pathogen surrogates in soil treatment units: Numerical modeling. Water (Switzerland). 2014;6(4):818–38.
29. Simunek J, Genuchten M Van, Sejna M. The HYDRUS software package for simulating the two-and three-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Technical manual. 2012.
30. Šejna M, van Genuchten MT, Šimůnek J. Development and Applications of the HYDRUS and STANMOD Software Packages and Related Codes. Vadose Zone Journal. 2008. p. 587.
31. Kirtman B, Power S, Adedoyin J, Boer G. Near-term climate change: projections and predictability. Clim Chang. 2013;953.
32. Holden S.A., Stolt M.H., Loomis G.W., Gold AJ. Seasonal variation in nitrogen leaching from shallow narrow drainfields. In: Mankin, editor. Proceedings of the Tenth International Symposium on Individual and Small Community Sewage Systems. St. Joseph, MI: ASAE; 2004. p. 432–40.
33. Holden SA. The effectiveness of shallow-narrow drainfields to treat domestic wastewater. University of Rhode Island; 2004.
34. Gill LW O’Súlleabháin C, Misstear BDR, Johnston PJ. The Treatment Performance of Different Subsoils in Ireland Receiving On-Site Wastewater Effluent. J Environ Qual. 2007;(6):1843–55. [PubMed]
35. Gill L O’Luanaigh N, Johnston P, Misstear BDR, O’Suilleabhain C. Nutrient loading on subsoils from on-site wastewater effluent, comparing septic tank and secondary treatment systems. Water Res. 2009;43(10):2739–49. doi: 10.1016/j.watres.2009.03.024 [PubMed]
36. van Genuchten MT. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils1. Soil Sci Soc Am J. 1980;44:892–8.
37. Mualem Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour Res. 1976;12(3):513–22.
38. Tyler K, Broadbent F. Nitrite transformations in California soils. Soil Sci Soc Am …. 1960.
39. Curds C, Hawkes H. Ecological aspects of used-water treatment, Vol. 2 Biological activities and treatment processes. 1983.
40. Tchobanoglous G, Burton FL, Stensel HD. Wastewater Engineering: Treatment and Reuse. Metcalf & Eddy, Inc. Engineering; 2003. 1819 p.
41. Grady CPL Jr., Daigger GT, Love NG, Filipe CDM. Biological Wastewater Treatment, Third Edition CRC Press; 2011. 1022 p.
42. Šimůnek J. van Genuchten MT. Numerical model for simulating multiple solute transport in variably-saturated soils In: Wrobel L., Latinopoulos P, editor. Water Pollution III: Modelling, Measurement, and Prediction. Southampton, UK: Computation Mechanics Publication, Ashurst Lodge, Ashurst; 1995. p. 21–30.
43. Selim HM, Schulin R, Flühler H. Transport and Ion Exchange of Calcium and Magnesium in an Aggregated Soil. Soil Sci Soc Am J. 1987;51(4):876.
44. van Genuchten MT, Wagenet RJ. Two-Site/Two-Region Models for Pesticide Transport and Degradation: Theoretical Development and Analytical Solutions. Soil Sci Soc Am J. 1989;53(5):1303.
45. Grant RF. Mathematical modelling of nitrous oxide evolution during nitrification. Soil Biol Biochem. 1995;27(9):1117–25.
46. Bollmann A, Conrad R. Acetylene blockage technique leads to underestimation of denitrificat1on rates in oxic soils due to scavenging of intermediate nitric oxide. Soil Biol Biochem. 1997;29(7):1067–77.
47. Millington RJ, Quirk JP. Transport in porous media. Trans 7th int Congr Soil Sci. 1960;1:97–106.
48. Youssef MA, Youssef MA, Skaggs RW, Skaggs RW, Chescheir GM, Chescheir GM, et al. Drainmod − n ii. Society. 2005;48(2):611–26.
49. Malhi SS, McGill WB. Nitrification in three Alberta soils: effects of temperature, moisture and substrate concentration. Soil Biol Biochem. 1982;14(4):393–9.
50. Prosser JI. Mathematical Modeling Of Nitrification Processes In: Marshall KC, editor. Advances in Microbial Ecology SE—7. Springer; US; 1990. p. 263–304.
51. Gelhar LW, Welty C, Rehfeldt KR. A critical review of data on field-scale dispersion in aquifers. Water Resour Res. 1992. July 9;28(7):1955–74.
52. Vanderborght J, Vereecken H. Review of Dispersivities for Transport Modeling in Soils. Vadose Zo J. 2007. February 1;6(1):29.
53. Cote CM, Bristow KL, Charlesworth PB, Cook FJ, Thorburn PJ. Analysis of soil wetting and solute transport in subsurface trickle irrigation. Irrig Sci. 2003. November 1;22(3–4):143–56.
54. Beggs R.A., Tchobanoglous G., Hills D., Crites RW. Modeling subsurface drip application of on-site wastewater treatment system effluent. In: tenth National Symposium on individual and small community sewage treatment. Sacramento, CA; 2004.
55. Weast R. CRC handbook of chemistry and physics: a ready-reference book of chemical and physical data. Boca Raton Fla: CRC Press; 1985.
56. Poeter EP, Hill MC. UCODE, a computer code for universal inverse modeling. Comput Geosci. 1999;25(4):457–62.
57. Schaap MG, Leij FJ, Van Genuchten MT. Rosetta: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J Hydrol. 2001;251(3–4):163–76.
58. McCray J., Geza M., Lowe K., Boving T., Radcliffe D., Tucholke M., Wunsch A., Roberts S., Amador J., Atoyan J., Drewes J., Kalen GL D.. Quantitative Tools to Determine the Expected Performance of Wastewater Soil Treatment Units: Guidance Manual. 2010.
59. Brooks R. H., Corey AT. Properties of porous media affecting fluid flow. J Irrig Drain Div. 1966;2(1992):61–90.
60. van Genuchten MT. A Closed-form Equation for Prediccting Hydraulic Conductivity of Unsaturated Soils. Soil Sci Soc Am J. 1980;44(5):892–8.
61. Brevé MA. Modeling the movement and fate of nitrogen in artificially drained soils. Ph.D. diss., North Carolina State University, Raleigh, N.C.; 1994.
62. Doran J. W., Mielke L.N., Stamatiadis S. Microbial activity and N cycling as regulated by soil water-filled pore space. In: Tillage and traffic in crop production Proc 11th Int Soil Tillage Res Organization. 1988. p. 49–54.
63. Linn DM, Doran JW. Effect of Water-Filled Pore Space on Carbon Dioxide and Nitrous Oxide Production in Tilled and Nontilled Soils1. Soil Science Society of America Journal. 1984. p. 1267.
64. Fischer EA. Nutrient Transformation and Fate During Intermittent Sand Filtration of Wastewater. MS Thesis, Colorado School of Mines, Golden, CO.; 1999.
65. Beach DNH. Infiltration of wastewater in sand columns. MS Thesis, Colorado School of Mines, Golden, CO; 2001.
66. Geza M, Lowe K, McCray J. STUMOD—a Tool for Predicting Fate and Transport of Nitrogen in Soil Treatment Units. Environ Model Assess. Springer International Publishing; 2014;19(3):243–56.
67. Tucholke M. B., McCray J. E., Thyne G. D., Waskom RM. Variability in denitrification rates: literature review and analysis. In: National Onsite Wastewater Recycling Association 16th Annual Technical Education & Exposition Conference. Baltimore, MD: NOWRA, Santa Cruz, CA; 2007.

Articles from PLoS ONE are provided here courtesy of Public Library of Science