Soil water processes

Soil water balances

The soil water processes in SPHY are modeled for three soil layers (Figure 2), being (i) the first soil layer (root zone), (ii) second soil layer (subzone), and (iii) third soil layer (groundwater layer). The water balance of the first soil layer is:

Equation 28

$$SW_{1,t}=SW_{1,t-1}+Pe_t-ET_{a,t}-RO_t-LF_{1,t}-Perc_{1,t}+Cap_t$$

with $SW_{l,t}$ and $SW_{l,t-1}$ the water content in the first soil layer on days $t$ and $t-1$, respectively, $Pe_t$ the effective precipitation on day $t$, $ET_{a,t}$ the actual evapotranspiration on day $t$, $RO_t$ the surface runoff on day $t$, $LF_{1,t}$ the lateral flow from the first soil layer on day $t$, $Perc_{1,t}$ the percolation from the first to the second soil layer on day $t$, and $Cap_t$ the capillary rise from the second to the first soil layer on day $t$. The second soil layer water balance is:

Equation 29

$$SW_{2,t}=SW_{2,t-1}+Perc_{1,t}-Perc_{2,t}-Cap_t$$

with $SW_{2,t}$ and $SW_{2,t-1}$ the water content in the second soil layer on day $t$ and $t-1$, respectively, and $Perc_{2,t}$ percolation from the second to the third soil layer on day $t$. The third soil layer water balance is given as:

Equation 30

$$SW_{3,t}=SW_{3,t-1}+Gchrg_{t}-BF_t$$

with $SW_{3,t}$ and $SW_{3,t-1}$ the water content in the third soil layer on day $t$ and $t-1$, respectively, $Gchrg_t$ groundwater recharge from the second to the third soil layer on day $t$, and $BF_t$ baseflow on day $t$. If the glacier module is used, then groundwater recharge consists of percolation from the second soil layer and percolated glacier melt; otherwise, only percolation from the second soil layer is taken into account.

The user can opt to run SPHY without the third soil layer (groundwater). This may be desirable if the user for example is mainly interested in simulating soil moisture conditions in the root zone, instead of evaluating for instance the contribution of baseflow to the total routed river flow. In that case, only the two upper soil layers are used where the bottom boundary of soil layer two is controlled by a seepage flux (positive outward), and instead of baseflow from the third soil layer, water leaves the second soil layer through lateral flow. With the groundwater module turned off, the water balance for the second soil layer is:

Equation 31

$$SW_{2,t}=SW_{2,t-1}+Perc_{1,t}-LF_{2,t}-Cap_{t}-Seep$$

with $LF_{2,t}$ lateral flow from the second soil layer, and Seep seepage in or out of the second soil layer (positive is outgoing). The units for all water balance terms are in mm.

Actual evapotranspiration

Evapotranspiration refers to both the transpiration from vegetation and the evaporation from soil or open water. As was mentioned in Section 2.3, the Kc accounts for both the crop transpiration and soil evaporation. The additional use of the dynamic vegetation module accounts for a time-variable vegetation cover, meaning that the role of evaporation becomes more dominant as soon as vegetation cover decreases.

Many limiting factors (e.g., salinity stress, water shortage, water excess, diseases) can cause a reduction in potential evapotranspiration$(ET_p)$, resulting in the actual evapotranspiration rate $(ET_a)$. Since SPHY is a water-balance model, SPHY only accounts for stresses related to water shortage or water excess. If there is too much water in the soil profile, then the plant is unable to extract water because of oxygen stress (Bartholomeus et al. 2008). The calculation of evapotranspiration reduction due to water excess (oxygen stress) is quite complex and requires a substantial number of plant and soil properties (e.g., soil temperature, root dry weight, plant respiration, and minimum gas filled soil porosity; (Bartholomeus et al. 2008)) that are generally not available for the spatial scale that SPHY is applied on. Therefore, SPHY uses an evapotranspiration reduction parameter $(ETred_{wet})$ that has a value of 0 if the soil is saturated, and otherwise it will have a value of 1. This parameter is used in the following equation to calculate the actual evapotranspiration:

Equation 32

$$ET_{a,t}=ET_{p,t} \cdot ETred_{wet} \cdot ETred_{dry}$$

with $ET_{a,t} (mm)$ the actual evapotranspiration on day $t$, $ET_{p,t} (mm)$ the potential evapotranspiration on day $t$, and $ETred_{wet}$ and $ETred_{dry}$ the reduction parameters for water excess and water shortage conditions, respectively. $ETred_{dry}$ is either calculated using the Feddes equation (Feddes et al., 1978) or the plant water stress method (Allen et al., 1998). The Feddes equation assumes a linear decline in rootwater uptake if the water pressure head drops below a critical value. This critical value can be determined using the soil water retention curve (pF curve), which relates the moisture content of the soil to its binding capacity. This relation is unique for each soil type. The binding capacity is a suction force (H) and is therefore often expressed in cm negative water column. The pF value is simply a conversion of the suction force (H), and is calculated as:

Equation 33

$$pF=log_{10}(-H)$$

Soils that are at field capacity generally have a pF of 2, meaning –100cm of water column, and soils that are at permanent wilting point have a pF of 4.2, or –16000cm of water column. The permanent wilting point is often referred to as the point where the crop dies. In SPHY it is assumed that the linear decline in rootwater uptake starts at a pF of 3 (–1000cm water column). Therefore, $ETred_{dry}$ (–) is calculated as:

Equation 34

$$ETred_{dry,t}= \frac{SW_{1,t}-SW_{1,pF4.2}}{SW_{1,pF3}-SW_{1,pF4.2}}$$

with $ETred_{dry,t} (-)$ the reduction in rootwater uptake due to water shortage on day $t$, $SW_{1,t} (mm)$ the actual soil water content in the first soil layer on day $t$, and $SW_{1,pF3} (mm)$ and $SW_{1,pF4.2} (mm)$ the soil water content in the first soil layer at pF3 and pF4.2, respectively. $ETred_{dry}$ can therefore have values ranging between 0 and 1, where a value of 1 represents optimal plant growing conditions, and 0 means no rootwater uptake at all.

The plant water stress method (Allen et al., 1998) accounts for plant-specific characteristics, in addition to soil hydraulic properties. With this method $ETred_{dry}$ (–) is determined with the following equation:

Equation 35

$$ETred_{dry,t}=\frac{TAW-D_r}{(1-p)TAW}$$

Where $ETred_{dry,t}$ is the reduction parameter for water shortage (-), TAW is the total available water in the rootzone (mm), the root zone depletion (mm) and p the depletion fraction (-). The total available TAW water is defined as:

Equation 36

$$TAW=SW_{l,fc}-SW_{1,pF4.2}$$

Where $SW_{l,fc}$ is the soil water content at field capacity (mm) and $SW_{1,pF4.2}$ the soil water content at wilting point (mm). The root zone depletion $D_r$ is defined as:

Equation 37

$$D_r=SW_{l,fc}-SW_{1,t}$$

Where $SW_{1,t}$ is the current soil water content of the first soil layer (mm). The depletion fraction p is defined as the fraction of TAW that a crop can extract from the root zone without suffering water stress, which is determined by the following equation:

Equation 38

$$p = p_{tab}+0.04(5-ET_{p,t})$$

Where $p_{tab}$ is a landuse-specific tabular value of the depletion fraction (-) and $ET_{p,t}$ is the potential evapotranspiration (mm). Values for the landuse-specific tabular value of the depletion fraction can be obtained from Allen et al., 1998 (Table 22).

$ETred_{dry}$ is eventually used in Equation 32 to calculate the $ET_a$.

Open-water evaporation

Open-water evaporation is determined in the open-water cells. In these cells all soil hydraulic processes are turned off and runoff equals precipitation minus open-water evaporation. Reservoir cells cannot dry up, i.e. we assume that there is always water present in the reservoir cells. Open-water evaporation is determined as follows:

Equation 39

$$ET_{open-water}=Kc_{open-water}ET_r$$

Where $Kc_{open-water}$ is the crop coefficient value for open-water evaporation (-) and $ET_r$ is the reference evapotranspiration (mm). We suggest to set $Kc_{open-water}$ to a value of 1.2, after Allen et al. (1998). In each time step the open-water evaporation is subtracted from the reservoir and lake storage.

Surface runoff

SPHY accounts for both infiltration excess (Hortonian runoff; Horton, 1933) and saturation excess surface runoff (Hewlettian runoff; Hewlett, 1961). Infiltration excess surface runoff occurs when the precipitation intensity exceeds the infiltration capacity of the first soil layer; hence, it is a sub-daily process. Since SPHY runs with a daily time step, we have developed a new infiltration excess surface runoff equation, which is inspired by the Green-Ampt formula (Heber Green and Ampt, 1911). We assumed a constant infiltration rate $f (mm \space hr^{-1})$, which is determined for each cell and each day by:

Equation 40

$$f=\frac{K_{eff}}{24} \Big[{1+\frac{SW_{1,sat}-SW_1}{SW_{1,sat}}}\Big]$$

where $K_{eff}$ is the effective hydraulic conductivity $(mm \space day^{-1})$, $SW_{1}(mm)$ the water content in the first soil layer, $SW_{1,sat}(mm)$ the saturated water content of the first soil layer and $\lambda$ is a calibration parameter (-). Bouwer (1969) suggested an approximation of $K_{eff} \approx 0.5 K_{sat}$.

Infiltration excess surface runoff occurs when the precipitation intensity exceeds the infiltration rate f (Beven, 2012). We assume that highest precipitation intensity is recorded in the first hour of the rain storm and decreases linearly until the end of the storm. Furthermore, we assume a triangular-shaped precipitation intensity $p(t) (mm \space hr^{-1})$ according to:

Equation 40

$$p(t)=\frac{1}{2} \alpha^{2}Pt+\alpha P$$

where $\alpha$ is the fraction of daily rainfall that occurs in the hour with the highest intensity (-), P is the daily rainfall (mm), and t is an hourly time step (hr). Daily infiltration excess surface runoff $Q_{surf}$ is determined as follows:

Equation 40

$$Q_{surf}=\begin{cases} \frac{(\alpha P-f)^{2}}{\alpha^2P} &\text{if } \alpha P>f \\ 0 &\text{if } \alpha P \le f \end{cases}$$

When the hourly precipitation intensity $\alpha P$ is higher than the infiltration rate f, surface runoff equals the triangular shaped area of the precipitation above the infiltration rate. The amount of precipitation below the infiltration rate will infiltrate into the rootzone.

Saturation excess surface runoff occurs when the first soil layer gets saturated and is calculated as:

Equation 40

$$RO = \begin{Bmatrix} SW_1 -SW_{1,sat} & \text{if }& SW_1>SW_{1,sat} \\ 0 & \text{if } & SW_1 \le SW_{1,sat} \end{Bmatrix}$$

with RO (mm) surface runoff, $SW_1 (mm)$ the water content in the first soil layer, and $SW_{1,sat}(mm)$ the saturated water content of the first soil layer.

Lateral flow

Lateral flow is substantial in catchments with steep gradients and soils with high hydraulic conductivities (Beven 1981; Beven and Germann 1982; Sloan and Moore 1984). In SPHY, it is assumed that only the amount of water exceeding field capacity can be used for lateral flow. Therefore, the drainable volume of water (excess water) needs to be calculated first:

Equation 41

$$W_{l,exc} = \begin{Bmatrix} SW_l -SW_{l,fc} & \text{if }& SW_l>SW_{l,fc} \\ 0 & \text{if } & SW_l \le SW_{l,fc} \end{Bmatrix}$$

with $W_{l,exc} (mm)$ the drainable volume of water from soil layer $l, SW_l (mm)$ the water content in soil layer $l$, and $SW_{l,fc} (mm)$ the field capacity of soil layer $l$. According to Sloan and Moore (1984), the lateral flow at the hillslope outlet can be calculated as:

Equation 42

$$LF^{*}_{l}=W_{l,excfrac} \cdot \upsilon_{lat,l}$$

with $LF_l^* (mm)$ lateral flow from soil layer layer $l$, $W_{l,excfrac}$(–) the drainable volume of water as a fraction of the saturated volume, and $\upsilon_{lat,l} (mm \cdot d^{-1})$ the flow velocity at the outlet. In SPHY, the drainable volume as a fraction of the saturated volume is calculated as:

Equation 43

$$W_{l,excfrac}=\frac{W_{l,exc}}{SW_{l,sat}-SW_{l,fc}}$$

The velocity of flow at the outlet, $\upsilon_{lat,l} (mm \cdot d^{-1})$, depends on both the saturated hydraulic conductivity $K_{sat,l} (mm \cdot d^{-1})$ and the slope of the hill slp (–), and is defined as:

Equation 44

$$\upsilon_{lat,l}=K_{sat,l} \cdot slp$$

The slope (slp) in SPHY is calculated for each grid cell as the increase in elevation per unit distance.

According to Neitsch et al. (2009), only a fraction of lateral flow will reach the main channel on the day it is generated if the catchment of interest has a time of concentration greater than 1 day. This concept is also implemented in the SPHY model, and uses a lateral flow travel time $TT_{lag,l} (d)$ to lag a portion of lateral flow release to the channel:

Equation 45

$$LF_l=(LF_l^*+LF^*_{l,t-1})\cdot \Big(1-exp \Big[\frac{-1}{TT_{lag,l}} \Big]\Big)$$

with $LF_{l}(mm)$the amount of lateral flow entering the channel on a given day,$LF_{l}^*(mm)$ the lateral flow (Equation 42) generated within the cell on a given day, and $LF_{l,t-1}^*(mm)$ the lateral flow lagged from the previous day. SPHY assumes the lateral flow travel time to be dependent on the field capacity $SW_{l,fc}(mm)$, saturated content $SW_{l,sat}(mm)$, and the saturated conductivity $K_{sat,l}(mm \cdot d^{-1})$, according to:

Equation 46

$$TT_{lag,l}=\frac{SW_{l,sat}-SW_{l,fc}}{K_{sat,l}}$$

A longer lateral flow travel time will result in a smoother streamflow hydrograph.

Percolation

If the groundwater module is used, then water can percolate from the first to the second soil layer, and from the second to the third soil layer. If the user decides to run SPHY without the groundwater module, percolation only occurs from the first to the second soil layer. In SPHY, water can only percolate if the water content exceeds the field capacity of that layer, and the water content of the underlying layer is not saturated. A similar approach has been used in the SWAT model (Neitsch et al. 2009). The water volume available for percolation to the underlying layer is calculated as:

Equation 47

$$W_{l,exc} = \begin{Bmatrix} 0, & \text{if }& SW_l \le SW_{l,fc} \text{or} SW_{l+1} \ge SW_{l+1,sat}\\ SW_{l+1,sat}-SW_{l+1},&\text{if} &SW_l-SW_{l,fc}>SW_{l+1,sat}-SW_{l+1} \\ SW_l-SW_{l,fc}, & else \end{Bmatrix}$$

with $W_{l,exc}(mm)$ the drainable volume of water from layer $l$, $SW_l (mm)$the water content in layer $l$, $SW_{l,fc} (mm)$ the field capacity of layer $l$, $SW_{l+1} (mm)$ the water content in layer $l+1$, and $SW_{l+1,sat} (mm)$ the saturated water content of layer $l+1$. Only a certain amount of $W_{l,exc}$ will percolate to the underlying soil layer, depending on the percolation travel time $TT_{perc,l}(d)$. This approach follows the storage routing methodology, which is also implemented in the SWAT model (Neitsch et al. 2009):

Equation 48

$$w_{l,perc}=W_{l,exc} \cdot \Big(1-exp \Big [ \frac{-1}{TT_{perc,l}} \Big] \Big)$$

with $w_{l,perc}(mm)$ the amount of water percolating to the underlying soil layer. Since the speed at which water can move through the soil is mainly dependent on the saturated hydraulic conductivity$(K_{sat})$, the travel time for percolation is calculated the same way as the travel time for lateral flow (Equation 46).

Groundwater recharge

Water that percolates from the second to the third soil layer will eventually reach the shallow aquifer. This process is referred to as groundwater recharge hereafter. If the glacier module is used as well, then glacier melt that percolates also contributes to the groundwater recharge. Groundwater recharge often does not occur instantaneously, but with a time lag that depends on the depth of the groundwater table and soil characteristics. SPHY uses the same exponential decay weighting function as proposed by Venetis (1969) and used by Sangrey, Harrop-Williams, and Klaiber (1984) in a precipitation groundwater response model. This approach has also been adopted in the SWAT model (Neitsch et al. 2009), using:

Equation 49

$$Gchrg_t= \Big(1-exp^{\frac{-1}{\delta_{gw}}}\Big) \cdot w_{2,perc}+exp^{\frac{-1}{\delta_{gw}}} \cdot Gchrg_{t-1}$$

with $Gchrg_t (mm)$ and $Gchrg_{t-1}(mm)$ the groundwater recharge on days $t$ and $t-1$, respectively. $\delta_{gw}(d)$ is the delay time and $w_{2,perc}(mm)$ is the amount of water that percolates from the second to the third layer on day $t$.

Baseflow

After groundwater recharge has been calculated, SPHY calculates baseflow, which is defined as the flow going from the shallow aquifer to the main channel. Baseflow only occurs when the amount of water stored in the third soil layer exceeds a certain threshold $(BF_{thresh})$ that can be specified by the user. Baseflow calculation in SPHY is based on the steady-state response of groundwater flow to recharge (Hooghoudt 1940) and the water table fluctuations that are a result of the non-steady response of groundwater flow to periodic groundwater recharge (Smedema and Rycroft 1983). The SWAT model (Neitsch et al. 2009) assumes a linear relation between the variation in groundwater flow (baseflow) and the rate of change in water table height, according to:

Equation 50

$$\frac{dBF}{dt}=10 \cdot \frac{K_{sat}}{\mu L^2_{gw}}\cdot(Gchrg-BF)=\alpha_{gw} \cdot (Gchrg-BF)$$

with BF (mm) the groundwater flow (baseflow) into the main channel on day $t$, $K_{sat} (mm \cdot d^{-1})$ the hydraulic conductivity of the shallow aquifer, $\mu(-)$ the specific yield of the shallow aquifer, $L_{gw}(m)$ the distance from the subbasin divide for the groundwater system to the main channel, $Gchrg(mm)$ the amount of groundwater (Equation 49) recharge entering the shallow aquifer on day $t$, and $\alpha_{gw} (-)$ the baseflow recession coefficient. Equation 50 can be integrated and rearranged to calculate baseflow, according to:

Equation 51

$$BF_t= \begin{Bmatrix} 0, & \text{if }& SW_3 \le BF_{thresh} \\ BF_{t-1}\cdot exp^{-\alpha_{gw}}+Gchrg_t \cdot (1-exp^{-\alpha_{gw}}), & \text{if } & SW_3>BF_{thresh} \end{Bmatrix}$$

with$BF_t (mm)$ the baseflow into the channel on day $t$, and $BF_{t-1} (mm)$ the baseflow into the channel on day $t-1$. Since this equation has proven its success in the SWAT model (Neitsch et al. 2009) throughout many applications worldwide, this equation has been adopted in the SPHY model as well.

The baseflow recession coefficient $(\alpha_{gw})$ is an index that relates the baseflow response to changes in groundwater recharge. Lower values for $\alpha_{gw}$ therefore correspond to areas that respond slowly to groundwater recharge, whereas higher values indicate areas that have a rapid response to groundwater recharge. The baseflow recession coefficient is generally used as a calibration parameter in the SPHY model, but a good first approximation of this coefficient can be calculated using the number of baseflow days (Neitsch et al. 2009):

Equation 52

$$\alpha_{gw}=\frac{2.3}{BFD}$$

with BFD (d) the number of baseflow days, which is defined as the number of days required for baseflow recession to decline.