The classification of a stream as ephemeral, intermittent or perennial is a function of the amount of groundwater contribution received by the stream. Ephemeral streams contain water during and immediately after a storm event and are dry the rest of the year. Intermittent streams are dry part of the year, but contain flow when the groundwater is high enough as well as during and after a storm event. Perennial streams receive continuous groundwater contributions and flow throughout the year.

During periods when a stream receives no groundwater contributions, it is possible for water to be lost from the channel via transmission through the side and bottom of the channel. Transmission losses are estimated with the equation

$tloss=K_{ch}*TT*P_{ch}*L_{ch}$ 7:1.5.1

where $tloss$ are the channel transmission losses (m$^3$ H$_2$O), $K_{ch}$ is the effective hydraulic conductivity of the channel alluvium (mm/hr), $TT$ is the flow travel time (hr), $P_{ch}$ is the wetted perimeter (m), and $L_{ch}$ is the channel length (km). Transmission losses from the main channel are assumed to enter bank storage or the deep aquifer.

Typical values for $K_{ch}$ for various alluvium materials are given in Table 7:1-4. For perennial streams with continuous groundwater contribution, the effective conductivity will be zero.

Table 7:1-5: SWAT+ input variables that pertain to transmission losses.

Open channel flow is defined as channel flow with a free surface, such as flow in a river or partially full pipe. SWAT+ uses Manning’s equation to define the rate and velocity of flow. Water is routed through the channel network using the variable storage routing method or the Muskingum river routing method. Both the variable storage and Muskingum routing methods are variations of the kinematic wave model. A detailed discussion of the kinematic wave flood routing model can be found in Chow et al. (1988).

The variable storage routing method was developed by Williams (1969) and used in the HYMO (Williams and Hann, 1973) and ROTO (Arnold et al., 1995) models.

For a given reach segment, storage routing is based on the continuity equation:

$V_{in}-V_{out}=\Delta V_{stored}$ 7:1.3.1

where $V_{in}$ is the volume of inflow during the time step (m$^3$ H$_2$O), $V_{out}$ is the volume of outflow during the time step (m$^3$ H$_2$O), and $\Delta V_{stored}$ is the change in volume of storage during the time step (m$^3$ H$_2$O). This equation can be written as

$\Delta t*(\frac{q_{in,1}+q_{in,2}}{2})-\Delta t*(\frac{q_{out,1}+q_{out,2}}{2})=V_{stored,2}-V_{stored,1}$ 7:1.3.2

where $\Delta t$ is the length of the time step (s), $q_{in,1}$ is the inflow rate at the beginning of the time step (m$^3$/s), $q_{in,2}$ is the inflow rate at the end of the time step (m$^3$/s), $q_{out,1}$ is the outflow rate at the beginning of the time step (m$^3$/s), $q_{out,2}$ is the outflow rate at the end of the time step (m$^3$/s), $V_{stored,1}$ is the storage volume at the beginning of the time step (m$^3$ H$_2$O), and $V_{stored,2}$ is the storage volume at the end of the time step (m$^3$ H$_2$O). Rearranging equation 7:1.3.2 so that all known variables are on the left side of the equation,

$q_{in,ave}+\frac{V_{stored,1}}{\Delta t}-\frac{q_{out,1}}{2}=\frac{V_{stored,2}}{\Delta t}+\frac{q_{out,2}}{2}$ 7:1.3.3

where $q_{in,ave}$ is the average inflow rate during the time step: $q_{in,ave}=\frac{q_{in,1}+q_{in,2}}{2}$.

Travel time is computed by dividing the volume of water in the channel by the flow rate.

$TT=\frac{V_{stored}}{q_{out}}=\frac{V_{stored,1}}{q_{out,1}}=\frac{V_{stored,2}}{q_{out,2}}$ 7:1.3.4

where $TT$ is the travel time (s), $V_{stored}$ is the storage volume (m$^3$ H$_2$O), and $q_{out}$ is the discharge rate (m$^3$/s).

To obtain a relationship between travel time and the storage coefficient, equation 7:1.3.4 is substituted into equation 7:1.3.3:

$q_{in,ave}+\frac{V_{stored,1}}{(\frac{\Delta t}{TT})*(\frac{V_{stored,1}}{q_{out,1}})}-\frac{q_{out,1}}{2}=\frac{V_{stored,2}}{(\frac{\Delta t}{TT})*(\frac{V_{stored,2}}{q_{out,2}})}+\frac{q_{out,2}}{2}$ 7:1.3.5

which simplifies to

$q_{out,2}=(\frac{2*\Delta t}{2*TT+\Delta t})*q_{in,ave}+(1-\frac{2*\Delta t}{2*TT+ \Delta t})*q_{out,1}$ 7:1.3.6

This equation is similar to the coefficient method equation

$q_{out,2}=SC*q_{in,ave}+(1-SC)*q_{out,1}$ 7:1.3.7

where $SC$ is the storage coefficient. Equation 7:1.3.7 is the basis for the SCS convex routing method (SCS, 1964) and the Muskingum method (Brakensiek, 1967; Overton, 1966). From equation 7:1.3.6, the storage coefficient in equation 7:1.3.7 is defined as

$SC=\frac{2*\Delta t}{2*TT+ \Delta t}$ 7:1.3.8

It can be shown that

$(1-SC)*q_{out}=SC*\frac{V_{stored}}{\Delta t}$ 7:1.3.9

Substituting this into equation 7:1.3.7 gives

$q_{out,2}=SC*(q_{in,ave}+\frac{V_{stored,1}}{\Delta t })$ 7:1.3.10

To express all values in units of volume, both sides of the equation are multiplied by the time step

$V_{out,2}=SC*(V_{in}+V_{stored,1})$ 7:1.3.11

Manning’s equation for uniform flow in a channel is used to calculate the rate and velocity of flow in a reach segment for a given time step:

$q_{ch}=\frac{A_{ch}*R_{ch}^{2/3}*slp_{ch}^{1/2}}{n}$ 7:1.2.1

$v_c=\frac{R_{ch}^{2/3}*slp_{ch}^{1/2}}{n}$ 7:1.2.2

where $q_{ch}$ is the rate of flow in the channel (m$^3$/s), $A_{ch}$ is the cross-sectional area of flow in the channel (m$^2$), $R_{ch}$ is the hydraulic radius for a given depth of flow (m), $slp_{ch}$ is the slope along the channel length (m/m), $n$ is Manning’s “n” coefficient for the channel, and $v_c$ is the flow velocity (m/s).

SWAT+ routes water as a volume. The daily value for cross-sectional area of flow, $A_{ch}$, is calculated by rearranging equation 7:1.1.7 to solve for the area:

$A_{ch}=\frac{V_{ch}}{1000*L_{ch}}$ 7:1.2.3

where $A_{ch}$ is the cross-sectional area of flow in the channel for a given depth of water (m$^2$), $V_{ch}$ is the volume of water stored in the channel (m$^3$), and $L_{ch}$ is the channel length (km). Equation 7:1.1.4 is rearranged to calculate the depth of flow for a given time step:

$depth=\sqrt{\frac{A_{ch}}{z_{ch}}+(\frac{W_{btm}}{2*z_{ch}})^2}-\frac{W_{btm}}{2*z_{ch}}$ 7:1.2.4

where $depth$ is the depth of flow (m), $A_{ch}$ is the cross-sectional area of flow in the channel for a given depth of water (m$^2$), $W_{btm}$ is the bottom width of the channel (m), and $z_{ch}$ is the inverse of the channel side slope. Equation 7:1.2.4 is valid only when all water is contained in the channel. If the volume of water in the reach segment has filled the channel and is in the flood plain, the depth is calculated:

$depth=depth_{bnkfull}+\sqrt{\frac{(A_{ch}-A_{ch,bnkfull})}{z_{fld}}+(\frac{W_{btm,fld}}{2*z_{fld}})^2}-\frac{W_{btm,fld}}{2*z_{fld}}$ 7:1.2.5

where $depth$ is the depth of flow (m), $depth_{bnkfull}$ is the depth of water in the channel when filled to the top of the bank (m), $A_{ch}$ is the cross-sectional area of flow in the channel for a given depth of water (m$^2$), $A_{ch,bnkfull}$ is the cross-sectional area of flow in the channel when filled to the top of the bank (m$^2$), $W_{btm,fld}$ is the bottom width of the flood plain (m), and $z_{fld}$ is the inverse of the flood plain side slope.

Once the depth is known, the wetting perimeter and hydraulic radius are calculated using equations 7:1.1.5 (or 7:1.1.10) and 7:1.1.6. At this point, all values required to calculate the flow rate and velocity are known and equations 7:1.2.1 and 7:1.2.2 can be solved.

Table 7:1-2: SWAT+ input variables that pertain to channel flow calculations.

SWAT+ assumes the main channels, or reaches, have a trapezoidal shape (Figure 7:1-1).

Users are required to define the width and depth of the channel when filled to the top of the bank as well as the channel length, slope along the channel length and Manning’s “n” value. SWAT+ assumes the channel sides have a 2:1 run to rise ratio ($z_{ch}$ = 2). The slope of the channel sides is then ½ or 0.5. The bottom width is calculated from the $bankfull$ width and depth with the equation:

$W_{btm}=W_{bnkfull}-2*z_{ch}*depth_{bnkfull}$ 7:1.1.1

where $W_{btm}$ is the bottom width of the channel (m), $W_{bnkfull}$ is the top width of the channel when filled with water (m), $z_{ch}$ is the inverse of the channel side slope, and $depth_{bnkfull}$ is the depth of water in the channel when filled to the top of the bank (m). Because of the assumption that $z_{ch}=2$, it is possible for the bottom width calculated with equation 7:1.1.1 to be less than or equal to zero. If this occurs, the model sets $W_{btm}=0.5*W_{bnkfull}$ and calculates a new value for the channel side slope run by solving equation 7:1.1.1 for $z_{ch}$:

$z_{ch}=\frac{(W_{bnkfull}-W_{btm})}{2*depth_{bnkfull}}$ 7:1.1.2

For a given depth of water in the channel, the width of the channel at water level is:

$W=W_{btm}+2*z_{ch}*depth$ 7:1.1.3

where $W$ is the width of the channel at water level (m), $W_{btm}$ is the bottom width of the channel (m), $z_{ch}$ is the inverse of the channel slope, and $depth$ is the depth of water in the channel (m). The cross-sectional area of flow is calculated:

$A_{ch}=(W_{btm}+z_{ch}*depth)*depth$ 7:1.1.4

where $A_{ch}$ is the cross-sectional area of flow in the channel (m$^2$), $W_{btm}$ is the bottom width of the channel (m), $z_{ch}$ is the inverse of the channel slope, and $depth$ is the depth of water in the channel (m). The wetted perimeter of the channel is defined as

$P_{ch}=W_{btm}+2*depth*\sqrt{1+z_{ch}^2}$ 7:1.1.5

where $P_{ch}$ is the wetted perimeter for a given depth of flow (m). The hydraulic radius of the channel is calculated

$R_{ch}=\frac{A_{ch}}{P_{ch}}$ 7:1.1.6

where $R_{ch}$ is the hydraulic radius for a given depth of flow (m), $A_{ch}$ is the cross-sectional area of flow in the channel (m$^2$), and $P_{ch}$ is the wetted perimeter for a given depth of flow (m). The volume of water held in the channel is

$V_{ch}=1000*L_{ch}*A_{ch}$ 7.1.1.7

where $V_{ch}$ is the volume of water stored in the channel (m$^3$), $L_{ch}$ is the channel length (km), and $A_{ch}$ is the cross-sectional area of flow in the channel for a given depth of water (m$^2$).

When the volume of water in the reach exceeds the maximum amount that can be held by the channel, the excess water spreads across the flood plain. The flood plain dimensions used by SWAT+ are shown in Figure 7:1-2.

The bottom width of the floodplain, $W_{btm,fld}$, is $W_{btm,fld}=5*W_{bnkfull}$. SWAT+ assumes the flood plain side slopes have a 4:1 run to rise ratio ($z_{fld}$ = 4). The slope of the flood plain sides is then ¼ or 0.25.

When flow is present in the flood plain, the calculation of the flow depth, cross-sectional flow area and wetting perimeter is a sum of the channel and floodplain components:

$depth=depth_{bnkfull}+depth_{fld}$ 7:1.1.8

$A_{ch}=(W_{btm}+z_{ch}*depth_{bnkfull})*depth_{bnkfull}+(W_{btm,fld}+z_{fld}+depth_{fld})*depth_{fld}$

7:1.1.9

$P_{ch}=W_{btm}+2*depth_{bnkfull}*\sqrt{1+z_{ch}^2}+4*W_{bnkfull}+2*depth_{fld}*\sqrt{1+z_{fld}^2}$

7:1.1.10

where $depth$ is the total depth of water (m), $depth_{bnkfull}$ is the depth of water in the channel when filled to the top of the bank (m), $depth_{fld}$ is the depth of water in the flood plain (m), $A_{ch}$ is the cross-sectional area of flow for a given depth of water (m$^2$), $W_{btm}$ is the bottom width of the channel (m), $z_{ch}$ is the inverse of the channel side slope, $W_{btm,fld}$ is the bottom width of the flood plain (m), $z_{fld}$ is the inverse of the flood plain side slope, $P_{ch}$ is the wetted perimeter for a given depth of flow (m), and $W_{bnkfull}$ is the top width of the channel when filled with water (m).

Table 7:1-1: SWAT+ input variables that pertain to channel dimension calculations.

The Muskingum routing method models the storage volume in a channel length as a combination of wedge and prism storages (Figure 7:1-3).

When a flood wave advances into a reach segment, inflow exceeds outflow and a wedge of storage is produced. As the flood wave recedes, outflow exceeds inflow in the reach segment and a negative wedge is produced. In addition to the wedge storage, the reach segment contains a prism of storage formed by a volume of constant cross-section along the reach length.

As defined by Manning’s equation (equation 7:1.2.1), the cross-sectional area of flow is assumed to be directly proportional to the discharge for a given reach segment. Using this assumption, the volume of prism storage can be expressed as a function of the discharge, $K*q_{out}$, where $K$ is the ratio of storage to discharge and has the dimension of time. In a similar manner, the volume of wedge storage can be expressed as $K*X*(q_{in}-q_{out})$, where $X$ is a weighting factor that controls the relative importance of inflow and outflow in determining the storage in a reach. Summing these terms gives a value for total storage

$V_{stored}=K*q_{out}+K*X*(q_{in}-q_{out})$ 7:1.4.1

where $V_{stored}$ is the storage volume (m$^3$ H$_2$O), $q_{in}$ is the inflow rate (m$^3$/s), $q_{out}$ is the discharge rate (m$^3$/s), $K$ is the storage time constant for the reach (s), and $X$ is the weighting factor. This equation can be rearranged to the form

$V_{stored}=K*(X*q_{in}+(1-X)*q_{out})$ 7:1.4.2

This format is similar to equation 7:1.3.7.

The weighting factor, $X$, has a lower limit of 0.0 and an upper limit of 0.5. This factor is a function of the wedge storage. For reservoir-type storage, there is no wedge and $X=0.0$. For a full-wedge, $X=0.5$. For rivers, $X$ will fall between 0.0 and 0.3 with a mean value near 0.2.

The definition for storage volume in equation 7:1.4.2 can be incorporated into the continuity equation (equation 7:1.3.2) and simplified to

$q_{out,2}=C_1*q_{in,2}+C_2*q_{in,1}+C_3*q_{out,1}$ 7:1.4.3

where $q_{in,1}$ is the inflow rate at the beginning of the time step (m$^3$/s), $q_{in,2}$ is the inflow rate at the end of the time step (m$^3$/s), $q_{out,1}$ is the outflow rate at the beginning of the time step (m$^3$/s), $q_{out,2}$ is the outflow rate at the end of the time step (m$^3$/s), and

$C_1=\frac{\Delta t-2*K*X}{2*K*(1-X)+\Delta t}$ 7:1.4.4

$C_2=\frac{\Delta t+2*K*X}{2*K*(1-X)+ \Delta t}$ 7:1.4.5

$C_3=\frac{2*K*(1-X)- \Delta t}{2*K*(1-X)+\Delta t}$ 7:1.4.6

where $C_1+C_2+C_3=1$. To express all values in units of volume, both sides of equation 7:1.4.3 are multiplied by the time step

$V_{out,2}=C_1*V_{in,2}+C_2*V_{in,1}+C_3*V_{out,1}$ 7:1.4.7

To maintain numerical stability and avoid the computation of negative outflows, the following condition must be met:

$2*K*X<\Delta t<2*K*(1-X)$ 7:1.4.8

The value for the weighting factor, $X$, is input by the user. The value for the storage time constant is estimated as:

$K=coef_1*K_{bnkfull}+coef_2*K_{0.1bnkfull}$ 7:1.4.9

where $K$ is the storage time constant for the reach segment (s), $coef_1$ and $coef_2$ are weighting coefficients input by the user, $K_{bnkfull}$ is the storage time constant calculated for the reach segment with bankfull flows (s), and $K_{0.1bnkfull}$ is the storage time constant calculated for the reach segment with one-tenth of the bankfull flows (s). To calculate $K_{bnkfull}$ and $K_{0.1bnkfull}$, an equation developed by Cunge (1969) is used:

$K=\frac{1000*L_{ch}}{c_k}$ 7:1.4.10

where $K$ is the storage time constant (s), $L_{ch}$ is the channel length (km), and $c_k$ is the celerity corresponding to the flow for a specified depth (m/s). Celerity is the velocity with which a variation in flow rate travels along the channel. It is defined as

$c_k=\frac{d}{dA_{ch}}(q_{ch})$ 7:1.4.11

where the flow rate, $q_{ch}$, is defined by Manning’s equation. Differentiating equation 7:1.2.1 with respect to the cross-sectional area gives

$c_k=\frac{5}{3}*(\frac{R_{ch}^{2/3}*slp_{ch}^{1/2}}{n})=\frac{5}{3}*v_c$ 7:1.4.12

where $c_k$ is the celerity (m/s), $R_{ch}$ is the hydraulic radius for a given depth of flow (m), $slp_{ch}$ is the slope along the channel length (m/m), n is Manning’s “$n$” coefficient for the channel, and $v_c$ is the flow velocity (m/s).

Table 7:1-3: SWAT+ input variables that pertain to Muskingum routing.

Evaporation losses from the reach are calculated:

7:1.6.1

where is the evaporation from the reach for the day (m HO), is an evaporation coefficient, is potential evaporation (mm HO), is the channel length (km), is the channel width at water level (m), and is the fraction of the time step in which water is flowing in the channel.

The evaporation coefficient is a calibration parameter for the user and is allowed to vary between 0.0 and 1.0.

The fraction of the time step in which water is flowing in the channel is calculated by dividing the travel time by the length of the time step.

Table 7:1-6: SWAT+ input variables that pertain to evaporation losses.

Water storage in the reach at the end of the time step is calculated:

7:1.8.1

where is the volume of water in the reach at the end of the time step (m HO), is the volume of water in the reach at the beginning of the time step (m HO), is the volume of water flowing into the reach during the time step (m HO), is the volume of water flowing out of the reach during the time step (m HO), is the volume of water lost from the reach via transmission through the bed (m HO), is the evaporation from the reach for the day (m HO), is the volume of water added or removed from the reach for the day through diversions (m HO), and is the volume of water added to the reach via return flow from bank storage (m HO).

SWAT+ treats the volume of outflow calculated with equation 7:1.3.11 or 7:1.4.7 as the net amount of water removed from the reach. As transmission losses, evaporation and other water losses for the reach segment are calculated, the amount of outflow to the next reach segment is reduced by the amount of the loss. When outflow and all losses are summed, the total amount will equal the value obtained from 7:1.3.11 or 7:1.4.7.

The amount of water entering bank storage on a given day is calculated:

7:1.7.1

where is the amount of water entering bank storage (m HO), are the channel transmission losses (m HO), and is the fraction of transmission losses partitioned to the deep aquifer.

Bank storage contributes flow to the main channel or reach within the subbasin. Bank flow is simulated with a recession curve similar to that used for groundwater. The volume of water entering the reach from bank storage is calculated:

7:1.7.2

where is the volume of water added to the reach via return flow from bank storage(m HO), is the total amount of water in bank storage (m HO), and is the bank flow recession constant or constant of proportionality.

Water may move from bank storage into an adjacent unsaturated zone. SWAT+ models the movement of water into adjacent unsaturated areas as a function of water demand for evapotranspiration. To avoid confusion with soil evaporation and transpiration, this process has been termed ‘revap’. This process is significant in watersheds where the saturated zone is not very far below the surface or where deep-rooted plants are growing. ‘Revap’ from bank storage is governed by the groundwater revap coefficient defined for the last HRU in the subbasin.

The maximum amount of water than will be removed from bank storage via ‘revap’ on a given day is:

7:1.7.3

where is the maximum amount of water moving into the unsaturated zone in response to water deficiencies (m HO), is the revap coefficient, is the potential evapotranspiration for the day (mm HO), is the channel length (km), and is the width of the channel at water level (m). The actual amount of revap that will occur on a given day is calculated:

if 7:1.7.4

if 7:1.7.5

where is the actual amount of water moving into the unsaturated zone in response to water deficiencies (m HO), is the maximum amount of water moving into the unsaturated zone in response to water deficiencies (m HO), and is the amount of water in bank storage at the beginning of day (m HO).

Table 7:1-7: SWAT+ input variables that pertain to bank storage.