Water that enters the soil may move along one of several different pathways. The water may be removed from the soil by plant uptake or evaporation. It can percolate past the bottom of the soil profile and ultimately become aquifer recharge. A final option is that water may move laterally in the profile and contribute to streamflow. Of these different pathways, plant uptake of water removes the majority of water that enters the soil profile.

Percolation is calculated for each soil layer in the profile. Water is allowed to percolate if the water content exceeds the field capacity water content for that layer and the layer below is not saturated. When the soil layer is frozen, no water flow out of the layer is calculated.

The volume of water available for percolation in the soil layer is calculated:

$SW_{ly,excess}=SW_{ly}-FC_{ly}$ if $SW_{ly}>FC_{ly}$ 2:3.2.1

$SW_{ly,excess}=0$ if $SW_{ly} \le FC_{ly}$ 2:3.2.2

where $SW_{ly,excess}$ is the drainable volume of water in the soil layer on a given day (mm H$_2$O), $SW_{ly}$ is the water content of the soil layer on a given day (mm H$_2$O) and $FC_{ly}$ is the water content of the soil layer at field capacity (mm H$_2$O).

The amount of water that moves from one layer to the underlying layer is calculated using storage routing methodology. The equation used to calculate the amount of water that percolates to the next layer is:

$w_{perc,ly}=SW_{ly,excess}*(1-exp[\frac{-\Delta t}{TT_{perc}}])$ 2:3.2.3

where $w_{perc,ly}$ is the amount of water percolating to the underlying soil layer on a given day (mm H$_2$O), $SW_{ly,excess}$ is the drainable volume of water in the soil layer on a given day (mm H$_2$O), $\Delta t$ is the length of the time step (hrs), and $TT_{perc}$ is the travel time for percolation (hrs). If the HRU has a seasonal high water table, percolation is not allowed when $SW_{ly+1} \le FC_{ly+1}+0.5*(SAT_{ly+1}-FC_{ly+1})$ where $SW_{ly+1}$ is the water content of the underlying soil layer (mm H$_2$O), $FC_{ly+1}$ is the water content of the underlying soil layer at field capacity (mm H$_2$O), and $SAT_{ly+1}$ is the amount of water in the underlying soil layer when completely saturated (mm H$_2$O). The water will instead stay ponded in the upper layer.

The travel time for percolation is unique for each layer. It is calculate

$TT_{perc}=\frac{SAT_{ly}-FC_{ly}}{K_{sat}}$ 2:3.2.4

where $TT_{perc}$ is the travel time for percolation (hrs), $SAT_{ly}$ is the amount of water in the soil layer when completely saturated (mm H$_2$O), $FC_{ly}$ is the water content of the soil layer at field capacity (mm H$_2$O), and $K_{sat}$ is the saturated hydraulic conductivity for the layer (mm∙h$^{-1}$).

Water that percolates out of the lowest soil layer enters the vadose zone. The vadose zone is the unsaturated zone between the bottom of the soil profile and the top of the aquifer. Movement of water through the vadose zone and into the aquifers is reviewed in Chapter 2:4.

Table 2:3-2: SWAT+ input variables used in percolation calculations.

For an HRU with a seasonal high water table, if the soil profile becomes saturated to the point that percolation for upper soil layers to lower soil layers is inhibited, water will pond in the soil profile and create a perched water table.

SWAT+ allows the user to define the depth to an impervious layer for the HRU. If the depth to the impervious layer is in the soil profile, no water is allowed to percolate out of the soil profile. If the impervious layer is defined below the soil profile, percolation out of the soil profile is adjusted from the value determined with equation 2:3.2.3 using:

2:3.4.1

where is the amount of water percolating out of the soil profile on a given day (mm HO), is the amount of water percolating out of the soil profile on a given day calculated with equation 2:3.2.3 (mm HO), and is the distance from the bottom of the soil profile to the impervious layer (m).

Water builds up in the soil profile from the bottom of the profile. After the bottom layer of the profile reaches saturation, any water exceeding the storage capacity of the bottom layer is allowed to fill the overlying layer. This continues upward until all the excess water has been distributed.

The height of the perched water table is calculated:

2:3.4.2

where is the height of the water table (mm), is the water content of the soil profile (mm HO), is the water content of the soil profile at field capacity (mm HO), is the porosity of the soil profile (mm), is the air-filled porosity expressed as a fraction, and is the depth to the impervious layer (mm).

Table 2:3-5: SWAT+ input variables used in perched water table calculations.

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

Once lateral flow is calculated, the amount of lateral flow released to the main channel is calculated:

2:3.5.10

where is the amount of lateral flow discharged to the main channel on a given day (mm HO), is the amount of lateral flow generated in the subbasin on a given day (mm HO), is the lateral flow stored or lagged from the previous day (mm HO), and is the lateral flow travel time (days).

The model will calculate lateral flow travel time or utilize a user-defined travel time. In the majority of cases, the user should allow the model to calculate the travel time. If drainage tiles are present in the HRU, lateral flow travel time is calculated:

2:3.5.11

where is the lateral flow travel time (days) and is the drain tile lag time (hrs). In HRUs without drainage tiles, lateral flow travel time is calculated:

2:3.5.12

where is the lateral flow travel time (days), is the hillslope length (m), and is the highest layer saturated hydraulic conductivity in the soil profile (mm/hr).

The expression in equation 2:3.5.10 represents the fraction of the total available water that will be allowed to enter the reach on any one day. Figure 2:3-5 plots values for this expression at different values of .

Figure 2:3-5: Influence of on fraction of lateral flow released.

The delay in release of lateral flow will smooth the streamflow hydrograph simulated in the reach.

Table 2:3-6: SWAT+ input variables used in lateral flow calculations.

Soil is comprised of three phases—solid, liquid and gas. The solid phase consists of minerals and/or organic matter that forms the matrix or skeleton of the soil. Between the solid particles, soil pores are formed that hold the liquid and gas phases. The soil solution may fill the soil pores completely (saturated) or partially (unsaturated). When the soil is unsaturated, the soil solution is found as thin films along particle surfaces, as annular wedges around contact points of particles and as isolated bodies in narrow pore passages.

The soil’s bulk density defines the relative amounts of pore space and soil matrix. Bulk density is calculated:

2:3.1.1

where is the bulk density (Mg m), is the mass of the solids (Mg), and is the total volume (m). The total volume is defined as

2:3.1.2

where is the volume of air (m), is the volume of water (m), and is the volume of solids (m). The relationship between soil porosity and soil bulk density is

2:3.1.3

where is the soil porosity expressed as a fraction of the total soil volume, is the bulk density (Mg m), and is the particle density (Mg m). The particle density, or density of the solid fraction, is a function of the mineral composition of the soil matrix. Based on research, a default value of 2.65 Mg m is used for particle density.

Storage, transport and availability of soil solution and soil air are not nearly as dependent on the total amount of porosity as they are on the arrangement of pore space. Soil pores vary in size and shape due to textural and structural arrangement. Based on the diameter of the pore at the narrowest point, the pores may be classified as macropores (narrowest diameter > 100 ), mesopores (narrowest diameter 30-100 ), and micropores (narrowest diameter < 30 ) (Koorevaar et al, 1983). Macropores conduct water only during flooding or ponding rain and drainage of water from these pores is complete soon after cessation of the water supply. Macropores control aeration and drainage processes in the soil. Mesopores conduct water even after macropores have emptied, e.g. during non-ponding rain and redistribution. Micropores retain soil solution or conduct it very slowly.

When comparing soils of different texture, clay soils contain a greater fraction of mesopores and micropores while sand soils contain mostly macropores. This is evident when the hydraulic conductivities of clay and sand soils are compared. The conductivity of a sand soil can be several orders of magnitude greater than that for a clay soil.

The water content of a soil can range from zero when the soil is oven dried to a maximum value () when the soil is saturated. For plant-soil interactions, two intermediate stages are recognized: field capacity and permanent wilting point. Field capacity is the water content found when a thoroughly wetted soil has drained for approximately two days. Permanent wilting point is the water content found when plants growing in the soil wilt and do not recover if their leaves are kept in a humid atmosphere overnight. To allow these two stages to be quantified more easily, they have been redefined in terms of tensions at which water is held by the soil. Field capacity is the amount of water held in the soil at a tension of 0.033 MPa and the permanent wilting point is the amount of water held in the soil at a tension of 1.5 MPa. The amount of water held in the soil between field capacity and permanent wilting point is considered to be the water available for plant extraction.

Table 2:3-1 lists the water content for three soils as a fraction of the total volume for different moisture conditions. Note that the total porosity, given by the water content at saturation, is lowest for the sand soil and highest for the clay soil.

The sand soil drains more quickly than the loam and clay. Only 15% of the water present in the sand soil at saturation remains at field capacity. 58% of the water present at saturation in the loam remains at field capacity while 68% of the water present at saturation in the clay soil remains at field capacity. The reduction of water loss with increase in clay content is cause by two factors. As mentioned previously, clay soils contain more mesopores and micropores than sand soils. Also, unlike sand and silt particles, clay particles possess a net negative charge. Due to the polar nature of water molecules, clay particles are able to attract and retain water molecules. The higher water retention of clay soils is also seen in the fraction of water present at permanent wilting point. In the soils listed in Table 2:3-1, the volumetric water content of the clay is 0.20 at the wilting point while the sand and loam have a volumetric water content of 0.02 and 0.05 respectively.

The plant available water, also referred to as the available water capacity, is calculated by subtracting the fraction of water present at permanent wilting point from that present at field capacity.

SWAT+ estimates the permanent wilting point volumetric water content for each soil layer as:

Water in the soil can flow under saturated or unsaturated conditions. In saturated soils, flow is driven by gravity and usually occurs in the downward direction. Unsaturated flow is caused by gradients arising due to adjacent areas of high and low water content. Unsaturated flow may occur in any direction.

SWAT+ directly simulates saturated flow only. The model records the water contents of the different soil layers but assumes that the water is uniformly distributed within a given layer. This assumption eliminates the need to model unsaturated flow in the horizontal direction. Unsaturated flow between layers is indirectly modeled with the depth distribution of plant water uptake (equation 5:2.2.1) and the depth distribution of soil water evaporation (equation 2:2.3.16).

Saturated flow occurs when the water content of a soil layer surpasses the field capacity for the layer. Water in excess of the field capacity water content is available for percolation, lateral flow or tile flow drainage unless the temperature of the soil layer is below 0°C. When the soil layer is frozen, no water movement is calculated.

Table 2:3-1: SWAT+ input variables used in percolation calculations.

Lateral flow will be significant in areas with soils having high hydraulic conductivities in surface layers and an impermeable or semipermeable layer at a shallow depth. In such a system, rainfall will percolate vertically until it encounters the impermeable layer. The water then ponds above the impermeable layer forming a saturated zone of water, i.e. a perched water table. This saturated zone is the source of water for lateral subsurface flow.

SWAT+ incorporates a kinematic storage model for subsurface flow developed by Sloan et al. (1983) and summarized by Sloan and Moore (1984). This model simulates subsurface flow in a two-dimensional cross-section along a flow path down a steep hillslope. The kinematic approximation was used in its derivation.

This model is based on the mass continuity equation, or mass water balance, with the entire hillslope segment used as the control volume. The hillslope segment has a permeable soil surface layer of depth $D_{perm}$ and length $L_{hill}$ with an impermeable soil layer or boundary below it as shown in Figure 2:3-3. The hillslope segment is oriented at an angle αhill to the horizontal.

Figure 2:3-3: Conceptual representation of the hillslope segment.

The kinematic wave approximation of saturated subsurface or lateral flow assumes that the lines of flow in the saturated zone are parallel to the impermeable boundary and the hydraulic gradient equals the slope of the bed.

Figure 2:3-4: Behavior of the water table as assumed in the kinematic storage model.

From Figure 2:3-4, the drainable volume of water stored in the saturated zone of the hillslope segment per unit area, $SW_{ly,excess}$, is

$SW_{ly,excess}=\frac{1000*H_o*\phi_d*L_{hill}}{2}$ 2:3.5.1

where $SW_{ly,excess}$ is the drainable volume of water stored in the saturated zone of the hillslope per unit area (mm H$_2$O), $H_o$ is the saturated thickness normal to the hillslope at the outlet expressed as a fraction of the total thickness (mm/mm), $\phi_d$ is the drainable porosity of the soil (mm/mm), $L_{hill}$ is the hillslope length (m), and 1000 is a factor needed to convert meters to millimeters. This equation can be rearranged to solve for $H_o$:

$H_o=\frac{2*SW_{ly,excess}}{1000*\phi_d*L_{hill}}$ 2:3.5.2

The drainable porosity of the soil layer is calculated:

$\phi_d=\phi_{soil}-\phi_{fc}$ 2:3.5.3

where $\phi_d$ is the drainable porosity of the soil layer (mm/mm), $\phi_{soil}$ is the total porosity of the soil layer (mm/mm), and $\phi_{fc}$ is the porosity of the soil layer filled with water when the layer is at field capacity water content (mm/mm).

A soil layer is considered to be saturated whenever the water content of the layer exceeds the layer’s field capacity water content. The drainable volume of water stored in the saturated layer is calculated:

$SW_{ly,excess}=SW_{ly}-FC_{ly}$ if $SW_{ly}> FC_{ly}$ 2:3.5.4

$SW_{ly,excess}=0$ if $SW_{ly} \le FC_{ly}$ 2:3.5.5

where $SW_{ly}$ is the water content of the soil layer on a given day (mm H$_2$O) and $FC_{ly}$ is the water content of the soil layer at field capacity (mm H$_2$O).

The net discharge at the hillslope outlet, $Q_{lat}$, is given by

$Q_{lat}=24*H_o*v_{lat}$ 2:3.5.6

where $Q_{lat}$ is the water discharged from the hillslope outlet (mm H$_2$O/day), $H_o$ is the saturated thickness normal to the hillslope at the outlet expressed as a fraction of the total thickness (mm/mm), $v_{lat}$ is the velocity of flow at the outlet (mm∙h$^{-1}$), and 24 is a factor to convert hours to days.

Velocity of flow at the outlet is defined as

$v_{lat}=K_{sat}*sin(\alpha_{hill})$ 2:3.5.7

where $K_{sat}$ is the saturated hydraulic conductivity (mm∙h$^{-1}$) and $\alpha_{hill}$ is the slope of the hillslope segment. The slope is input to SWAT+ as the increase in elevation per unit distance ($slp$) which is equivalent to $^{tan(\alpha_{hill})}$. Because $^{tan(\alpha_{hill}) \cong sin(\alpha_{hill})}$ , equation 2:3.5.3 is modified to use the value for the slope as input to the model:

$v_{lat}=K_{sat}*tan(\alpha_{hill})=K_{sat}*slp$ 2:3.5.8

Combining equations 2:3.5.2 and 2:3.5.8 with equation 2:3.5.6 yields the equation

$Q_{lat}=0.024*(\frac{2*SW_{ly,excess}*K_{sat}*slp}{\phi_d*L_{hill}})$ 2:3.5.9

where all terms are previously defined.

One of the most unique soil orders is the Vertisols. These soils are characterized by a propensity to shrink when dried and swell when moistened. When the soil is dry, large cracks form at the soil surface. This behavior is a result of the type of soil material present and the climate. Vertisols contain at least 30% clay with the clay fraction dominated by smectitic mineralogy and occur in areas with cyclical wet and dry periods.

Vertisols are found worldwide (Figure 2:3-1). They have a number of local names, some of which are listed in Table 2:3-3.

Figure 2:3-1: Soil associations of Vertisols (After Dudal and Eswaran, 1988)

Table 2:3-3: Alternative names for Vertisols or soils with Vertic properties (Dudal and Eswaran, 1988).

One criteria used to classify a soil as a Vertisol is the formation of shrinkage cracks in the dry season that penetrate to a depth of more than 50 cm and are at least 1 cm wide at 50 cm depth. The cracks can be considerably wider at the surface—30 cm cracks at the surface are not unusual although 6-15 cm cracks are more typical.

To accurately predict surface runoff and infiltration in areas dominated by soils that exhibit Vertic properties, the temporal change in soil volume must be quantified. Bouma and Loveday (1988) identified three soil moisture conditions for which infiltration needs to be defined (Figure 2:3-2).

Figure 2:3-2: Diagram showing the effect of wetting and drying on cracking in Vertisols (After Bouma and Loveday, 1988)

Traditional models of infiltration are applicable to soils in which cracks have been closed by swelling and the soil acts as a relatively homogenous porous medium (Condition 3 in Figure 2:3-2). Condition 1 in Figure 2:3-2 represents the driest state with cracks at maximum width, a condition present at the end of the dry season/beginning of the rainy season. Condition 2 in Figure 2:3-2 represents the crack development typical with an actively growing crop requiring multiple irrigation or rainfall events to sustain growth. Bypass flow, the vertical movement of free water along macropores through unsaturated soil horizons, will occur in conditions 1 and 2. Bypass flow (finf,2 in Figure 2:3-2) occurs when the rate of rainfall or irrigation exceeds the vertical infiltration rate into the soil peds (finf,1 in Figure 2:3-2).

When bypass flow is modeled, SWAT+ calculates the crack volume of the soil matrix for each day of simulation by layer. On days in which precipitation events occur, infiltration and surface runoff is first calculated for the soil peds (finf,1 in Figure 2:3-2) using the curve number or Green & Ampt method. If any surface runoff is generated, it is allowed to enter the cracks. A volume of water equivalent to the total crack volume for the soil profile may enter the profile as bypass flow. Surface runoff in excess of the crack volume remains overland flow.

Water that enters the cracks fills the soil layers beginning with the lowest layer of crack development. After cracks in one layer are filled, the cracks in the overlying layer are allowed to fill.

The crack volume initially estimated for a layer is calculated:

$crk_{ly,i}=crk_{max,ly}*\frac{coef_{crk}*FC_{ly}-SW_{ly}}{coef_{crk}*FC_{ly}}$ 2:3.3.1

where $crk_{ly,i}$ is the initial crack volume calculated for the soil layer on a given day expressed as a depth (mm), $crk_{max,ly}$ is the maximum crack volume possible for the soil layer (mm), $coef_{crk}$ is an adjustment coefficient for crack flow, $FC_{ly}$ is the water content of the soil layer at field capacity (mm H$_2$O), and $SW_{ly}$ is the water content of the soil layer on a given day (mm H$_2$O). The adjustment coefficient for crack flow, $coef_{crk}$, is set to 0.10.

When the moisture content of the entire profile falls below 90% of the field capacity water content for the profile during the drying stage, the crack volume for a given day is a function of the crack volume estimated with equation 2:3.3.1 and the crack volume of the layer on the previous day. When the soil is wetting and/or when the moisture content of the profile is above 90% of the field capacity water content, the crack volume for a given day is equal to the volume calculated with equation 2:3.3.1.

$crk_{ly}=\Box_{crk}*crk_{ly,d-1}+(1.0-\Box_{crk})*crk_{ly,i}$

when $SW<0.90*FC$ and $crk_{ly,i} > crk_{ly,d-1}$ 2:3.3.2

$crk_{ly}=crk_{ly,i}$ when $SW \ge 0.90*FC$ or $crk_{ly,i} \le crk_{ly,d-1}$ 2:3.3.3

where $crk_{ly}$ is the crack volume for the soil layer on a given day expressed as a depth (mm),$\Box_{crk}$ is the lag factor for crack development during drying, $crk_{ly,d-1}$ is the crack volume for the soil layer on the previous day (mm), $crk_{ly,i}$ is the initial crack volume calculated for the soil layer on a given day using equation 2:3.3.1 (mm), $SW$ is the water content of the soil profile on a given day (mm H$_2$O), and $FC$ is the water content of the soil profile at field capacity (mm H$_2$O).

As the tension at which water is held by the soil particles increases, the rate of water diffusion slows. Because the rate of water diffusion is analogous to the coefficient of consolidation in classical consolidation theory (Mitchell, 1992), the reduction in diffusion will affect crack formation. The lag factor is introduced during the drying stage to account for the change in moisture redistribution dynamics that occurs as the soil dries. The lag factor, $\Box_{crk}$, is set to a value of 0.99.

The maximum crack volume for the layer, $crk_{max,ly}$, is calculated:

$crk_{max,ly}=0.916*crk_{max}*exp\lfloor-0.0012*z_{l,ly}\rfloor*depth_{ly}$ 2:3.3.4

where $crk_{max,ly}$ is the maximum crack volume possible for the soil layer (mm), $crk_{max}$ is the potential crack volume for the soil profile expressed as a fraction of the total volume, $z_{l,ly}$ is the depth from the soil surface to the bottom of the soil layer (mm), and $depth_{ly}$ is the depth of the soil layer (mm). The potential crack volume for the soil profile, $crk_{max}$, is input by the user. Those needing information on the measurement of this parameter are referred to Bronswijk (1989; 1990).

Once the crack volume for each layer is calculated, the total crack volume for the soil profile is determined.

$crk=\sum^n_{ly=1} crk_{ly}$ 2:3.3.5

where $crk$ is the total crack volume for the soil profile on a given day (mm), $crk_{ly}$ is the crack volume for the soil layer on a given day expressed as a depth (mm), $ly$ is the layer, and $n$ is the number of layers in the soil profile.

After surface runoff is calculated for rainfall events using the curve number or Green & Ampt method, the amount of runoff is reduced by the volume of cracks present that day:

$Q_{surf}=Q_{surf,i}-crk$ if $Q_{surf,i}>crk$ 2:3.3.6

$Q_{surf}=0$ if $Q_{surf,i} \le crk$ 2:3.3.7

where $Q_{surf}$ is the accumulated runoff or rainfall excess for the day (mm H$_2$O), $Q_{surf,i}$ is the initial accumulated runoff or rainfall excess determined with the Green & Ampt or curve number method (mm H$_2$O), and $crk$ is the total crack volume for the soil profile on a given day (mm). The total amount of water entering the soil is then calculated:

$w_{inf}=R_{day}-Q_{surf}$ 2:3.3.8

where $w_{inf}$ is the amount of water entering the soil profile on a given day (mm H$_2$O), $R_{day}$ is the rainfall depth for the day adjusted for canopy interception (mm H$_2$O), and $Q_{surf}$ is the accumulated runoff or rainfall excess for the day (mm H$_2$O).

Bypass flow past the bottom of the profile is calculated:

$w_{crk,btm}=0.5*crk*(\frac{crk_{ly=nn}}{depth_{ly=nn}})$ 2:3.3.9

where $w_{crk,btm}$ is the amount of water flow past the lower boundary of the soil profile due to bypass flow (mm H$_2$O), $crk$ is the total crack volume for the soil profile on a given day (mm), $crk_{ly=nn}$ is the crack volume for the deepest soil layer ($ly=nn$) on a given day expressed as a depth (mm), and $depth_{ly=nn}$ is the depth of the deepest soil layer ($ly=nn$) (mm).

After $w_{crk,btm}$ is calculated, each soil layer is filled to field capacity water content beginning with the lowest layer and moving upward until the total amount of water entering the soil, $w_{inf}$, has been accounted for.

Table 2:3-4: SWAT+ input variables used in bypass flow calculations.