On-site wastewater systems (OWSs) are considered the cause of significant non-point source pollution to water bodies, especially in rural areas and suburban areas. Quantifying their effects on water quality is important (McCray et al., 2005). Among the several conceptual and (or) mathematical models, a biozone algorithm proposed by Siegrist et al. (2005) was adapted to current work due to the fact that the algorithm fully describes biozone processes, validated at watershed scale, and provides coefficients specifically developed for biozone processes that were calibrated to field scale experimental data. In this algorithm, septic tank effluent directly drains into subsurface soil layer (infiltrative surfaces) affecting soil moisture content and the percolation of soil water through the vadose zone. The amount of percolation is regulated by unsaturated hydraulic conductivity as a function of soil moisture content. Total suspended solid (TSS), biochemical oxygen demand (BOD) and nutrients in the septic tank effluent (STE) trigger complex bio-physical processes in the biozone which, in turn, affect the hydraulics of soil water. The accumulation of TSS and plaque (i.e., dead bodies of biomass) in the pore space causes clogging and eventually hydraulic failure of the onsite wastewater system. The concentration of nutrients, BOD, and Fecal Coliform is estimated by a first-order decay equation based on the reaction coefficients specifically developed for each constituent by Siegrist et al. (2005). The biozone is considered as a control volume either on top of soil layer 1 or wedged in between soil layer 1 and soil layer 2 depending on the soil thickness. This control volume will accept the OWS effluent as well as infiltration from above. It can output flow vertically below and can contribute to detention storage on the surface (if the biozone layer is fully saturated).

Soil clogging decreases porosity in the biozone with inert and biological materials. With the reduced soil porosity, the hydraulic conductivity of soil decreases with time (USEPA, 1980). A field-scale experiment suggests the reduction in hydraulic conductivity in the biozone is primarily influenced by STE loading rate and the type of infiltrative surface (Bumgarner and McCray, 2007). Weintraub et al. (2002) proposed a relationship between the biozone hydraulic conductivity and soil moisture contents.

$K_{bz}=K_{sat}\frac{\theta -\theta_f}{\theta_s-\theta_f}$ (7)

where $K_{bz}$ is biozone hydraulic conductivity (mm/hr), $K_{sat}$ is saturated hydraulic conductivity of soil (mm/hr), and $\theta$ is moisture content of biozone (mm). An advantage of this model is that, $K_b$ is directly related to $K_s$ and the soil moisture content. The theoretical basis on the formulation of the equation is not presented in the literature; instead, the formula is indirectly validated by comparing percolation to the subsoil layer as a function of time.

Unlike natural soil conditions, the field capacity of a biozone dynamically changes with time due to the development of filamentous material of live bacterial biomass allowing the biozone layer to retain additional water. Therefore, temporal change in the biozone field capacity is related to the amount of biomass in the layer. This is shown by the equation:

$\theta_f^t=\theta_f^{t-1}+\Phi*(\theta_s^{t-1}-\theta_f^{t-1})^{\xi}\frac{[\frac{d(Bio)}{dt}]}{10*\rho_{bm}}$ (6)

where $\theta_f^t$ is field capacity at the end of the day (mm), $\theta_f^{t-1}$ is field capacity at the beginning of the day (mm), $\theta_s^{t-1}$ is saturated moisture content at the beginning of the day (mm), $\rho_{bm}$ is the density of live bacterial biomass (~1000 kg/m3$\Box \Box \Box$), Φ is field capacity coefficient 1 (unitless), and $\xi$ is field capacity coefficient 2 (unitless).

The biozone layer is formed as a biologically active layer in the soil absorption system near the infiltrative surface by the growth of microorganisms feeding on the organic matter (BOD) of the septic tank effluent (See Figure 6:4-1). The amount of live bacteria biomass in the biozone is estimated using a mass balance equation assuming the biozone as a control volume. The mass balance equation of microorganisms (live bacteria biomass) in the control volume is then estimated by:

(1)

where is the amount of live bacteria biomass in biozone (kg/ha), ,in is the concentration in the (mg/L), is the concentration in biozone (mg/L), is gram of live bacteria growth to gram of in STE (conversion factor), is the flow rate of (m/day), Ip is the amount of percolation out of the biozone (m/day), is the amount of respiration of bacteria (kg/ha), is the amount of mortality of bacteria (kg/ha), and is the amount of sloughed off bacteria (kg/ha).

Fate and transport of biomass, including respiration, mortality, and slough-off, are estimated based on empirical relationships. Respiration and mortality rates are functions of the amount of the biomass. These values are normalized by the unit area (1/ha) so that these equations are applicable in different scales of simulations without unintended amplification due to a higher mass of biomass. For each time step, a portion of live biomass is removed during respiration and death. The reaction for bacterial respiration is calculated as follows.

$R_{resp}=\gamma *Bio$ (2)

where $\gamma$ is a respiration rate coefficient (unitless). The reaction for bacterial mortality is calculated by

$R_{mort}=\phi*Bio$ (3)

where $\phi$ is mortality rate coefficient (unitless). Bacterial biomass can be washed off to the subsoil layer by a high velocity of infiltrating water.

$R_{slough}=\eta*v_p^\delta$ (4)

where $v_p$ is the pore velocity in the biozone layer (mm/day), $\eta$ is a linear coefficient (kg/ha), and $\delta$ is an exponential coefficient (unitless). Equations (2) to (4) are highly dependent on empirical calibration coefficients and the nature of these processes makes it difficult to validate the model equations.

A portion of dead body of biomass becomes plaque. Total solids in the STE may contribute to increasing plaque accumulation in the pore space. Plaque can be sloughed off from biozone by high pore velocity of infiltrating water. As the amount of live biomass increases in the biozone, plaque also increases. The rate of change in plaque is computed by a mass balance equation.

$\frac{d(plaque)}{dt}=R_{mort}+\frac{\sigma *\sum Q_{STE}*TS}{1000*A_d}-R_{slough}$ (5)

where $plaque$ is the amount of dead bacteria biomass and residue (kg/ha), $\sigma$ is a calibration parameter to convert total solids in $STE$ to $plaque$ (unitless), $TS$ is the total solids contained in STE (mg/l), and $A_d$ is the area of drain field (ha).

Phosphorus adsorption takes place in the soil media below the biozone. The concentration of P in the biozone is often in the linear range of reported nonlinear isotherms (McCray et al., 2005). A linear isotherm is represented by the equation:

$S=K_DC$ (15)

where $S$ is the mass of solute sorbed per unit dry weight of solid (mg/kg), $C$ is the concentration of the solute in solution in equilibrium with the mass of solute sorbed onto the solid (mg/L), and $K_D$ is a linear distribution coefficient (L/kg). McCray et al. (2005) recommends $K_D$= 15.1 L/kg, the linear sorption isotherm constant as median value, but the value may vary from the 10th percentile ($K_D$= 5 L/kg) to 90th percentile ($K_D$= 128 L/kg) for modeling purpose. Similarly, a median value of $S_{max}$= 237 mg/kg is recommended for the maximum $P$ sorption capacity. This value may underestimate the P sorption capacity of the soil in some cases. A larger value (~800mg/kg) can be used (Zanini et al., 1998) when the $P$ sorption capacity is underestimated. The concentration of P in the biozone is often reported low; thus, only the linear portion of a nonlinear isotherm is enough for estimation.

Phosphorus sorption isotherm described in Equation (15) gives an estimate of $P$ sorption capacity given the $P$ concentration and the distribution coefficient. According to this equation, effluent $P$ concentration leaching to sub-soil layer should be zero until the soil is saturated with$P$; however, small amount of soluble $P$ leaches to sub-soil layer with daily inflow of $P$ to the biozone. The effluent $P$ concentration is estimated by a linear relationship suggested by Bond et al. (2006) in which the outflow P concentration is proportional to the total amount P in the soil layer based on soil type as depicted in Figure 6:4-2.

After STE effluent passes through the biozone, it is discharged into the soil layers below, where the constituents are subject to additional transport/fate processes that are expected to occur in natural soils. Siegrist et al. (2005) have made the following assumptions to represent the physical system of biozone algorithm development based on laboratory observations and available knowledge.

1. Typical biozone thickness is 2-5 cm.

2. The biozone receives a continuous daily loading of STE. Intermittent dosing over the course of one day is not considered.

3. Depending on the climatic region, bermuda grass or similar type of grass is the typical vegetation seen above drainfield.

4. No STE inflow occurs if soil temperature goes below the freezing point and the biozone processes become inactive.

5. After a hydraulic failure, the model starts counting the failure days once the STE saturates the upper soil layers completely. Therefore, the initial several days after the hydraulic failure are not counted as failing period.

6. A failing system automatically turns to a fresh new active system after a user specified failing period (typically 2-3 months).

7. Total solids and TDS concentrations of STE are the same.

The biozone algorithm was integrated into the latest version of (SWAT+ version 2009 and the GIS interface version ArcSWAT+ 2009). Each unit of onsite septic systems is represented by a septic HRU. A septic HRU represents either one OWS or a hypothetical system that aggregates several OWSs that have similar HRU properties such as soil type, landuse and slope. In the latter case, septic HRUs in a subbasin are not spatially referenced and only lumped simulation is available. Septic HRUs are generated during the GIS interface processing based on the septic land use type defined by the users. Biozone and other properties of a septic system specific to a HRU is listed in septic input files (*.sep files). The septic water quality database developed based on literature lists STE loading rate and pollutants concentration in STE of the conventional system 25 kinds of advanced systems and an untreated effluent system (septwq.dat). Pollutant characteristics included in the water quality database are per capita septic effluent flow rate, Bio-chemical Oxygen Demand (BOD), Total Suspended Solid (TSS), Total Nitrogen (TN), ammonium, nitrite, nitrate, organic nitrogen, Total Phosphorus (TP), phosphate, organic P, and Fecal Coliform. A generic type for conventional system and advanced system are listed as the first and the second in the water quality database list, so the user can select these types with limited information. The hydrologic linkage to SWAT+ allows daily STE directly to be introduced to the biozone layer by increasing the soil moisture content in the subroutine percmain. The movement of soil water through percolation in and out of the biozne layer is then simulated by a SWAT+ soil water model in percmicro. The movement of nutrients through soil profile from the biozone layer is simulated by a combination of nutrient modules in SWAT+ and biozone. SWAT+ septic input parameters are presented in the updated SWAT+ Input/Output File Documentation (Neitsch et al., 2010).

Soil porosity (or saturated moisture content) is generally constant in natural soil; however, the porosity of biozone changes with time. The actual porosity of biozone decreases as the suspended solids from STE accumulate in the pore space and the mineralized biomass (dead body) increases in the biozone.

$\theta_s=\theta_{si}-\frac{plaque}{\rho_{bm}}$ (8)

where $\theta_{si}$ is initial soil porosity with zero plaque (mm). The moisture content at each time step is estimated using the mass balance of water within the biozone.

$\theta^t=\theta^{t-1}+\frac{Q_{STE}}{10*A_d}-I_p-ET-Q_{lat}$ (9)

where $ET$ is evaportranspiration from biozone (mm/day) and $Q_{lat}$ is lateral flow (mm/day). Percolation to a subsoil layer is triggered if moisture content exceeds the field capacity in the biozone layer. Potential percolation is the maximum amount of water that can percolate during the time interval.

$I_{p,pot}=K_{bz}*\Delta t$ (10)

where $I_{p,pot}$ is the potential amount of percolation (mm/day). The amount of water percolating to the sub-soil layer is calculated using storage routing methodology (Neitsch et al., 2005).

$I_{p,excess}=(\theta -\theta_f)(1-exp[\frac{-\Delta t}{TT_{perc}}])$ if $\theta > \theta_f$

$I_{p,excess}=0$ if $\theta \le \theta_f$ (11)

where $I_{p,excess}$ is the minimum amount of percolation (mm/day) and $TT_{perc}$ is travel time for percolation $(TT_{perc}=(\theta_s-\theta_f)/K_{bz})$ in hour. The actual percolation is the smaller of the potential and minimum percolation.

$I_p=min(I_{p,pot},I_{p,excess})$ (12)

In SWAT+ modeling, active systems are septic systems that are in operation or functioning as per guidelines and failing systems are septic systems that are subject to hydraulic failure and the effluent discharged differ from the standards. A typical service life span of an OWS ranges from 10 to 25 years depending on maintenance, pollutant loading rate, soil conditions and other factors. A septic HRU starting as an active system becomes a failing system as the biozone gets clogged by TSS and plaque of biomass. In the SWAT+ biozone module, biozone clogging, or hydraulic failure, is the main cause of system failure. There are other types of failing, but they are difficult to model in the SWAT+ structure and thus hydraulic failure was modeled as the only type of failing.

The schematic of the SWAT+ biozone algorithm is described in Figure 6:4-3. The biozone processes in septic HRUs are simulated on a daily basis as a subbasin-level. The biozone subroutine is called within the HRU iteration loop whenever the current HRU is septic. Each septic HRU is simulated based on a maintenance plan in which a failing system is assumed to be fixed and re-activated in two to three months time. In an active system, saturated water content gets lower as plaque fills in the soil pore space and field capacity increases as the amount of live bacteria biomass grows as filamentous biomass soaks up soil water. System failure occurs when the soil achieves saturated water content and field capacity. The model starts counting the number of days as the system fails and remains as a failing system. STE migrates to upper soil layers as water does not percolate through the bottom of the biozone layer. Depending on the thickness of soil above the biozone layer, it may take a few days to months until failure creates STE surface ponding. As STE migrates to upper soil layers, the nutrients in the STE are transported along with the STE. The amount of nutrients that transports to the upper soil layers is estimated based on the nutrient concentration in STE and the amount of water that migrates to the upper layer. There are no special treatment processes that apply to the nutrients in failing septic systems. If the number of days a septic system remained as failing exceeds the designed time, the failing system is updated to an active system and related properties are reinitialized as a fresh active system. An active system simulates soil water hydraulics and pollutants decay by executing the biozone equations described above.

Transformation and removal of pollutants in the biozone is directly related with the population of live bacteria biomass and bio-physical processes in the biozone layer. The fate of pollutants including Nitrogen, BOD, and Fecal Coliform is estimated by a first order reaction equation:

$C_{k,end}=C_{k,i}*e^{-K_k \Delta t}$ (13)

where $C_{k,end}$ is concentration of k constituent in the biozone at the end of the day (mg/L), $C_{k,i}$ is concentration of $k$ constituent in the biozone at the beginning of the day (mg/L), and $K_k$ is a first order reaction rate (1/day), which is a function of the total biomass of live bacteria and a reaction rate coefficient.

$K_k=[\frac{K_{1,k}*A_d}{\theta_{si}*V_{bz}}](Bio)$ (14)

where $K_{1,k}$ is the reaction rate calibration parameter for each constituent $k$ (m$^3$/kg) and $V_{bz}$ is the volume of the biozone (m$^3$). The various constituents included are meant for the primary reactions/processes that occur in the biozone layer such as nitrification, denitrification, BOD decay, and fecal coliform decay.

In Equation (14), the reaction rate is normalized with respect to the volume of bacterial biomass (pore volume,$\theta_{si} V_{bz}$) in the biozone layer as in the case of mortality and respiration equations. This normalization is done to avoid scaling issues in applying the algorithm to watershed scale simulations with parameters calibrated to small scale results (lab column tests).