Quantifying the impact of land management and land use on water supply and quality is a primary focus of environmental modeling. SWAT+ allows very detailed management information to be incorporated into a simulation.

The following five chapters review the methodology used by SWAT+ to simulate water management, tillage, urban areas, septic systems, and grassed waterways.

The plant operation initiates plant growth. This operation can be used to designate the time of planting for agricultural crops or the initiation of plant growth in the spring for a land cover that requires several years to reach maturity (forests, orchards, etc.).

The plant operation will be performed by SWAT+ only when no land cover is growing in an HRU. Before planting a new land cover, the previous land cover must be removed with a kill operation or a harvest and kill operation. If two plant operations are placed in the management file and the first land cover is not killed prior to the second plant operation, the second plant operation is ignored by the model.

Information required in the plant operation includes the timing of the operation (month and day or fraction of base zero potential heat units), the total number of heat units required for the land cover to reach maturity, and the specific land cover to be simulated in the HRU. If the land cover is being transplanted, the leaf area index and biomass for the land cover at the time of transplanting must be provided. Also, for transplanted land covers, the total number of heat units for the land cover to reach maturity should be from the period the land cover is transplanted to maturity (not from seed generation). Heat units are reviewed in Chapter 5:1.

The user has the option of varying the curve number in the HRU throughout the year. New curve number values may be entered in a plant operation, tillage operation and harvest and kill operation. The curve number entered for these operations are for moisture condition II. SWAT+ adjusts the entered value daily to reflect change in water content or plant evapotranspiration.

For simulations where a certain amount of crop yield and biomass is required, the user can force the model to meet this amount by setting a harvest index target and a biomass target. These targets are effective only if a harvest and kill operation is used to harvest the crop.

The grazing operation simulates plant biomass removal and manure deposition over a specified period of time. This operation is used to simulate pasture or range grazed by animals.

Information required in the grazing operation includes the time during the year at which grazing begins (month and day or fraction of plant potential heat units), the length of the grazing period, the amount of biomass removed daily, the amount of manure deposited daily, and the type of manure deposited. The amount of biomass trampled is an optional input.

Biomass removal in the grazing operation is similar to that in the harvest operation. However, instead of a fraction of biomass being specified, an absolute amount to be removed every day is given. In some conditions, this can result in a reduction of the plant biomass to a very low level that will result in increased erosion in the HRU. To prevent this, a minimum plant biomass for grazing may be specified (BIO_MIN). When the plant biomass falls below the amount specified for BIO_MIN, the model will not graze, trample, or apply manure in the HRU on that day.

If the user specifies an amount of biomass to be removed daily by trampling, this biomass is converted to residue.

Nutrient fractions and bacteria content of the manure applied during grazing must be stored in the fertilizer database. The manure nutrient and bacteria loadings are added to the topmost 10 mm of soil. This is the portion of the soil with which surface runoff interacts.

After biomass is removed by grazing and/or trampling, the plant’s leaf area index and accumulated heat units are set back by the fraction of biomass removed.

The harvest operation will remove plant biomass without killing the plant. This operation is most commonly used to cut hay or grass.

The only information required by the harvest operation is the date. However, a harvest index override and a harvest efficiency can be set.

When no harvest index override is specified, SWAT+ uses the plant harvest index from the plant growth database to set the fraction of biomass removed. The plant harvest index in the plant growth database is set to the fraction of the plant biomass partitioned into seed for agricultural crops and a typical fraction of biomass removed in a cutting for hay. If the user prefers a different fraction of biomass to be removed, the harvest index override should be set to the desired value.

A harvest efficiency may also be defined for the operation. This value specifies the fraction of harvested plant biomass removed from the HRU. The remaining fraction is converted to residue on the soil surface. If the harvest efficiency is left blank or set to zero, the model assumes this feature is not being used and removes 100% of the harvested biomass (no biomass is converted to residue).

After biomass is removed in a harvest operation, the plant’s leaf area index and accumulated heat units are set back by the fraction of biomass removed. Reducing the number of accumulated heat units shifts the plant’s development to an earlier period in which growth is usually occurring at a faster rate.

Management operations that control the plant growth cycle, the timing of fertilizer and pesticide and the removal of plant biomass are explained in this chapter. Water management and the simulation of urban areas are summarized in subsequent chapters.

A new operations input file was added that allows user to schedule management by Julian day and calendar year without considering cropping rotations and without using heat unit scheduling. Operations have been added for contouring, terracing, subsurface drains, filter strips, fire, and grass waterways.

Previous versions of SWAT+ only allowed the growth of one plant species at a time to be simulated. Algorithms from the ALMANAC model (Kiniry et al., 1992, Johnson et al., 2009) have been added to simulate multiple plant species growing and competing within a plant community. Plant communities that have been simulated include: crops and weeds, trees and grasses, different tree species in a boreal forest, and grasses and shrubs in rangeland communities. A data file is developed prior to simulation that describes the various plants within each community.

The harvest and kill operation stops plant growth in the HRU. The fraction of biomass specified in the land cover’s harvest index (in the plant growth database) is removed from the HRU as yield. The remaining fraction of plant biomass is converted to residue on the soil surface.

The only information required by the harvest and kill operation is the timing of the operation (month and day or fraction of plant potential heat units). The user also has the option of updating the moisture condition II curve number in this operation.

The kill operation stops plant growth in the HRU. All plant biomass is converted to residue.

The only information required by the kill operation is the timing of the operation (month and day or fraction of plant potential heat units).

Biological mixing is the redistribution of soil constituents as a result of the activity of biota in the soil (e.g. earthworms, etc.). Studies have shown that biological mixing can be significant in systems where the soil is only infrequently disturbed. In general, as a management system shifts from conventional tillage to conservation tillage to no-till there will be an increase in biological mixing. SWAT+ allows biological mixing to occur to a depth of 300 mm (or the bottom of the soil profile if it is shallower than 300 mm). The efficiency of biological mixing is defined by the user. The redistribution of nutrients by biological mixing is calculated using the same methodology as that used for a tillage operation.

The tillage operation redistributes residue, nutrients, pesticides and bacteria in the soil profile. Information required in the tillage operation includes the timing of the operation (month and day or fraction of base zero potential heat units), and the type of tillage operation.

The user has the option of varying the curve number in the HRU throughout the year. New curve number values may be entered in a plant operation, tillage operation and harvest and kill operation. The curve number entered for these operations are for moisture condition II. SWAT+ adjusts the entered value daily to reflect change in water content.

The mixing efficiency of the tillage implement defines the fraction of a residue/nutrient/pesticide/bacteria pool in each soil layer that is redistributed through the depth of soil that is mixed by the implement. To illustrate the redistribution of constituents in the soil, assume a soil profile has the following distribution of nitrate.

If this soil is tilled with a field cultivator, the soil will be mixed to a depth of 100 mm with 30% efficiency. The change in the distribution of nitrate in the soil is:

Because the soil is mixed to a depth of 100 mm by the implement, only the nitrate in the surface layer and layer 1 is available for redistribution. To calculated redistribution, the depth of the layer is divided by the tillage mixing depth and multiplied by the total amount of nitrate mixed. To calculate the final nitrate content, the redistributed nitrate is added to the unmixed nitrate for the layer.

All nutrient/pesticide/bacteria/residue pools are treated in the same manner as the nitrate example above. Bacteria mixed into layers below the surface layer is assumed to die.

Fertilization in an HRU may be scheduled by the user or automatically applied by SWAT+. When the user selects auto-application of fertilizer in an HRU, a nitrogen stress threshold must be specified. The nitrogen stress threshold is a fraction of potential plant growth. Anytime actual plant growth falls below this threshold fraction due to nitrogen stress, the model will automatically apply fertilizer to the HRU. The user specifies the type of fertilizer, the fraction of total fertilizer applied to the soil surface, the maximum amount of fertilizer that can be applied during the year, the maximum amount of fertilizer that can be applied in any one application, and the application efficiency.

To determine the amount of fertilizer applied, an estimate of the amount of nitrogen that will be removed in the yield is needed. For the first year of simulation, the model has no information about the amount of nitrogen removed from the soil by the plant. The nitrogen yield estimate is initially assigned a value using the following equations:

$yld_{est,N}=350*fr_{N,yld}*RUE$ if $HI_{opt}<1.0$ 6:1.8.1

$yld_{est,N}=1000*fr_{N,yld}*RUE$ if $HI_{opt}\ge 1.0$ 6:1.8.2

where $yld_{est,N}$ is the nitrogen yield estimate (kg N/ha), $fr_{N,yld}$ is the fraction of nitrogen in the yield, $RUE$ is the radiation-use efficiency of the plant (kg/ha⋅(MJ/m$^2$)$^{-1}$ or 10$^{-1}$g/MJ), and $HI_{opt}$ is the potential harvest index for the plant at maturity given ideal growing conditions. The nitrogen yield estimate is updated at the end of every simulation year using the equation:

$yld_{est,N}=\frac{yld_{est,Nprev}*yr_{sim}+yld_{yr,N}}{yr_{sim}+1}$ 6:1.8.3

where $yld_{est,N}$ is the nitrogen yield estimate update for the current year (kg N/ha), $yld_{est,Nprev}$ is the nitrogen yield estimate from the previous year (kg N/ha), $yr_{sim}$ is the year of simulation, $yld_{yr,N}$ is the nitrogen yield target for the current year (kg N/ha). The nitrogen yield target for the current year is calculated at the time of harvest using the equation:

$yld_{yr,N}=bio_{ag}*fr_N*fert_{eff}$ 6:1.8.4

where $yld_{yr,N}$ is the nitrogen yield target for the current year (kg N/ha), $bio_{ag}$ is the aboveground biomass on the day of harvest (kg ha$^{-1}$), $fr_N$ is the fraction of nitrogen in the plant biomass calculated with equation 5:2.3.1, and $fert_{eff}$ is the fertilizer application efficiency assigned by the user. The fertilizer application efficiency allows the user to modify the amount of fertilizer applied as a function of plant demand. If the user would like to apply additional fertilizer to adjust for loss in runoff, $fert_{eff}$ will be set to a value greater than 1. If the user would like to apply just enough fertilizer to meet the expected demand, $fert_{eff}$ will be set to 1. If the user would like to apply only a fraction of the demand, $fert_{eff}$ will be set to a value less than 1.

The optimal amount of mineral nitrogen to be applied is calculated:

$minN_{app}=yld_{est,N}-(NO3+NH4)-bio_N$ 6:1.8.5

where $minN_{app}$ is the amount of mineral nitrogen applied (kg N/ha), $yld_{est,N}$ is the nitrogen yield estimate (kg N/ha), $NO3$ is the nitrate content of the soil profile (kg NO$_3$-N/ha), NH4 is the ammonium content of the soil profile (kg NH$_4$-N/ha), and $bio_N$ is the actual mass of nitrogen stored in plant material (kg N/ha). If the amount of mineral nitrogen calculated with equation 6:1.8.5 exceeds the maximum amount allowed for any one application, $minN_{app}$ is reset to the maximum value ($minN_{app}=minN_{app,mx}$). The total amount of nitrogen applied during the year is also compared to the maximum amount allowed for the year.

Once the amount applied reaches the maximum amount allowed for the year ($minN_{app,mxyr}$), SWAT+ will not apply any additional fertilizer regardless of nitrogen stress. Once the amount of mineral nitrogen applied is determined, the total amount of fertilizer applied is calculated by dividing the mass of mineral nitrogen applied by the fraction of mineral nitrogen in the fertilizer:

$fert=\frac{minN_{app}}{fert_{minN}}$ 6:1.8.6

where $fert$ is the amount of fertilizer applied (kg/ha), $minN_{app}$ is the amount of mineral nitrogen applied (kg N/ha), and $fert_{minN}$ is the fraction of mineral nitrogen in the fertilizer.

The type of fertilizer applied in the HRU is specified by the user. In addition to mineral nitrogen, organic nitrogen and phosphorus and mineral phosphorus are applied to the HRU. If a manure is applied, bacteria loadings to the HRU are also determined. The amount of each type of nutrient and bacteria is calculated from the amount of fertilizer and fraction of the various nutrient types in the fertilizer as summarized in Section 6:1.7.

While the model does not allow fertilizer to be applied as a function of phosphorus stress, the model does monitor phosphorus stress in the auto-fertilization subroutine. If phosphorus stress causes plant growth to fall below 75% of potential growth, the model ignores the fraction of mineral phosphorus in the fertilizer and applies an amount of mineral phosphorus equal to ($\frac{1}{7}*minN_{app}$).

Table 6:1-8: SWAT+ input variables that pertain to auto-fertilization.

Accurately reproducing water management practices can be one of the most complicated portions of data input for the model. Because water management affects the hydrologic balance, it is critical that the model is able to accommodate a variety of management practices. Water management options modeled by SWAT+ include irrigation, tile drainage, impounded/depressional areas, water transfer, consumptive water use, and loadings from point sources.

Irrigation in an HRU may be scheduled by the user or automatically applied by SWAT+ in response to a water deficit in the soil. In addition to specifying the timing and application amount, the user must specify the source of irrigation water.

Water applied to an HRU is obtained from one of five types of water sources: a reach, a reservoir, a shallow aquifer, a deep aquifer, or a source outside the watershed. In addition to the type of water source, the model must know the location of the water source (unless the source is outside the watershed). For the reach, shallow aquifer or deep aquifer, SWAT+ needs to know the reach number or subbasin number, respectively, in which the source is located. If a reservoir is used to supply water, SWAT+ must know the reservoir number.

If the source of the irrigation water is a reach, SWAT+ allows additional input parameters to be set. These parameters are used to prevent flow in the reach from being reduced to zero as a result of irrigation water removal. Users may define a minimum in-stream flow, a maximum irrigation water removal amount that cannot be exceeded on any given day, and/or a fraction of total flow in the reach that is available for removal on a given day.

For a given irrigation event, SWAT+ determines the amount of water available in the source. The amount of water available is compared to the amount of water specified in the irrigation operation. If the amount available is less than the amount specified, SWAT+ will only apply the available water.

A manual irrigation application can be scheduled by date or by heat units. Irrigation amount (mm), input by the user, is the amount of water applied that reaches the soil. An irrigation efficiency factor is applied to account for losses from the source to the soil including conveyance loss and evaporative loss. The surface runoff ratio is the fraction of water applied that leaves the field as surface runoff. The remainder infiltrates into the soil and is subject to the soil water routing algorithms described in Section 2, Chapter 3. This allows for more realistic simulation of the soil water profile and application of excess irrigation for leaching salts.

A primary mechanism of disposal for manure generated by intensive animal operations such as confined animal feedlots is the land application of waste. In this type of a land management system, waste is applied every few days to the fields. Using the continuous fertilization operation allows a user to specify the frequency and quantity of manure applied to an HRU without the need to insert a fertilizer operation in the management file for every single application.

The continuous fertilizer operation requires the user to specify the beginning date of the continuous fertilization period, the total length of the fertilization period, and the number of days between individual fertilizer/manure applications. The amount of fertilizer/manure applied in each application is specified as well as the type of fertilizer/manure.

Nutrients and bacteria in the fertilizer/manure are applied to the soil surface. Unlike the fertilization operation or auto-fertilization operation, the continuous fertilization operation does not allow the nutrient and bacteria loadings to be partitioned between the surface 10 mm and the part of the 1st soil layer underlying the top 10 mm. Everything is added to the top 10 mm, making it available for transport by surface runoff.

Nutrient and bacteria loadings to the HRU are calculated using the equations reviewed in Section 6:1.7.

Table 6:1-9: SWAT+ input variables that pertain to continuous fertilization.

The pesticide operation applies pesticide to the HRU.

Information required in the pesticide operation includes the timing of the operation (month and day or fraction of plant potential heat units), the type of pesticide applied, and the amount of pesticide applied.

Field studies have shown that even on days with little or no wind, a portion of pesticide applied to the field is lost. The fraction of pesticide that reaches the foliage or soil surface is defined by the pesticide’s application efficiency. The amount of pesticide that reaches the foliage or ground is:

$pest'=ap_{ef}*pest$ 6:1.10.1

where $pest'$ is the effective amount of pesticide applied (kg pst/ha), $ap_{ef}$ is the pesticide application efficiency, and pest is the actual amount of pesticide applied (kg pst/ha).

The amount of pesticide reaching the ground surface and the amount of pesticide added to the plant foliage is calculated as a function of ground cover. The ground cover provided by plants is:

$gc=\frac{1.99532-erfc[1.333*LAI-1]}{2.1}$ 6:1.10.2

where $gc$ is the fraction of the ground surface covered by plants, $erfc$ is the complementary error function, and $LAI$ is the leaf area index.

The complementary error function frequently occurs in solutions to advective-dispersive equations. Values for $erfc(\beta)$ and $erf(\beta)$ (erf is the error function for $\beta$), where $\beta$ is the argument of the function, are graphed in Figure 6:1-1. The figure shows that $erf(\beta)$ ranges from –1 to +1 while $erfc(\beta)$ ranges from 0 to +2. The complementary error function takes on a value greater than 1 only for negative values of the argument.

Once the fraction of ground covered by plants is known, the amount of pesticide applied to the foliage is calculated:

$pest_{fol}=gc*pest'$ 6:1.10.3

and the amount of pesticide applied to the soil surface is

$pest_{surf}=(1-gc)*pest'$ 6:1.10.4

where $pest_{fol}$ is the amount of pesticide applied to foliage (kg pst/ha), $pest_{surf}$ is the amount of pesticide applied to the soil surface (kg pst/ha), $gc$ is the fraction of the ground surface covered by plants, and $pest'$ is the effective amount of pesticide applied (kg pst/ha).

Table 6:1-10: SWAT+ input variables that pertain to pesticide application.

The fertilizer operation applies fertilizer or manure to the soil.

Information required in the fertilizer operation includes the timing of the operation (month and day or fraction of plant potential heat units), the type of fertilizer/manure applied, the amount of fertilizer/manure applied, and the depth distribution of fertilizer application.

SWAT+ assumes surface runoff interacts with the top 10 mm of soil. Nutrients contained in this surface layer are available for transport to the main channel in surface runoff. The fertilizer operation allows the user to specify the fraction of fertilizer that is applied to the top 10 mm. The remainder of the fertilizer is added to the first soil layer defined in the HRU .sol file.

In the fertilizer database, the weight fraction of different types of nutrients and bacteria are defined for the fertilizer. The amounts of nutrient added to the different pools in the soil are calculated:

$NO3_{fert}=fert_{minN}*(1-fert_{NH4})*fert$ 6:1.7.1

$NH4_{fert}=fert_{minN}*fert_{NH4}*fert$ 6:1.7.2

$orgN_{frsh,fert}=0.5*fert_{orgN}*fert$ 6:1.7.3

$orgN_{act,fert}=0.5*fert_{orgN}*fert$ 6:1.7.4

$P_{solution,fert}=fert_{minP}*fert$ 6:1.7.5

$orgP_{frsh,fert}=0.5*fert_{orgP}*fert$ 6:1.7.6

$orgP_{hum,fert}=0.5*fert_{orgP}*fert$ 6:1.7.7

where $NO3_{fert}$ is the amount of nitrate added to the soil in the fertilizer (kg N/ha), $NH4_{fert}$ is the amount of ammonium added to the soil in the fertilizer (kg N/ha), $orgN_{frsh,fert}$ is the amount of nitrogen in the fresh organic pool added to the soil in the fertilizer (kg N/ha), $orgN_{act,fert}$ is the amount of nitrogen in the active organic pool added to the soil in the fertilizer (kg N/ha), $P_{solution,fert}$ is the amount of phosphorus in the solution pool added to the soil in the fertilizer (kg P/ha), $orgP_{frsh,fert}$ is the amount of phosphorus in the fresh organic pool added to the soil in the fertilizer (kg P/ha), $orgP_{hum,fert}$ is the amount of phosphorus in the humus organic pool added to the soil in the fertilizer (kg P/ha), $fert_{minN}$ is the fraction of mineral N in the fertilizer, $fert_{NH4}$ is the fraction of mineral $N$ in the fertilizer that is ammonium, $fert_{orgN}$ is the fraction of organic $N$ in the fertilizer, $fert_{minP}$ is the fraction of mineral $P$ in the fertilizer, $fert_{orgP}$ is the fraction of organic $P$ in the fertilizer, and $fert$ is the amount of fertilizer applied to the soil (kg/ha).

If manure is applied, the bacteria in the manure may become attached to plant foliage or be incorporated into the soil surface layer during application. The amount of bacteria reaching the ground surface and the amount of bacteria adhering to the plant foliage is calculated as a function of ground cover. The ground cover provided by plants is:

$gc=\frac{1.99532-erfc[1.333*LAI-2]}{2.1}$ 6:1.7.8

where $gc$ is the fraction of the ground surface covered by plants, $erfc$ is the complementary error function, and $LAI$ is the leaf area index.

The complementary error function frequently occurs in solutions to advective-dispersive equations. Values for $erfc(\beta)$ and $erf(\beta)$ ($erf$ is the error function for $\beta$), where $\beta$ is the argument of the function, are graphed in Figure 6:1-1. The figure shows that $erf(\beta)$ ranges from –1 to +1 while $erfc(\beta)$ ranges from 0 to +2. The complementary error function takes on a value greater than 1 only for negative values of the argument.

Once the fraction of ground covered by plants is known, the amount of bacteria applied to the foliage is calculated:

$bact_{lp,fol}=\frac{gc*fr_{active}*fert_{lpbact}*fert}{10}$ 6:1.7.9

$bact_{p,fol}=\frac{gc*fr_{active}*fert_{pbact}*fert}{10}$ 6:1.7.10

and the amount of bacteria applied to the soil surface is

$bact_{lpsol,fert}=\frac{(1-gc)*fr_{active}*fert_{lpbact}*k_{bact}*fert}{10}$ 6:1.7.11

$bact_{lpsorb,fert}=\frac{(1-gc)*fr_{active}*fert_{lpbact}*(1-k_{bact})*fert}{10}$ 6:1.7.12

$bact_{psol,fert}=\frac{(1-gc)*fr_{active}*fert_{pbact}*k_{bact}*fert}{10}$ 6:1.7.13

$bact_{psorb,fert}=\frac{(1-gc)*fr_{active}*fert_{pbact}*(1-k_{bact})*fert}{10}$ 6:1.7.14

where $bact_{lp,fol}$ is the amount of less persistent bacteria attached to the foliage (# cfu/m$^2$), $bact_{p,fol}$ is the amount of persistent bacteria attached to the foliage (# cfu/m$^2$), $bact_{lpsol,fert}$ is the amount of less persistent bacteria in the solution pool added to the soil (# cfu/m$^2$), $bact_{lpsorb,fert}$ is the amount of less persistent bacteria in the sorbed pool added to the soil (# cfu/m$^2$), $bact_{psol,fert}$ is the amount of persistent bacteria in the solution pool added to the soil (# cfu/m$^2$), $bact_{psorb,fert}$ is the amount of persistent bacteria in the sorbed pool added to the soil (# cfu/m$^2$), $gc$ is the fraction of the ground surface covered by plants, $fr_{active}$ is the fraction of the manure containing active colony forming units, $fert_{lpbact}$ is the concentration of less persistent bacteria in the fertilizer (# cfu/g manure), $fert_{pbact}$ is the concentration of persistent bacteria in the fertilizer (# cfu/g manure), $k_{bact}$ is the bacterial partition coefficient, and $fert$ is the amount of fertilizer/manure applied to the soil (kg/ha).

Table 6:1-7: SWAT+ input variables that pertain to fertilizer application.

Impounded/depressional areas are simulated as a water body overlying a soil profile in an HRU. This type of ponded system is needed to simulate the growth of rice, cranberries or any other plant that grows in a waterlogged system. The simulation and management operations pertaining to impounded/depressional areas are reviewed in Chapter 8:1.

SWAT+ directly simulates the loading of water, sediment and other constituents off of land areas in the watershed. To simulate the loading of water and pollutants from sources not associated with a land area (e.g. sewage treatment plants), SWAT+ allows point source information to be read in at any point along the channel network. The point source loadings may be summarized on a daily, monthly, yearly, or average annual basis.

Files containing the point source loads are created by the user. The loads are read into the model and routed through the channel network using rechour, recday, recmon, recyear, or reccnst commands in the watershed configuration file. SWAT+ will read in water, sediment, organic N, organic P, nitrate, soluble P, ammonium, nitrite, metal, and bacteria data from the point source files. Chapter 2 in the SWAT+ User’s Manual reviews the format of the command lines in the watershed configuration file while Chapter 31 in the SWAT+ User’s Manual reviews the format of the point source files.

When the user selects auto-application of irrigation, the application can be triggered by a water stress threshold or a soil water deficit threshold. The water stress threshold is a fraction of potential plant growth. Any day actual plant stress falls below this threshold fraction due to water stress, the model will automatically apply water up to a maximum amount per application as input by the user. An irrigation application can also be triggered by a soil water deficit threshold. When total soil water in the profile falls below field capacity by more than the soil water deficit threshold, an irrigation application occurs. As with a manual application, a maximum irrigation amount, irrigation efficiency and surface runoff ratio are applied to each application.

Table 6:2-1: SWAT+ input variables that pertain to irrigation.

Consumptive water use is a management tool that removes water from the basin. Water removed for consumptive use is considered to be lost from the system. SWAT+ allows water to be removed from the shallow aquifer, the deep aquifer, the reach or the pond within any subbasin in the watershed. Water also may be removed from reservoirs for consumptive use.

Consumptive water use is allowed to vary from month to month. For each month in the year, an average daily volume of water removed from the source is specified. For reservoirs, the user may also specify a fraction of the water removed that is lost during removal. The water lost in the removal process becomes outflow from the reservoir.

Table 6:2-4: SWAT+ input variables that pertain to consumptive water use.

While water is most typically removed from a water body for irrigation purposes, SWAT+ also allows water to be transferred from one water body to another. This is performed with a transfer command in the watershed configuration file.

The transfer command can be used to move water from any reservoir or reach in the watershed to any other reservoir or reach in the watershed. The user must input the type of water source, the location of the source, the type of water body receiving the transfer, the location of the receiving water body, and the amount of water transferred.

Three options are provided to specify the amount of water transferred: a fraction of the volume of water in the source; a volume of water left in the source; or the volume of water transferred. The transfer is performed every day of the simulation.

The transfer of water from one water body to another can be accomplished using other methods. For example, water could be removed from one water body via consumptive water use and added to another water body using point source files.

Table 6:2-3: SWAT+ input variables that pertain to water transfer.

To simulate tile drainage in an HRU, the user must specify the depth from the soil surface to the drains, the amount of time required to drain the soil to field capacity, and the amount of lag between the time water enters the tile till it exits the tile and enters the main channel.

Tile drainage occurs when the perched water table rises above the depth at which the tile drains are installed. The amount of water entering the drain on a given day is calculated:

$tile_{wtr}=\frac{h_{wtbl}-h_{drain}}{h_{wtbl}}*(SW-FC)*(1-exp[\frac{-24}{t_{drain}}])$ if $h_{wtbl}>h_{drain}$ 6:2.2.1

where $tile_{wtr}$ is the amount of water removed from the layer on a given day by tile drainage (mm H$_2$O), $h_{wtbl}$ is the height of the water table above the impervious zone (mm), $h_{drain}$ is the height of the tile drain above the impervious zone (mm), $SW$ is the water content of the profile on a given day (mm H$_2$O), $FC$ is the field capacity water content of the profile (mm H$_2$O), and $t_{drain}$ is the time required to drain the soil to field capacity (hrs).

Water entering tiles is treated like lateral flow. The flow is lagged using equations reviewed in Chapter 2:3.

Table 6:2-2: SWAT+ input variables that pertain to tile drainage.

The linear regression models incorporated into SWAT+ are those described by Driver and Tasker (1988). The regression models were developed from a national urban water quality database that related storm runoff loads to urban physical, land use, and climatic characteristics. USGS developed these equations to predict loadings in ungaged urban watersheds.

The regression models calculate loadings as a function of total storm rainfall, drainage area and impervious area. The general equation is

$Y=\frac{\beta_0*(R_{day}/25.4)^{\beta_1}*(DA*imp_{tot}/2.59)^{\beta_2}*(imp_{tot}*100+1)^{\beta_3}*\beta_4}{2.205}$ 6:3.3.1

where $Y$ is the total constituent load (kg), $R_{day}$ is precipitation on a given day (mm H$_2$O), $DA$ is the HRU drainage area (km$^2$), $imp_{tot}$ is the fraction of the total area that is impervious, and the $\beta$ variables are regression coefficients. The regression equations were developed in English units, so conversion factors were incorporated to adapt the equations to metric units: 25.4 mm/inch, 2.59 km2/mi2, and 2.205 lb/kg.

USGS derived three different sets of regression coefficients that are based on annual precipitation. Category I coefficients are used in watersheds with less than 508 mm of annual precipitation. Category II coefficients are used in watersheds with annual precipitation between 508 and 1016 mm. Category III coefficients are used in watersheds with annual precipitation greater than 1016 mm. SWAT+ determines the annual precipitation category for each subbasin by summing the monthly precipitation totals provided in the weather generator input file.

Regression coefficients were derived to estimate suspended solid load, total nitrogen load, total phosphorus load and carbonaceous oxygen demand (COD). SWAT+ calculates suspended solid, total nitrogen, and total phosphorus loadings (the carbonaceous oxygen demand is not currently calculated). Regression coefficients for these constituents are listed in Table 6:3-3.

Once total nitrogen and phosphorus loads are calculated, they are partitioned into organic and mineral forms using the following relationships from Northern Virginia Planning District Commission (1979). Total nitrogen loads consist of 70 percent organic nitrogen and 30 percent mineral (nitrate). Total phosphorus loads are divided into 75 percent organic phosphorus and 25 percent orthophosphate.

Table 6:3-4: SWAT+ input variables that pertain to urban modeling with regression equations.

Most large watersheds and river basins contain areas of urban land use. Estimates of the quantity and quality of runoff in urban areas are required for comprehensive management analysis. SWAT+ calculates runoff from urban areas with the SCS curve number method or the Green & Ampt equation. Loadings of sediment and nutrients are determined using one of two options. The first is a set of linear regression equations developed by the USGS (Driver and Tasker, 1988) for estimating storm runoff volumes and constituent loads. The other option is to simulate the buildup and washoff mechanisms, similar to SWMM - Storm Water Management Model (Huber and Dickinson, 1988).

Urban areas differ from rural areas in the fraction of total area that is impervious. Construction of buildings, parking lots and paved roads increases the impervious cover in a watershed and reduces infiltration. With development, the spatial flow pattern of water is altered and the hydraulic efficiency of flow is increased through artificial channels, curbing, and storm drainage and collection systems. The net effect of these changes is an increase in the volume and velocity of runoff and larger peak flood discharges.

Impervious areas can be differentiated into two groups: the area that is hydraulically connected to the drainage system and the area that is not directly connected. As an example, assume there is a house surrounded by a yard where runoff from the roof flows into the yard and is able to infiltrate into the soil. The rooftop is impervious but it is not hydraulically connected to the drainage system. In contrast, a parking lot whose runoff enters a storm water drain is hydraulically connected. Table 6:3-1 lists typical values for impervious and directly connected impervious fractions in different urban land types.

During dry periods, dust, dirt and other pollutants build up on the impervious areas. When precipitation events occur and runoff from the impervious areas is generated, the runoff will carry the pollutants as it moves through the drainage system and enters the channel network of the watershed.

In an impervious area, dust, dirt and other constituents are built up on street surfaces in periods of dry weather preceding a storm. Build up may be a function of time, traffic flow, dry fallout and street sweeping. During a storm runoff event, the material is then washed off into the drainage system. Although the build up/wash off option is conceptually appealing, the reliability and credibility of the simulation may be difficult to establish without local data for calibration and validation (Huber and Dickinson, 1988).

When the build up/wash off option is used in SWAT+, the urban hydrologic response unit (HRU) is divided into pervious and impervious areas. Management operations other than sweep operations are performed in the pervious portion of the HRU. Sweep operations impact build up of solids in the impervious portion of the HRU. For the pervious portion of the HRU, sediment and nutrient loadings are calculated using the methodology summarized in Chapters 4:1 and 4:2. The impervious portion of the HRU uses the build up/wash off algorithm to determine sediment and nutrient loadings.

The build up/wash off algorithm calculates the build up and wash off of solids. The solids are assumed to possess a constant concentration of organic and mineral nitrogen and phosphorus where the concentrations are a function of the urban land type.

Build up of solids is simulated on dry days with a Michaelis-Menton equation:

$SED=\frac{SED_{mx}*td}{(t_{half}+td)}$ 6:3.4.1

where $SED$ is the solid build up (kg/curb km) $td$ days after the last occurrence of $SED=0$ kg/curb km, $SED_{mx}$ is the maximum accumulation of solids possible for the urban land type (kg/curb km), and $t_{half}$ is the length of time needed for solid build up to increase from 0 kg/curb km to $\frac{1}{2} SED_{mx}$(days). A dry day is defined as a day with surface runoff less than 0.1 mm. An example build-up curve is shown in Figure 6:3-1. As can be seen from the plot, the Michaelis-Menton function will initially rise steeply and then approach the asymptote slowly.

The two parameters that determine the shape of this curve are $SED_{mx}$ and $t_{half}$. These parameters are a function of the urban land type.

Wash off is the process of erosion or solution of constituents from an impervious surface during a runoff event. An exponential relationship is used to simulate the wash off process (Huber and Dickinson, 1988):

$Y_{sed}=SED_0*(1-e^{-kk*t})$ 6:3.4.2

where $Y_{sed}$ is the cumulative amount of solids washed off at time $t$ (kg/curb km), $SED_0$ is the amount of solids built up on the impervious area at the beginning of the precipitation event (kg/curb km), and $kk$ is a coefficient.

The coefficient, $kk$, may be estimated by assuming it is proportional to the peak runoff rate:

$kk=urb_{coef}*q_{peak}$ 6:3.4.3

where $urb_{coef}$ is the wash off coefficient (mm$^{-1}$) and $q_{peak}$ is the peak runoff rate (mm/hr).

The original default value for $urb_{coef}$ was calculated as 0.18 mm-1 by assuming that 13 mm of total runoff in one hour would wash off 90% of the initial surface load. Later estimates of $urb_{coef}$ gave values ranging from 0.002-0.26 mm$^{-1}$. Huber and Dickinson (1988) noted that values between 0.039 and 0.390 mm-1 for $urb_{coef}$ give sediment concentrations in the range of most observed values. They also recommended using this variable to calibrate the model to observed data.

To convert the sediment loading from units of kg/curb km to kg/ha, the amount of sediment removed by wash off is multiplied by the curb length density. The curb length density is a function of the urban land type. Nitrogen and phosphorus loadings from the impervious portion of the urban land area are calculated by multiplying the concentration of nutrient by the sediment loading.

Street cleaning is performed in urban areas to control buildup of solids and trash. While it has long been thought that street cleaning has a beneficial effect on the quality of urban runoff, studies by EPA have found that street sweeping has little impact on runoff quality unless it is performed every day (U.S. Environmental Protection Agency, 1983).

SWAT+ performs street sweeping operations only when the build up/wash off algorithm is specified for urban loading calculations. Street sweeping is performed only on dry days, where a dry day is defined as a day with less than 0.1 mm of surface runoff. The sweeping removal equation (Huber and Dickinson, 1988) is:

$SED=SED_0*(1-fr_{av}*reff)$ 6:3.4.4

where $SED$ is amount of solids remaining after sweeping (kg/curb km), $SED_0$ is the amount of solids present prior to sweeping (kg/curb km), $fr_{av}$ is the fraction of the curb length available for sweeping (the availability factor), and $reff$ is the removal efficiency of the sweeping equipment. The availability factor and removal efficiency are specified by the user.

The availability factor, $fr_{av}$, is the fraction of the curb length that is sweepable. The entire curb length is often not available for sweeping due to the presence of cars and other obstacles.

The removal efficiency of street sweeping is a function of the type of sweeper, whether flushing is a part of the street cleaning process, the quantity of total solids, the frequency of rainfall events and the constituents considered. Removal efficiency can vary depending on the constituent being considered, with efficiencies being greater for particulate constituents. The removal efficiencies for nitrogen and phosphorus are typically less than the solid removal efficiency (Pitt, 1979). Because SWAT+ assumes a set concentration of nutrient constituents in the solids, the same removal efficiency is in effect used for all constituents. Table 6:3-5 provides removal efficiencies for various street cleaning programs.

Table 6:3-6: SWAT+ input variables that pertain to build up/wash off.

In urban areas, surface runoff is calculated separately for the directly connected impervious area and the disconnected impervious/pervious area. For directly connected impervious areas, a curve number of 98 is always used. For disconnected impervious/pervious areas, a composite curve number is calculated and used in the surface runoff calculations. The equations used to calculate the composite curve number for disconnected impervious/pervious areas are (Soil Conservation Service Engineering Division, 1986):

$CN_c=CN_p+imp_{tot}*(CN_{imp}-CN_p)*(1-\frac{imp_{dcon}}{2*imp_{tot}})$ if $imp_{tot}<0.30$ 6:3.2.1

$CN_c=CN_p+imp_{tot}*(CN_{imp}-CN_p)$ if $imp_{tot}>0.30$ 6:3.2.2

where $CN_c$ is the composite moisture condition II curve number, $CN_p$ is the pervious moisture condition II curve number, $CN_{imp}$ is the impervious moisture condition II curve number, $imp_{tot}$ is the fraction of the HRU area that is impervious (both directly connected and disconnected), and $imp_{dcon}$ is the fraction of the HRU area that is impervious but not hydraulically connected to the drainage system.

The fraction of the HRU area that is impervious but not hydraulically connected to the drainage system, $imp_{dcon}$, is calculated

$imp_{dcon}=imp_{tot}-imp_{con}$ 6:3.2.3

where $imp_{tot}$ is the fraction of the HRU area that is impervious (both directly connected and disconnected), and $imp_{con}$ is the fraction of the HRU area that is impervious and hydraulically connected to the drainage system.

Table 6:3-2: SWAT+ input variables that pertain to surface runoff calculations in urban areas.

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).

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).

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:

$\frac{d(Bio)}{dt}=\alpha*[\sum Q_{STE}*C_{BOD,in}-I_p*C_{BOD}]-R_{resp}-R_{mort}-R_{slough}$ (1)

where $Bio$ is the amount of live bacteria biomass in biozone (kg/ha), $C_{BOD}$,in is the $BOD$concentration in the $STE$ (mg/L), $C_{BOD}$ is the $BOD$ concentration in biozone (mg/L), $\alpha$ is gram of live bacteria growth to gram of $BOD$ in STE (conversion factor), $Q_{STE}$ is the flow rate of $STE$ (m$^3$/day), Ip is the amount of percolation out of the biozone (m$^3$/day), $R_{resp}$ is the amount of respiration of bacteria (kg/ha), $R_{mort}$ is the amount of mortality of bacteria (kg/ha), and $R_{slough}$ is the amount of sloughed off bacteria (kg/ha).

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:

(6)

where is field capacity at the end of the day (mm), is field capacity at the beginning of the day (mm), is saturated moisture content at the beginning of the day (mm), is the density of live bacterial biomass (~1000 kg/m3), Φ is field capacity coefficient 1 (unitless), and is field capacity coefficient 2 (unitless).

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.

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)

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).

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.

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.

A widely used conservation practice to remove agricultural and urban pollutants before reaching nearby water bodies is the vegetative filter strip (VFS). A VFS is a strip of dense vegetation located to intercept runoff from upslope pollutant sources and filter it. The previous version of SWAT+ contained a VFS algorithm, but is had some limitations. It used the same filtering efficiency for sediment and all nutrient forms. Differing trapping efficiencies have been observed between soluble and particulate nutrients (Goel et al., 2004). In the previous version of SWAT+ the VFS model does not consider the effects of flow concentration apparent at the field and watershed scales. Due to widespread use of the SWAT+ to simulate VFSs (Chu et al., 2005; Arabi et al., 2008; Parajuli et al., 2008) improvements in these routines were needed.

A model to predict filter strip effectiveness under ideal uniform sheet flow conditions was developed from Vegetative Filter Strip MODel (VFSMOD) (Muñoz-Carpena, 1999) and measured data derived from 22 published studies. These studies were identified from a general search of the literature and other published summaries of VFS or riparian buffer effectiveness (Wenger, 1999; Helmers, 2003; Parkyn, 2004; Krutz et al., 2005; Mayer et al., 2005; Dorioz et al., 2006; Koelsch et al., 2006).

The filter strip model was adapted to operate at the field scale by considering the effects of flow concentration generally absent from plot scale experiments. Flow distribution through ten hypothetical filer strips was evaluated using high resolution (2m) topographical data and multipath flow accumulation (Quinn et al., 1991). Significant flow concentration was predicted at all sites, on average 10% of the filter strip received half of the field runoff. As implemented in SWAT+, the filter strip model contains two sections, a large section receiving relatively modest flow densities and a smaller section treating more concentrated flow. The combined model was incorporated into SWAT+ and verified for proper function. A full description of the filter strip model is presented below. A theoretical approach due to a lack of measured effectiveness data is used for grassed waterways. The model includes separate algorithms for the submerged and unsubmerged portions of the waterway. Particulate trapping in the submerged portion is based on the same sediment transport capacity algorithms employed in SWAT+’s channel reaches. The unsubmerged portion of the waterway is treated as a simplified filter strip. Runoff which enters laterally along the length of the waterway is subjected to this additional filtering effect. Larger events submerge a larger fraction of the waterway leaving less area to filter incoming runoff. Channel geometry for grassed waterways is defined as trapezoidal with 8:1 side slopes. Length, width, depth, and slope are required to simulate waterways in SWAT+.

To evaluate the effectiveness of VFSs under ideal conditions, a model was developed from a combination of measured data derived from literature and Vegetative Filter Strip MODel (VFSMOD) (Muñoz-Carpena, 1999) simulations. VFSMOD was selected for this application over other VFS-related models due to its process-based nature, abundant documentation, and ease of use. The algorithms used in VFSMOD are complex, requiring iterative solutions and significant computational resources. A watershed scale model, which may simulate many hundreds of VFS daily for decades, requires a less computationally intensive solution. VFSMOD was, therefore, not a candidate for incorporation into SWAT+.

VFSMOD model and its companion program, UH, were used to generate a database of 1650 VFS simulations. The UH utility uses the curve number approach (USDA-SCS, 1972), unit hydrograph and the Modified Universal Soil Loss Equation (MUSLE) (Williams, 1975) to generate synthetic sediment and runoff loads from a source area upslope of the VFS (Muñoz-Carpena and Parsons, 2005). This simulation database contained 3 h rainfall events ranging from 10 mm to 100 mm, on a cultivated field with a curve number of 85 and a C factor of 0.1. Field dimensions were fixed at 100 m by 10 m with a 10 m wide VFS at the downslope end. Width of the VFS ranged from 1 m to 20 m yielding drainage area to VFS area ratios from 5 to 100. Slopes of 2%, 5% and 10% were simulated on 11 soil textural classes. This database was generated via software, which provided input parameters to both UH and VFSMOD then executed each program in turn. This database and a variety of other VFSMOD simulations were used to evaluate the sensitivity of various parameters and correlations between model inputs and predictions.

An empirical model for runoff reduction by VFSs was developed based on VFSMOD simulations. The model was derived from runoff loading and saturated hydraulic conductivity using the statistical package Minitab 15 (Minitab-Inc., 2006). Saturated hydraulic conductivity is available in SWAT+, and runoff loading can be calculated from HRU-predicted runoff volume and drainage area to VFS area ratio (DAFSratio). Both independent variables were transformed to improve the regression. The final form is given below:

$R_R=75.8-10.8(ln(R_L)+25.9 ln(K_{SAT})$ 6:5.1.1

where $R_R$ is the runoff reduction (%); $R_L$ is the runoff loading (mm); and $K_{SAT}$ is the saturated hydraulic conductivity (mm/hr). The regression was able to explain the majority of the variability (R$^2$ = 0.76; n = 1,650) in the simulated runoff reduction. The resulting model (Figure 6:5-1) produced runoff reduction efficiencies from -30% to 160%. Reductions greater than 100% are not possible; these were an artifact of the regression model. VFSs in SWAT+ were not allowed to generate additional runoff or pollutants; the model was limited to a range of 0% to 100%. The comparison between the empirical model and VFSMOD simulations improved (R$^2$ = 0.84) when the range was limited.

The sediment reduction model developed for SWAT+ was based on measured VFS data. A VFS removes sediment by reducing runoff velocity due to increased resistance of the vegetative media and enhanced infiltration in the VFS area (Barfield et al., 1998). Both result in a reduction in transport capacity and the deposition of sediment. Both the filtering and infiltration aspects are represented in the model. Similar to the runoff loading approach used earlier, sediment loading per unit VFS area was found to correlate with measured sediment reduction. Dosskey et al. (2002) hypothesized that sediment trapping efficiency decreases as the load per unit of buffer area increases. Sediment loading was calculated as the mass of sediment originating from the upslope area per unit of VFS area express as kg m$^{-2}$. The infiltration aspect was represented in the model by incorporating the runoff reduction as a percentage. Sixty-two experiments reported in the literature were used to develop this model.

$S_R(\%)=79.0-1.04 S_L+0.213R_R$ 6:5.1.2

where $S_R$ is the predicted sediment reduction (%); $S_L$ is sediment loading (kg/m$^2$); and $R_R$ is the runoff reduction (%). Sediment loading alone was correlated with sediment reduction ($R^2$ = 0.41) (Figure 6:5-2). The addition of runoff reduction allowed the regression model to explain most of the variability ($R^2$ = 0.64) in the measured data.

The total nitrogen model was based on sediment reduction only. Much of the nitrogen lost in runoff from agricultural fields travels with sediments. Harmel et al. (2006) found that approximately 75% of the nitrogen lost from conventional tilled fields was in particulate forms. They also found that dissolved nitrogen forms, such as nitrate, were more dominant in no-till treatments. The vast majority of VFS data derived from literature were designed to simulate higher erosion conditions where particulate forms would represent the majority of nitrogen losses.

The total nitrogen model was based on sediment reduction from 44 observations reported in the literature. Two trials were censored during the development of the model. These experiments from Magette et al. (1989) yielded significant increases in total nitrogen exiting the VFS. The authors attributed this phenomenon to flushing of fine particulates captured in the VFS from prior experimental trials. Both the slope and the intercept were significant (P < 0.01). The model is given below and shown in Figure 6:5-3.

$TN_R=0.036 S_R^{1.69}$ 6:5.1.3

where $TN_R$ is the total nitrogen reduction (%); and $S_R$ is the sediment reduction (%). Although this model was developed from total nitrogen, which includes both soluble and particulate forms, it was applied only to particulate forms in the SWAT+ model.

The nitrate nitrogen model was developed from 42 observations. Four observations from Dillaha et al. (1989) had negative runoff reduction values due to additional runoff generated in the VFS. Because the VFS SWAT+ sub model is not allowed to generate additional loads, these observations were censored. All nutrient models initially included both runoff and sediment reductions as independent variables, but the nitrate nitrogen model was the only model where both were significant ($\alpha$=0.05). Nitrate is soluble and should not be associated with sediments, yet they were statistically correlated in the measured data. It is likely that the relationship between nitrate and sediment is an artifact of cross- correlation between sediment and runoff reductions (as demonstrated by Equation (2)). The nitrate model was based only on runoff reduction; both the slope and the intercept were significant (p<0.01). The nitrate nitrogen model is given below:

$NN_R=39.4+0.584R_R$ 6:5.1.4

where $NN_R$ is the nitrate nitrogen reduction (%); and $R_R$ is the runoff reduction (%). Because both the slope and the intercept were significant, there is a minimum reduction of 39.4% in nitrate, even if there is no reduction in runoff due to the VFS. This outcome may be unexpected, but it is supported by the measured data. Dillaha et al. (1989) observed nitrate reductions of 52% and 32% with only 0% and 7% reductions in runoff volume. Lee et al. (2000) also found significant reductions in nitrate (61%) with low runoff reductions (7%). One possible explanation is that sufficient runoff can be generated in the VFS such that there is little net reduction in runoff, but significant infiltration may still occur. Another possibility is foliar uptake of nitrates by vegetation within the strip.

The model for total phosphorus was based on sediment reduction. Although total phosphorus is composed of both soluble and particulate forms, particulate forms represent the bulk of phosphorus lost from conventionally cultivated fields. The total phosphorus model was developed from 63 observations; more data than any other nutrient model. The intercept was not significant. Sediment reduction was able to explain 43% of the variability. The model was applied to all particulate forms of phosphorus in the SWAT+ VFS submodel. The model is given below:

$TP_R=0.90 S_R$ 6:5.1.5

where $TP_R$ is the total phosphorus reduction (%); and $S_R$ is the sediment reduction (%).

The soluble phosphorus model was based on runoff reduction only. The observations censored from (Dillaha et al., 1989) for the nitrates were also censored for soluble phosphorus. The soluble phosphorus model has the weakest relationship of all the nutrient models (R$^2$ = 0.27), yet both the slope and intercept were significant (p = 0.01).

$DP_R=29.3+0.51R_R$ 6:5.1.6

where $DP_R$ is the dissolved phosphorus reduction (%); and $R_R$ is the runoff reduction (%). Like nitrate, there is a significant reduction in soluble phosphorus (29.3%) with no corresponding runoff reduction. Experimental observations of soluble phosphorus reduction at low runoff reductions are highly variable. Dillaha et al. (1989) found reductions in soluble phosphorus ranging from 43% to -31% with near zero runoff reduction. The minimum reduction predicted by equation 6 could be the result of mechanisms similar to those cited for the removal of nitrates, or simply an artifact of experimental variability.

Support for grass waterways was added to SWAT+. Waterways are treated as trapezoidal channels; the deposition of sediment and organic nutrients is calculated in the same manner as SWAT+ subbasin tributary channels. The primary user inputs are waterway width and length.

The sediment transport capacity is defined as:

$Scap=Spcon*v^{1.5}$ 6:5.2.1

where $Scap$ is the sediment transport capacity in (mg/m$^3$) , $Spcon$ is the sediment transport coefficient and $v$ is flow velocity in the waterway (m/s).

Unsubmerged portions of the waterway act as filter strips, and may trap both soluble and organic nutrients. These equations are simplified forms of those used by White and Arnold (2009) in the simulation of filter strips. Removal of soluble pollutants from the unsubmerged portion is calculates as:

$SolR=75.8-10.8log(SD)+25.9log(SolK)$ 6:5.2.2

where $SolR$ is soluble pollutant removal (%), $SD$ is runoff depth over unsubmerged waterway area in (mm/day), and $SolK$ is the saturated hydrologic conductivity of the soil surface (mm/hr). Removal of particulate pollutant and sediment in the unsubmerged area is calculated as:

$SedR=79.0-1.04(SedL)+0.213*(SolR)$ 6:5.2.3

where $SedR$ is the sediment and particulate pollutant removal (%) and $SedL$ is the sediment load per unit area of unsubmerged waterway in (kg/ha/day).

VFSs were implemented at the HRU level in SWAT+. Three additional model parameters were added as SWAT+ inputs: the drainage area to VFS area ratio ($DAFS_{ratio}$), the fraction of the field drained by the most heavily loaded 10% of the VFS ($DF_{con}$), and the fraction of the flow through the most heavily loaded 10% of the VFS which is fully channelized ($CF_{frac}$), all are specified in the HRU (.hru) file. A two-segment VFS was used. Section one represents the bulk of the VFS area (90%) which receives the least flow. Section two is the remaining 10% of the buffer which receives between 25% and 75% of the field runoff (Figure 6:5-4).

The fraction of flow through section two which is channelized is not subject to the VFS model; all sediment and nutrient are conservatively delivered to the tributary channel. $DAFS_{ratio}$ for sections one and two are calculated from $DF_{con}$ using the following equations:

$DAFS_{ratio1}=DAFS_{ratio}(1-DF_{con})/0.9$ 6:5.1.7

$DAFS_{ratio2}=DAFS_{ratio}(1-CF_{frac})DF_{con}/0.1$ 6:5.1.8

where $DAFS_{ratio1}$ is the drainage area to VFS area ratio for section 1; $DAFS_{ratio2}$ is the drainage area to the VFS area ratio for section 2; $DAFS_{ratio}$ is the average drainage area to the VFS area ratio for the entire HRU (user input).; $CF_{frac}$ is fraction of the flow through the most heavily loaded 10% of the VFS which is fully channelized (user input); and $DF_{con}$ is the fraction of the field drained by the most heavily loaded 10% of the VFS (user input). Sediment, runoff, and nutrient loadings are calculated assuming all are generated uniformly across the HRU. The $DAFS_{ratio}$ for each VFS section is combined with SWAT+ HRU level runoff and sediment yield predictions to calculate the runoff and sediment loadings. Equations 6:5.1-6 are applied to predict sediment and nutrient transport through the VFS. The fraction of runoff retained in VFS is calculated for the purposes of estimating the retention of other constituents only. It is beyond the scope of this research to predict the aspects of a VFS’s hydrologic budget needed to represent that component in the SWAT+ model. In addition, the area occupied by VFSs within a SWAT+ HRU is not removed from that HRU for simplicity. For these reasons, the VFS routine in SWAT+ is not used to predict changes in runoff delivered to streams. The VFS SWAT+ sub model also includes pesticides and bacteria. Due to a lack of measured data, these models are based on assumptions. The pesticide model assumes that pesticides sorbed to sediments are captured with the sediment, and soluble pesticides are captured with runoff. Similarly, bacteria, which are attached to sediment, are captured with sediment and unattached bacteria are captured with runoff. These assumptions are common in the structure of other SWAT+ model components.

In addition to the mechanisms by which sediment and runoff are captured, nutrients may be adsorbed onto vegetation, surface residues, or the soil surface (Barfield et al., 1998). For the sake of simplicity, nutrient reduction was considered to be a function of sediment or runoff reduction only. Only nitrogen and phosphorus were considered. All nutrient models were developed from measured VFS data; the current version of VFSMOD does not account for nutrients.

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).