Surface runoff occurs whenever the rate of water application to the ground surface exceeds the rate of infiltration. When water is initially applied to a dry soil, the infiltration rate is usually very high. However, it will decrease as the soil becomes wetter. When the application rate is higher than the infiltration rate, surface depressions begin to fill. If the application rate continues to be higher than the infiltration rate once all surface depressions have filled, surface runoff will commence.

SWAT+ provides two methods for estimating surface runoff: the SCS curve number procedure (SCS, 1972) and the Green & Ampt infiltration method (1911).

The SCS curve number is a function of the soil’s permeability, land use and antecedent soil water conditions. Typical curve numbers for moisture condition II are listed in tables 2:1-1, 2:1-2 and 2:1-3 for various land covers and soil types (SCS Engineering Division, 1986). These values are appropriate for a 5% slope.

Table 2:1-1: Runoff curve numbers for cultivated agricultural lands (from SCS Engineering Division, 1986)

[1] Poor: < 50% ground cover or heavily grazed with no mulch; Fair: 50 to 75% ground cover and not heavily grazed; Good: > 75% ground cover and lightly or only occasionally grazed

[2] Poor: < 50% ground cover; Fair: 50 to 75% ground cover; Good: > 75% ground cover

[3] Poor: Forest litter, small trees, and brush are destroyed by heavy grazing or regular burning; Fair: Woods are grazed but not burned, and some forest litter covers the soil; Good: Woods are protected from grazing, and litter and brush adequately cover the soil.

[1] Poor: < 50% ground cover or heavily grazed with no mulch; Fair: 50 to 75% ground cover and not heavily grazed; Good: > 75% ground cover and lightly or only occasionally grazed

[2] Poor: < 50% ground cover; Fair: 50 to 75% ground cover; Good: > 75% ground cover

[3] Poor: Forest litter, small trees, and brush are destroyed by heavy grazing or regular burning; Fair: Woods are grazed but not burned, and some forest litter covers the soil; Good: Woods are protected from grazing, and litter and brush adequately cover the soil.

[2] Poor: < 50% ground cover or heavily grazed with no mulch;

[3] Fair: 50 to 75% ground cover and not heavily grazed;

[4] Good: > 75% ground cover and lightly or only occasionally grazed

The land phase of the hydrologic cycle is based on the water balance equation:

The SCS runoff equation is an empirical model that came into common use in the 1950s. It was the product of more than 20 years of studies involving rainfall-runoff relationships from small rural watersheds across the U.S. The model was developed to provide a consistent basis for estimating the amounts of runoff under varying land use and soil types (Rallison and Miller, 1981).

The SCS curve number equation is (SCS, 1972):

$Q_{surf}=\frac{(R_{day}-I_a)^2}{(R_{day}-I_a+S)}$ 2:1.1.1

where $Q_{surf}$ is the accumulated runoff or rainfall excess (mm H$_2$O), $R_{day}$ is the rainfall depth for the day (mm H$_2$O), $I_a$ is the initial abstractions which includes surface storage, interception and infiltration prior to runoff (mm H$_2$O), and $S$ is the retention parameter (mm H$_2$O). The retention parameter varies spatially due to changes in soils, land use, management and slope and temporally due to changes in soil water content. The retention parameter is defined as:

$S=25.4(\frac{1000}{CN}-10)$ 2:1.1.2

where $CN$ is the curve number for the day. The initial abstractions, $I_a$, is commonly approximated as $0.2S$ and equation 2:1.1.1 becomes

$Q_{surf}=\frac{(R_{day}-0.2S)^2}{(R_{day}+0.8S)}$ 2:1.1.3

Runoff will only occur when $R_{day}> I_a$. A graphical solution of equation 2:1.1.3 for different curve number values is presented in Figure 2:1-1.

The U.S. Natural Resource Conservation Service (NRCS) classifies soils into four hydrologic groups based on infiltration characteristics of the soils. NRCS Soil Survey Staff (1996) defines a hydrologic group as a group of soils having similar runoff potential under similar storm and cover conditions. Soil properties that influence runoff potential are those that impact the minimum rate of infiltration for a bare soil after prolonged wetting and when not frozen.

These properties are depth to seasonally high water table, saturated hydraulic conductivity, and depth to a very slowly permeable layer.

Soil may be placed in one of four groups, A, B, C, and D, or three dual classes, A/D, B/D, and C/D.

Definitions of the classes are:

A: (Low runoff potential). The soils have a high infiltration rate even when thoroughly wetted. They chiefly consist of deep, well drained to excessively drained sands or gravels. They have a high rate of water transmission.

B: The soils have a moderate infiltration rate when thoroughly wetted. They chiefly are moderately deep to deep, moderately well-drained to well-drained soils that have moderately fine to moderately coarse textures. They have a moderate rate of water transmission.

C: The soils have a slow infiltration rate when thoroughly wetted. They chiefly have a layer that impedes downward movement of water or have moderately fine to fine texture. They have a slow rate of water transmission.

D. (High runoff potential). The soils have a very slow infiltration rate when thoroughly wetted. They chiefly consist of clay soils that have a high swelling potential, soils that have a permanent water table, soils that have a claypan or clay layer at or near the surface, and shallow soils over nearly impervious material. They have a very slow rate of water transmission.

Dual hydrologic groups are given for certain wet soils that can be adequately drained. The first letter applies to the drained condition, the second to the undrained. Only soils that are rated D in their natural condition are assigned to dual classes.

SCS defines three antecedent moisture conditions:

I—dry (wilting point), II—average moisture, and III—wet (field capacity). The moisture condition I curve number is the lowest value the daily curve number can assume in dry conditions. The curve numbers for moisture conditions I and III are calculated with the equations:

$CN_1=CN_2-\frac{20*(100-CN_2)}{(100-CN_2+exp[2.533-0.0636*(100-CN_2)])}$ 2:1.1.4

$CN_3=CN_2*exp[0.00673*(100-CN_2)]$ 2:1.1.5

where $CN_1$ is the moisture condition I curve number, $CN_2$ is the moisture condition II curve number, and $CN_3$ is the moisture condition III curve number.

With SWAT+, users are allowed to select between two methods for calculating the retention parameter. The traditional method is to allow the retention parameter to vary with soil profile water content. An alternative added in SWAT+ allows the retention parameter to vary with accumulated plant evapotranspiration. Calculation of the daily CN value as a function of plant evapotranspiration was added because the soil moisture method was predicting too much runoff in shallow soils. By calculating daily CN as a function of plant evapotranspiration, the value is less dependent on soil storage and more dependent on antecedent climate.

When the retention parameter varies with soil profile water content, the following equation is used:

$S=S_{max}*(1-\frac{SW}{[SW+exp(w_1-w_2*SW)]})$ 2:1.1.6

where $S$ is the retention parameter for a given day (mm), $S_{max}$ is the maximum value the retention parameter can achieve on any given day (mm), $SW$ is the soil water content of the entire profile excluding the amount of water held in the profile at wilting point (mm H$_2$O), and $w_1$ and $w_2$ are shape coefficients. The maximum retention parameter value, $S_{max}$, is calculated by solving equation 2:1.1.2 using $CN_1$.

The shape coefficients are determined by solving equation 2:1.1.6 assuming that

1. the retention parameter for moisture condition I curve number corresponds to wilting point soil profile water content,

2. the retention parameter for moisture condition III curve number corresponds to field capacity soil profile water content, and

3. the soil has a curve number of 99 (S = 2.54) when completely saturated.

$w_1=ln[\frac{FC}{1-S_3*S_{max}^{-1}}-FC]+w_2*FC$ 2.1.1.7

$w_2=\frac{(ln[\frac{FC}{1-S_3*S_{max}^{-1}}-FC]-ln[\frac{SAT}{1-2.54*S_{max}^{-1}}-SAT])}{(SAT-FC)}$ 2.1.1.8

where $w_1$ is the first shape coefficient, $w_2$ is the second shape coefficient, $FC$ is the amount of water in the soil profile at field capacity (mm H$_2$O), $S_3$ is the retention parameter for the moisture condition III curve number, $S_{max}$ is the retention parameter for the moisture condition I curve number, $SAT$ is the amount of water in the soil profile when completely saturated (mm H$_2$O), and 2.54 is the retention parameter value for a curve number of 99.

When the retention parameter varies with plant evapotranspiration, the following equation is used to update the retention parameter at the end of every day:

$S=S_{prev}+E_o *exp(\frac{-cncoef*S_{prev}}{S_{max}})-R_{day}+Q_{surf}$ 2:1.1.9

where $S$ is the retention parameter for a given day (mm), $S_{prev}$ is the retention parameter for the previous day (mm), $E_o$ is the potential evapotranspiration for the day (mm d$^{-1}$), $cncoef$ is the weighting coefficient used to calculate the retention coefficient for daily curve number calculations dependent on plant evapotranspiration, $S_{max}$ is the maximum value the retention parameter can achieve on any given day (mm), Rday is the rainfall depth for the day (mm H$_2$O), and $Q_{surf}$ is the surface runoff (mm H$_2$O). The initial value of the retention parameter is defined as

$S=0.9*S_{max}$

When the top layer of the soil is frozen, the retention parameter is modified using the following equation:

$S_{frz}=S_{max}*[1-exp(-0.000862*S)]$ 2:1.1.10

where $S_{frz}$ is the retention parameter adjusted for frozen conditions (mm), $S_{max}$ is the maximum value the retention parameter can achieve on any given day (mm), and $S$ is the retention parameter for a given moisture content calculated with equation 2:1.1.6 (mm).

The daily curve number value adjusted for moisture content is calculated by rearranging equation 2:1.1.2 and inserting the retention parameter calculated for that moisture content:

$CN=\frac{25400}{(S+254)}$ 2:1.1.11

where $CN$ is the curve number on a given day and $S$ is the retention parameter calculated for the moisture content of the soil on that day.

The time of concentration is the amount of time from the beginning of a rainfall event until the entire subbasin area is contributing to flow at the outlet. In other words, the time of concentration is the time for a drop of water to flow from the remotest point in the subbasin to the subbasin outlet. The time of concentration is calculated by summing the overland flow time (the time it takes for flow from the remotest point in the subbasin to reach the channel) and the channel flow time (the time it takes for flow in the upstream channels to reach the outlet):

$t_{conc}=t_{ov}+t_{ch}$ 2:1.3.2

where $t_{conc}$ is the time of concentration for a subbasin (hr), $t_{ov}$ is the time of concentration for overland flow (hr), and $t_{ch}$ is the time of concentration for channel flow (hr).

The peak runoff rate is the maximum runoff flow rate that occurs with a given rainfall event. The peak runoff rate is an indicator of the erosive power of a storm and is used to predict sediment loss. SWAT+ calculates the peak runoff rate with a modified rational method.

The rational method is widely used in the design of ditches, channels and storm water control systems. The rational method is based on the assumption that if a rainfall of intensity $i$ begins at time $t=0$ and continues indefinitely, the rate of runoff will increase until the time of concentration, $t=t_{conc}$, when the entire subbasin area* is contributing to flow at the outlet. The rational formula is:

$q_{peak}=\frac{C*i*Area}{3.6}$ 2:1.3.1

where $q_{peak}$ is the peak runoff rate ($m^3 s^{-1}$), $C$ is the runoff coefficient, $i$ is the rainfall intensity (mm/hr), Area is the subbasin area (km$^2$) and 3.6 is a unit conversion factor.

The moisture condition II curve numbers provided in the tables are assumed to be appropriate for 5% slopes. Williams (1995) developed an equation to adjust the curve number to a different slope:

$CN_{2s}=\frac{(CN_3-CN_2)}{3}*[1-2*exp(-13.86*slp)]+CN_2$ 2:1.1.12

where $CN_{2s}$ is the moisture condition II curve number adjusted for slope, $CN_3$ is the moisture condition III curve number for the default 5% slope, $CN_2$ is the moisture condition II curve number for the default 5% slope, and $slp$ is the average fraction slope of the subbasin. SWAT+ does not adjust curve numbers for slope. If the user wishes to adjust the curve numbers for slope effects, the adjustment must be done prior to entering the curve numbers in the management input file.

Table 2:1-1: SWAT+ input variables that pertain to surface runoff calculated with the SCS curve number method.

The overland flow time of concentration, tov, can be computed using the equation

$t_{ov}=\frac{L_{slp}}{3600*v_{ov}}$ 2:1.3.3

where $L_{slp}$ is the subbasin slope length (m), $v_{ov}$ is the overland flow velocity (m s$^{-1}$) and 3600 is a unit conversion factor.

The overland flow velocity can be estimated from Manning’s equation by considering a strip 1 meter wide down the sloping surface:

$v_{ov}=\frac{q_{ov}^{0.4}*slp^{0.3}}{n^{0.6}}$ 2:1.3.4

where $q_{ov}$ is the average overland flow rate ($m^3 s^{-1}$), $slp$ is the average slope in the subbasin (m $m^{-1}$), and $n$ is Manning’s roughness coefficient for the subbasin. Assuming an average flow rate of 6.35 mm/hr and converting units

$v_{ov}=\frac{0.005*L_{slp}^{0.4}*slp^{0.3}}{n^{0.6}}$ 2:1.3.5

Substituting equation 2:1.3.5 into equation 2:1.3.3 gives

$t_{ov}=\frac{L_{slp}^{0.6}*n^{0.6}}{18*slp^{0.3}}$ 2:1.3.6

The Green & Ampt equation was developed to predict infiltration assuming excess water at the surface at all times (Green & Ampt, 1911). The equation assumes that the soil profile is homogenous and antecedent moisture is uniformly distributed in the profile. As water infiltrates into the soil, the model assumes the soil above the wetting front is completely saturated and there is a sharp break in moisture content at the wetting front. Figure 2:1-2 graphically illustrates the difference between the moisture distribution with depth modeled by the Green & Ampt equation and what occurs in reality.

Mein and Larson (1973) developed a methodology for determining ponding time with infiltration using the Green & Ampt equation. The Green-Ampt Mein-Larson excess rainfall method was incorporated into SWAT+ to provide an alternative option for determining surface runoff. This method requires sub-daily precipitation data supplied by the user.

The Green-Ampt Mein-Larson infiltration rate is defined as:

$f_{inf,t}=K_e*(1+\frac{\Psi_{wf}*\Delta\theta_v}{F_{inf,t}})$ 2:1.2.1

where $f_{inf}$ is the infiltration rate at time $t$ (mm/hr), $K_e$ is the effective hydraulic conductivity (mm/hr), $\Psi_{wf}$ is the wetting front matric potential (mm), $\Delta\theta_v$ is the change in volumetric moisture content across the wetting front (mm/mm) and $F_{inf}$ is the cumulative infiltration at time $t$ (mm H$_2$O).

When the rainfall intensity is less than the infiltration rate, all the rainfall will infiltrate during the time period and the cumulative infiltration for that time period is calculated:

$F_{inf,t}=F_{inf,t-1}+R_{\Delta t}$ 2:1.2.2

where $F_{inf,t}$ is the cumulative infiltration for a given time step (mm H$_2$O), $F_{inf,t-1}$ is the cumulative infiltration for the previous time step(mm H$_2$O), and $R_{\Delta t}$ is the amount of rain falling during the time step (mm H$_2$O).

The infiltration rate defined by equation 2:1.2.1 is a function of the infiltrated volume, which in turn is a function of the infiltration rates in previous time steps. To avoid numerical errors over long time steps, $f_{inf}$ is replaced by $dF_{inf}/dt$ in equation 2:1.2.1 and integrated to obtain

$F_{inf,t}=F_{inf,t-1}+K_e*\Delta t+ \Psi_{wf}*\Delta\theta_v*ln[\frac{F_{inf,t}+\Psi_{wf}*\Delta\theta_v}{F_{inf,t-1}+\Psi_{wf}*\Delta\theta_v}]$ 2:1.2.3

Equation 2:1.2.3 must be solved iteratively for $F_{inf,t}$, the cumulative infiltration at the end of the time step. A successive substitution technique is used.

The Green-Ampt effective hydraulic conductivity parameter, $K_e$, is approximately equivalent to one-half the saturated hydraulic conductivity of the soil, $K_{sat}$ (Bouwer, 1969). Nearing et al. (1996) developed an equation to calculate the effective hydraulic conductivity as a function of saturated hydraulic conductivity and curve number. This equation incorporates land cover impacts into the calculated effective hydraulic conductivity. The equation for effective hydraulic conductivity is:

$K_e=\frac{56.82*K_{sat}^{0.286}}{1+0.051*exp(0.062*CN)}-2$ 2:1.2.4

where $K_e$ is the effective hydraulic conductivity (mm/hr), $K_{sat}$ is the saturated hydraulic conductivity (mm/hr), and $CN$ is the curve number.

Wetting front matric potential, $\Psi_{wf}$, is calculated as a function of porosity, percent sand and percent clay (Rawls and Brakensiek, 1985):

$\Psi_{wf}=10*exp[6.5309-7.32561*\Phi_{soil}+0.001583*m_c^2+3.809479*\Phi_{soil}^2+0.000344*m_s*m_c-0.049837*m_s*\Phi_{soil}+0.001608*m_s^2*\Phi_{soil}^2+0.001602*m_c^2*\Phi_{soil}^2-0.0000136*m_s^2*m_c-0.003479*m_c^2*\Phi_{soil}-0.000799*m_s^2*\Phi_{soil}]$

2:1.2.5

where $\Phi_{soil}$ is the porosity of the soil (mm/mm), $m_c$ is the percent clay content, and $m_s$ is the percent sand content.

The change in volumetric moisture content across the wetting front is calculated at the beginning of each day as:

$\Delta\theta_v=(1-\frac{SW}{FC})*(0.95*\Phi_{soil})$ 2:1.2.6

where $\Delta\theta_v$ is the change in volumetric moisture content across the wetting front (mm/mm), $SW$ is the soil water content of the entire profile excluding the amount of water held in the profile at wilting point (mm H$_2$O), $FC$ is the amount of water in the soil profile at field capacity (mm H$_2$O), and $\Phi_{soil}$ is the porosity of the soil (mm/mm). If a rainfall event is in progress at midnight, $\Delta\theta_v$ is then calculated:

$\Delta\theta_v=0.001*(0.95*\Phi_{soil})$ 2:1.2.7

For each time step, SWAT+ calculates the amount of water entering the soil. The water that does not infiltrate into the soil becomes surface runoff.

Table 2:1-2: SWAT+ input variables that pertain to Green & Ampt infiltration calculations.

The channel flow time of concentration, , can be computed using the equation

2:1.3.7

where is the average flow channel length for the subbasin (km), is the average channel velocity (m s), and 3.6 is a unit conversion factor.

The average channel flow length can be estimated using the equation

2:1.3.8

where is the channel length from the most distant point to the subbasin outlet (km), and is the distance along the channel to the subbasin centroid (km). Assuming , the average channel flow length is

2:1.3.9

The average velocity can be estimated from Manning’s equation assuming a trapezoidal channel with 2:1 side slopes and a 10:1 bottom width-depth ratio.

2:1.3.10

where is the average channel velocity (m s), is the average channel flow rate (), is the channel slope (m m), and is Manning’s roughness coefficient for the channel. To express the average channel flow rate in units of mm/hr, the following expression is used

2.1.3.11

where is the average channel flow rate (mm hr), is the subbasin area (km), and 3.6 is a unit conversion factor. The average channel flow rate is related to the unit source area flow rate (unit source area = 1 ha)

2:1.3.12

where is the unit source area flow rate (mm hr), is the subbasin area (km), and 100 is a unit conversion factor. Assuming the unit source area flow rate is 6.35 mm/hr and substituting equations 2:1.3.11 and 2:1.3.12 into 2:1.3.10 gives

2:1.3.13

Substituting equations 2:1.3.9 and 2:1.3.13 into 2:1.3.7 gives

2:1.3.14

Although some of the assumptions used in developing equations 2:1.3.6 and 2:1.3.14 may appear liberal, the time of concentration values obtained generally give satisfactory results for homogeneous subbasins. Since equations 2:1.3.6 and 2:1.3.14 are based on hydraulic considerations, they are more reliable than purely empirical equations.

The runoff coefficient is the ratio of the inflow rate, $i*Area$, to the peak discharge rate, $q_{peak}$. The coefficient will vary from storm to storm and is calculated with the equation:

$C=\frac{Q_{surf}}{R_{day}}$ 2:1.3.15

where $Q_{surf}$ is the surface runoff (mm H$_2$O) and $R_{day}$ is the rainfall for the day (mm H$_2$O).

Evapotranspiration is a collective term that includes all processes by which water at the earth’s surface is converted to water vapor. It includes evaporation from the plant canopy, transpiration, sublimation and evaporation from the soil.

Evapotranspiration is the primary mechanism by which water is removed from a watershed. Roughly 62% of the precipitation that falls on the continents is evapotranspired. Evapotranspiration exceeds runoff in most river basins and on all continents except Antarctica (Dingman, 1994).

The difference between precipitation and evapotranspiration is the water available for human use and management. An accurate estimation of evapotranspiration is critical in the assessment of water resources and the impact of climate and land use change on those resources.

The rainfall intensity is the average rainfall rate during the time of concentration. Based on this definition, it can be calculated with the equation:

$i=\frac{R_{tc}}{t_{conc}}$ 2:1.3.16

where $i$ is the rainfall intensity (mm/hr), $R_{tc}$ is the amount of rain falling during the time of concentration (mm H$_2$O), and $t_{conc}$ is the time of concentration for the subbasin (hr).

An analysis of rainfall data collected by Hershfield (1961) for different durations and frequencies showed that the amount of rain falling during the time of concentration was proportional to the amount of rain falling during the 24-hr period.

$R_{tc}=\alpha_{tc}*R_{day}$ 2:1.3.17

where $R_{tc}$ is the amount of rain falling during the time of concentration (mm H$_2$O), $\alpha_{tc}$ is the fraction of daily rainfall that occurs during the time of concentration, and $R_{day}$ is the amount of rain falling during the day (mm H$_2$O).

For short duration storms, all or most of the rain will fall during the time of concentration, causing $\alpha_{tc}$ to approach its upper limit of 1.0. The minimum value of $\alpha_{tc}$ would be seen in storms of uniform intensity ($i_{24}=i$). This minimum value can be defined by substituting the products of time and rainfall intensity into equation 2:1.3.17

$\alpha_{tc,min}=\frac{R_{tc}}{R_{day}}=\frac{i*t_{conc}}{i_{24}*24}=\frac{t_{conc}}{24}$ 2:1.3.18

Thus, $\alpha_{tc}$ falls in the range $t_{conc}/24 \le \alpha_{tc} \le1.0$

SWAT+ estimates the fraction of rain falling in the time of concentration as a function of the fraction of daily rain falling in the half-hour of highest intensity rainfall.

$\alpha_{tc}=1-exp[2*t_{conc}*ln(1-\alpha_{0.5})]$ 2:1.3.19

where $\alpha_{0.5}$ is the fraction of daily rain falling in the half-hour highest intensity rainfall, and $t_{conc}$ is the time of concentration for the subbasin (hr). The determination of a value for $\alpha_{0.5}$ is discussed in Chapters 1:2 and 1:3.

Many semiarid and arid watersheds have ephemeral channels that abstract large quantities of streamflow (Lane, 1982). The abstractions, or transmission losses, reduce runoff volume as the flood wave travels downstream. Chapter 19 of the SCS Hydrology Handbook (Lane, 1983) describes a procedure for estimating transmission losses for ephemeral streams which has been incorporated into SWAT+. This method was developed to estimate transmission losses in the absence of observed inflow-outflow data and assumes no lateral inflow or out-of-bank flow contributions to runoff.

The prediction equation for runoff volume after transmission losses is

$vol_{Qsurf,f}=\begin {cases} 0 & vol{Qsurf,i} \le vol_{thr} \\ a_x+b_x*vol_{Qsurf,i} & vol_{Qsurf,i} > vol_{thr} \end{cases}$ 2:1.5.1

where $vol_{Qsurf,f}$ is the volume of runoff after transmission losses ($m^3$), $a_x$ is the regression intercept for a channel of length $L$ and width $W$ ($m^3$), $b_x$ is the regression slope for a channel of length $L$ and width $W$, $vol_{Qsurf,i}$,$i$ is the volume of runoff prior to transmission losses ($m^3$), and $vol_{thr}$ is the threshold volume for a channel of length $L$ and width $W$ ($m^3$). The threshold volume is

$vol_{thr}=-\frac{a_x}{b_x}$ 2:1.5.2

The corresponding equation for peak runoff rate is

$q_{peak,f}=\frac{1}{(3600*dur_{flw})}*[a_x-(1-b_x)*vol_{Qsurf,i}]+b_x*q_{peak,i}$ 2:1.5.3

where $q_{peak,f}$ is the peak rate after transmission losses ($m^3$/s), $dur_{flw}$ is the duration of flow (hr), $a_x$ is the regression intercept for a channel of length $L$ and width $W$ ($m^3$), $b_x$ is the regression slope for a channel of length $L$ and width $W$, $vol_{Qsurf,i}$ is the volume of runoff prior to transmission losses ($m^3$), $q_{peak,i}$ is the peak rate before accounting for transmission losses ($m^3$/s). The duration of flow is calculated with the equation:

$dur_{flw}=\frac{Q_{surf}*Area}{3.6*q_{peak}}$ 2:1.5.4

where $dur_{flw}$ is the duration of runoff flow (hr),$Q_{surf}$ is the surface runoff (mm H$_2$O), $Area$ is the area of the subbasin (km$^2$), $q_{peak}$ is the peak runoff rate (m$^3$/s), and 3.6 is a conversion factor.

In order to calculate the regression parameters for channels of differing lengths and widths, the parameters of a unit channel are needed. A unit channel is defined as a channel of length $L$= 1 km and width $W$= 1 m. The unit channel parameters are calculated with the equations:

$k_r=-2.22*ln[1-2.6466*\frac{K_{ch}*dur_{flw}}{vol_{Qsurf,i}}]$ 2:1.5.5

$a_r=-0.2258*K_{ch}*dur_{flw}$ 2:1.5.6

$b_r=exp[-0.4905*k_r]$ 2:1.5.7

where $k_r$ is the decay factor ($m^{-1}$ k$m^{-1}$), $a_r$ is the unit channel regression intercept ($m^3$), $b_r$ is the unit channel regression slope, $K_{ch}$ is the effective hydraulic conductivity of the channel alluvium (mm/hr), $dur_{flw}$ is the duration of runoff flow (hr), and $vol_{Qsurf,i}$ is the initial volume of runoff ($m^3$). The regression parameters are

$b_x=exp[-k_r*L*W]$ 2:1.5.8

$a_x=\frac{a_r}{(1-b_r)}*(1-b_x)$ 2:1.5.9

where $a_x$ is the regression intercept for a channel of length $L$ and width $W$ ($m^3$), $b_x$ is the regression slope for a channel of length $L$ and width$W$,$k_r$ is the decay factor ($m^{-1}$ k$m^{-1}$), $L$ is the channel length from the most distant point to the subbasin outlet (km), $W$ is the average width of flow, i.e. channel width (m) $a_r$ is the unit channel regression intercept ($m^3$), and $b_r$ is the unit channel regression slope.

Transmission losses from surface runoff are assumed to percolate into the shallow aquifer.

Table 2:1-7: SWAT+ input variables that pertain to transmission loss calculations.

In large subbasins with a time of concentration greater than 1 day, only a portion of the surface runoff will reach the main channel on the day it is generated. SWAT+ incorporates a surface runoff storage feature to lag a portion of the surface runoff release to the main channel.

Once surface runoff is calculated with the Curve Number or Green & Ampt method, the amount of surface runoff released to the main channel is calculated:

$Q_{surf}=(Q'_{surf}+Q_{stor,i-1})*(1-exp[\frac{-surlag}{t_{conc}}])$ 2:1.4.1

where $Q_{surf}$ is the amount of surface runoff discharged to the main channel on a given day (mm H$_2$O), is the amount of surface runoff generated in the subbasin on a given day (mm H$_2$O), $Q_{stor,i-1}$ is the surface runoff stored or lagged from the previous day (mm H$_2$O), $surlag$ is the surface runoff lag coefficient, and $t_{conc}$ is the time of concentration for the subbasin (hrs).

The expression$(1-exp[\frac{-surlag}{t_{conc}}])$ in equation 2:1.4.1 represents the fraction of the total available water that will be allowed to enter the reach on any one day. Figure 2:1-3 plots values for this expression at different values for $surlag$and $t_{conc}$.

Figure 2:1-3: Influence of $surlag$ and $t_{conc}$ on fraction of surface runoff released.

Note that for a given time of concentration, as $surlag$ decreases in value more water is held in storage. The delay in release of surface runoff will smooth the streamflow hydrograph simulated in the reach.

Table 2:1-6: SWAT+ input variables that pertain to surface runoff lag calculations.

The modified rational formula used to estimate peak flow rate is obtained by substituting equations 2:1.3.15, 2:1.3.16, and 2:1.3.17 into equation 2:1.3.1

2:1.3.20

where is the peak runoff rate (), is the fraction of daily rainfall that occurs during the time of concentration, is the surface runoff (mm HO), is the subbasin area (km), is the time of concentration for the subbasin (hr) and 3.6 is a unit conversion factor.

Table 2:1-5: SWAT+ input variables that pertain to peak rate calculations.

The plant canopy can significantly affect infiltration, surface runoff and evapotranspiration. As rain falls, canopy interception reduces the erosive energy of droplets and traps a portion of the rainfall within the canopy. The influence the canopy exerts on these processes is a function of the density of plant cover and the morphology of the plant species.

When calculating surface runoff, the SCS curve number method lumps canopy interception in the term for initial abstractions. This variable also includes surface storage and infiltration prior to runoff and is estimated as 20% of the retention parameter value for a given day (see Chapter 2:1). When the Green and Ampt infiltration equation is used to calculate surface runoff and infiltration, the interception of rainfall by the canopy must be calculated separately.

SWAT+ allows the maximum amount of water that can be held in canopy storage to vary from day to day as a function of the leaf area index:

2:2.1.1

where is the maximum amount of water that can be trapped in the canopy on a given day (mm HO), is the maximum amount of water that can be trapped in the canopy when the canopy is fully developed (mm HO), is the leaf area index for a given day, and is the maximum leaf area index for the plant.

When precipitation falls on any given day, the canopy storage is filled before any water is allowed to reach the ground:

and

when 2:2.1.2

and

when 2:2.1.3

where is the initial amount of free water held in the canopy on a given day (mm HO), is the final amount of free water held in the canopy on a given day (mm HO), is the amount of precipitation on a given day before canopy interception is removed (mm HO), is the amount of precipitation on a given day that reaches the soil surface (mm HO), and is the maximum amount of water that can be trapped in the canopy on a given day (mm HO).

Table 2:2-1: SWAT+ input variables used in canopy storage calculations.

Potential evapotranspiration (PET) was a concept originally introduced by Thornthwaite (1948) as part of a climate classification scheme. He defined PET is the rate at which evapotranspiration would occur from a large area uniformly covered with growing vegetation that has access to an unlimited supply of soil water and that was not exposed to advection or heat storage effects. Because the evapotranspiration rate is strongly influenced by a number of vegetative surface characteristics, Penman (1956) redefined PET as “the amount of water transpired ... by a short green crop, completely shading the ground, of uniform height and never short of water”. Penman used grass as his reference crop, but later researchers (Jensen et al., 1990) have suggested that alfalfa at a height of 30 to 50 cm may be a more appropriate choice.

Numerous methods have been developed to estimate PET. Three of these methods have been incorporated into SWAT+: the Penman-Monteith method (Monteith, 1965; Allen, 1986; Allen et al., 1989), the Priestley-Taylor method (Priestley and Taylor, 1972) and the Hargreaves method (Hargreaves et al., 1985). The model will also read in daily PET values if the user prefers to apply a different potential evapotranspiration method.

The three PET methods included in SWAT+ vary in the amount of required inputs. The Penman-Monteith method requires solar radiation, air temperature, relative humidity and wind speed. The Priestley-Taylor method requires solar radiation, air temperature and relative humidity. The Hargreaves method requires air temperature only.

Soil heat storage or release can be significant over a few hours, but is usually small from day to day because heat stored as the soil warms early in the day is lost when the soil cools late in the day or at night. Since the magnitude of daily soil heat flux over a 10- to 30-day period is small when the soil is under a crop cover, it can normally be ignored for most energy balance estimates. SWAT+ assumes the daily soil heat flux, $G$, is always equal to zero.

The aerodynamic resistance to sensible heat and vapor transfer, $r_a$, is calculated:

$r_a=\frac{ln[(z_w-d)/z_{om}]ln\lfloor(z_p-d)/z_{ov}\rfloor}{k^2u_z}$ 2:2.2.3

where $z_w$ is the height of the wind speed measurement (cm), $z_p$ is the height of the humidity (psychrometer) and temperature measurements (cm), $d$ is the zero plane displacement of the wind profile (cm), $z_{om}$ is the roughness length for momentum transfer (cm), $z_{ov}$ is the roughness length for vapor transfer (cm), $k$ is the von Kármán constant, and $u_z$ is the wind speed at height $z_w$ (m s$^{-1}$).

The von Kármán constant is considered to be a universal constant in turbulent flow. Its value has been calculated to be near 0.4 with a range of 0.36 to 0.43 (Jensen et al., 1990). A value of 0.41 is used by SWAT+ for the von Kármán constant.

Brutsaert (1975) determined that the surface roughness parameter, $z_o$, is related to the mean height ($h_c$) of the plant canopy by the relationship $h_c/z_o$ = $3e$ or 8.15 where e is the natural log base. Based on this relationship, the roughness length for momentum transfer is estimated as:

$z_{om} = h_c/8.15 =0.123*h_c$ when $h_c \le 200cm$ 2:2.2.4

$z_{om}= 0.058*(h_c)^{1.19}$ when $h_c>200cm$ 2:2.2.5

where mean height of the plant canopy ($h_c$) is reported in centimeters.

The roughness length for momentum transfer includes the effects of bluff-body forces. These forces have no impact on heat and vapor transfer, and the roughness length for vapor transfer is only a fraction of that for momentum transfer. To estimate the roughness length for vapor transfer, Stricker and Brutsaert (1978) recommended using:

$z_{ov} =0.1*z_{om}$ 2:2.2.6

The displacement height for a plant can be estimated using the following relationship (Monteith, 1981; Plate, 1971):

$d=2/3*h_c$ 2:2.2.7

The height of the wind speed measurement, $z_w$, and the height of the humidity (psychrometer) and temperature measurements, $z_p$, are always assumed to be 170 cm.

For wind speed in m s$^{-1}$, Jensen et al. (1990) provided the following relationship to calculate $K_1*0.622\lambda\rho/P$:

$K_1 *0.622*\lambda*\rho/P=1710-6.85*\overline T_{av}$ 2:2.2.19

where $\overline T_{av}$ is the mean air temperature for the day (˚C).

To calculate potential evapotranspiration, the Penman-Monteith equation must be solved for a reference crop. SWAT+ uses alfalfa at a height of 40 cm with a minimum leaf resistance of 100 s m$^{-1}$ for the reference crop. Using this canopy height, the equation for aerodynamic resistance (2:2.2.3) simplifies to:

$r_a=\frac{114.}{u_z}$ 2:2.2.20

The equation for canopy resistance requires the leaf area index. The leaf area index for the reference crop is estimated using an equation developed by Allen et al. (1989) to calculate $LAI$ as a function of canopy height. For nonclipped grass and alfalfa greater than 3 cm in height:

$LAI=1.5*ln(h_c)-1.4$ 2:2.2.21

where $LAI$ is the leaf area index and $h_c$ is the canopy height (cm). For alfalfa with a 40 cm canopy height, the leaf area index is 4.1. Using this value, the equation for canopy resistance simplifies to:

$r_c=49/(1.4-0.4*\frac{CO_2}{330})$ 2:2.2.22

The most accurate estimates of evapotranspiration with the Penman-Monteith equation are made when evapotranspiration is calculated on an hourly basis and summed to obtain the daily values. Mean daily parameter values have been shown to provide reliable estimates of daily evapotranspiration values and this is the approach used in SWAT+. However, the user should be aware that calculating evapotranspiration with the Penman-Monteith equation using mean daily values can potentially lead to significant errors. These errors result from diurnal distributions of wind speed, humidity, and net radiation that in combination create conditions which the daily averages do not replicate.

The Penman-Monteith equation combines components that account for energy needed to sustain evaporation, the strength of the mechanism required to remove the water vapor and aerodynamic and surface resistance terms. The Penman-Monteith equation is:

$\lambda E=\frac{\Delta*(H_{net}-G)+\rho_{air}*c_p*[e^o_z-e_z]/r_a}{\Delta+\gamma*(1+r_c/r_a)}$ 2:2.2.1

where $\lambda E$ is the latent heat flux density (MJ m$^{-2}$ d$^{-1}$), $E$ is the depth rate evaporation (mm d$^{-1}$), $\Delta$ is the slope of the saturation vapor pressure-temperature curve, $de/dT$ (kPa ˚C$^{-1}$), $H_{net}$ is the net radiation (MJ m$^{-2}$ d$^{-1}$), $G$ is the heat flux density to the ground (MJ m$^{-2}$ d$^{-1}$), $\rho_{air}$ is the air density (kg m$^{-3}$), $c_p$ is the specific heat at constant pressure (MJ kg$^{-1}$ ˚C$^{-1}$), is the saturation vapor pressure of air at height $z$ (kPa), $e_z$ is the water vapor pressure of air at height $z$ (kPa), $\gamma$ is the psychrometric constant (kPa ˚C$^{-1}$), $r_c$ is the plant canopy resistance (s m$^{-1}$), and $r_a$ is the diffusion resistance of the air layer (aerodynamic resistance) (s m$^{-1}$).

For well-watered plants under neutral atmospheric stability and assuming logarithmic wind profiles, the Penman-Monteith equation may be written (Jensen et al., 1990):

$\lambda E_t=\frac{\Delta*(H_{net}-G)+\gamma*K_1*(0.622*\gamma*\rho_{air}/P)*(e^o_z-e_z)/r_a}{\Delta+\gamma*(1+r_c/r_a)}$ 2:2.2.2

where $\lambda$ is the latent heat of vaporization (MJ kg$^{-1}$), $E_t$ is the maximum transpiration rate (mm d$^{-1}$), $K_1$ is a dimension coefficient needed to ensure the two terms in the numerator have the same units (for $u_z$ in m s$^{-1}$, $K_1$ = 8.64 x 104), and $P$ is the atmospheric pressure (kPa).

The calculation of net radiation, $H_{net}$, is reviewed in Chapter 1:1. The calculations for the latent heat of vaporization, $\lambda$, the slope of the saturation vapor pressure-temperature curve, $\Delta$, the psychrometric constant, $\gamma$, and the saturation and actual vapor pressure, $e^o_z$and $e_z$, are reviewed in Chapter 1:2. The remaining undefined terms are the soil heat flux, $G$, the combined term $K_1 0.622 \lambda\rho/P$, the aerodynamic resistance, $r_a$, and the canopy resistance, $r_c$.

Studies in canopy resistance have shown that the canopy resistance for a well-watered reference crop can be estimated by dividing the minimum surface resistance for a single leaf by one-half of the canopy leaf area index (Jensen et al., 1990):

$r_c=r_l/(0.5*LAI)$ 2:2.2.8

where $r_c$ is the canopy resistance (s m$^{-1}$), $r_{l}$ is the minimum effective stomatal resistance of a single leaf (s m$^{-1}$), and $LAI$ is the leaf area index of the canopy.

The distribution of stomates on a plant leaf will vary between species. Typically, stomates are distributed unequally on the top and bottom of plant leaves. Plants with stomates located on only one side are classified as hypostomatous while plants with an equal number of stomates on both sides of the leaf are termed amphistomatous. The effective leaf stomatal resistance is determined by considering the stomatal resistance of the top (adaxial) and bottom (abaxial) sides to be connected in parallel (Rosenburg et al., 1983). When there are unequal numbers of stomates on the top and bottom, the effective stomatal resistance is calculated:

$r_{l}=\frac{r_{l -ad}*r_{l-ab}}{r_{l-ab}+r_{l-ad}}$ 2:2.2.9

where $r_{l}$ is the minimum effective stomatal resistance of a single leaf (s m$^{-1}$), $r_{l-ad}$ is the minimum adaxial stomatal leaf resistance (s m$^{-1}$), $r_{l-ab}$ and is the minimum abaxial stomatal leaf resistance (s m$^{-1}$). For amphistomatous leaves, the effective stomatal resistance is:

$r_{l}=\frac{r_{l-ad}}{2}=\frac{r_{l-ab}}{2}$ 2:2.2.10

For hypostomatous leaves the effective stomatal resistance is:

$r_{l}=r_{l-ad}=r_{l-ab}$ 2:2.2.11

Leaf conductance is defined as the inverse of the leaf resistance:

$g_{l}=\frac{1}{r_{l}}$ 2:2.2.12

where $g_{l}$ is the maximum effective leaf conductance (m s$^{-1}$). When the canopy resistance is expressed as a function of leaf conductance instead of leaf resistance, equation 2:2.2.8 becomes:

$r_c=(0.5*g_l*LAI)^{-1}$ 2:2.2.13

where $r_c$ is the canopy resistance (s m$^{-1}$), $g_{l}$ is the maximum conductance of a single leaf (m s$^{-1}$), and $LAI$ is the leaf area index of the canopy.

For climate change simulations, the canopy resistance term can be modified to reflect the impact of change in CO$_2$ concentration on leaf conductance. The influence of increasing CO$_2$ concentrations on leaf conductance was reviewed by Morison (1987). Morison found that at CO$_2$ concentrations between 330 and 660 ppmv, a doubling in CO$_2$ concentration resulted in a 40% reduction in leaf conductance. Within the specified range, the reduction in conductance is linear (Morison and Gifford, 1983). Easterling et al. (1992) proposed the following modification to the leaf conductance term for simulating carbon dioxide concentration effects on evapotranspiration:

$g_{l,CO_2}=g_{l}*[1.4-0.4*(CO_2/330)]$ 2:2.2.14

where $g_{l,CO_2}$ is the leaf conductance modified to reflect CO$_2$ effects (m s$^{-1}$) and CO$_2$ is the concentration of carbon dioxide in the atmosphere (ppmv).

Incorporating this modification into equation 2:2.2.8 gives

$r_c=r_l*[(0.5*LAI)*(1.4-0.4*\frac{CO_2}{330})]^{-1}$ 2:2.2.15

SWAT+ will default the value of CO$_2$ concentration to 330 ppmv if no value is entered by the user. With this default, the term $(1.4-0.4*\frac{CO_2}{330})$reduces to 1.0 and the canopy resistance equation becomes equation 2:2.2.8.

When calculating actual evapotranspiration, the canopy resistance term is modified to reflect the impact of high vapor pressure deficit on leaf conductance (Stockle et al, 1992). For a plant species, a threshold vapor pressure deficit is defined at which the plant’s leaf conductance begins to drop in response to the vapor pressure deficit. The adjusted leaf conductance is calculated:

$g_{l}=g_{l,mx}*\lfloor1-\Delta g_{l,dcl}(vpd-vpd_{thr})\rfloor$ if $vpd>vpd_{thr}$ 2:2.2.16

$g_{l}=g_{l,mx}$ if $vpd \le vpd_{thr}$ 2:2.2.17

where $g_{l}$ is the conductance of a single leaf (m s$^{-1}$), $g_{l,mx}$ is the maximum conductance of a single leaf (m s$^{-1}$), $\Delta g_{l,dcl}$ is the rate of decline in leaf conductance per unit increase in vapor pressure deficit (m s$^{-1}$ kPa$^{-1}$), $vpd$ is the vapor pressure deficit (kPa), and $vpd_{thr}$ is the threshold vapor pressure deficit above which a plant will exhibit reduced leaf conductance (kPa). The rate of decline in leaf conductance per unit increase in vapor pressure deficit is calculated by solving equation 2:2.2.16 using measured values for stomatal conductance at two different vapor pressure deficits:

$\Delta g_{l,dcl}=\frac{(1-fr_{g,mx})}{(vpd_{fr}-vpd_{thr})}$ 2:2.2.18

where $\Delta g_{l,dcl}$ is the rate of decline in leaf conductance per unit increase in vapor pressure deficit (m s$^{-1}$ kPa$^{-1}$), $fr_{g,mx}$ is the fraction of the maximum stomatal conductance, $g_{l,mx}$, achieved at the vapor pressure deficit $vpd_{fr}$, and $vpd_{thr}$ is the threshold vapor pressure deficit above which a plant will exhibit reduced leaf conductance (kPa). The threshold vapor pressure deficit is assumed to be 1.0 kPa for all plant species.