Chapter 5 Forest hydrology

Forest hydrological processes are key for soil and plant water balances. This chapter details the design and implementation of most of the vertical hydrological processes included in package medfate. These processes determine the value of the water flows mentioned in eq. (3.1). The only process that is not described here is woody plant transpiration (\(Tr_w\)), as it will be covered in detail in chapter 6.

Although each process can be simulated in specific sub-model functions, function hydrology_verticalInputs() allows simulating soil water input processes altogether, including canopy water interception, snow accumulation/melt, soil infiltration and percolation.

5.1 Snow and rainfall

Precipitation (\(P\)) is considered be snow precipitation (\(Ps\)) when \(T_{mean}<0\), and is considered rainfall (\(Pr\)) otherwise. Thus, we have: \[\begin{equation} P = Pr + Ps \end{equation}\]

Interception of snow by the canopy is neglected, and all snow is assumed to accumulate in a single storage compartment \(S_{snow}\) over the soil (i.e. canopy snow storage capacity is neglected). Rainfall interception is described in section 5.3.

5.2 Snow pack dynamics

A very simple snow submodel is used for snow pack dynamics (accumulation and melt), taken from Kergoat (1998). When mean air temperature is above 0 Celsius (\(T_{mean}>0\)), a simple energy budget relates snow melt, \(Sm\) (mm), to air temperature and soil-level radiation (see function hydrology_snowMelt()): \[\begin{equation} Sm = \frac{Rad\cdot L^{SWR}_{ground}\cdot (1-\alpha_{ice}) + \tau_{day} \cdot T_{mean} \cdot \rho_{air} \cdot C_p/r_{s}}{\lambda_{ice}} \end{equation}\] where \(Rad\) is solar radiation (\(MJ \cdot m^{-2}\)), \(L^{SWR}_{ground}\) is the fraction of (short-wave) radiation reaching the ground, \(\alpha_{ice} = 0.9\) is the albedo of snow, \(\tau_{day} = 86400\) is the day duration in seconds, \(\rho_{air}\) is the air density (\(kg \cdot m^{-3}\)), depending on temperature and elevation (see utility functions of the meteoland reference manual), \(C_{p} = 1013.86 \cdot 10^{-6}\, MJ \cdot kg^{-1} \cdot C^{-1}\) is the specific heat capacity of the air and \(r_{s} = 100\,s \cdot m^{-1}\) is the snow aerodynamic resistance and \(\lambda_{ice} = 0.33355\, MJ \cdot kg^{-1}\) is the (latent) heat of fusion of snow.

5.3 Rainfall interception loss

As mentioned above, interception loss is only modelled for liquid precipitation (i.e. snow interception is not modelled). Rainfall interception loss, \(In\), is estimated using either the Gash et al. (1995) or Liu (2001) models.

5.3.1 Gash (1995) model

In the Gash et al. (1995) analytical interception model for sparse canopies, rain is assumed to fall in a single event during the day. First, the amount of rainfall needed to saturate the canopy is calculated: \[\begin{equation} P_G = - \frac{S_{canopy}/C_{canopy}}{ER_{ratio}} \cdot \ln(1-ER_{ratio}) \end{equation}\] where \(S_{canopy}\) is the canopy water storage capacity (in mm) – i.e. the minimum amount of water needed to saturate the canopy –, \(C_{canopy}\) is the canopy cover and \(ER_{ratio}\) is the ratio of evaporation rate to rainfall intensity during the rainfall event (\(R_{int}\); in \(mm \cdot h^{-1}\)). Although interception models are normally applied to single-canopy stands, we apply the sparse Gash model to the whole stand (including shrubs and the herbaceous layer). Moreover, in our implementation stem interception is lumped with crown interception, so that \(S_{canopy}\) represents both.

Following Watanabe & Mizutani (1996) we estimate \(S_{canopy}\), the canopy water storage capacity, from adjusted LAI values: \[\begin{equation} S_{canopy}=LAI_{herb}+\sum_{i}{s_{water,i}\cdot LAI_{i}^{\phi}} \end{equation}\] where \(s_{water,i}\) is the depth of water that can be retained by leaves, branches and stems of cohort \(i\) per unit of leaf area index (\(mm \cdot LAI^{-1}\)). To estimate the stand cover, \(C_{canopy}\), we use the complement of the percentage of PAR that reaches the ground, i.e. \(C_{canopy} = 1 - L^{PAR}_{ground}\) (Deguchi et al. 2006).

The amount of water evaporated from interception (\(In\), in mm), is calculated as: \[\begin{eqnarray} In = C_{canopy}\cdot P_G+C_{canopy}\cdot ER_{ratio}\cdot(Pr-P_G) \: {if}\: Pr > P_G \\ In = C_{canopy}\cdot Pr\: {if}\: Pr \leq P_G \end{eqnarray}\] where \(Pr\) is the daily gross rainfall (in mm). Net rainfall (\(Pr_{net}\), also in mm) includes throughfall and stemflow, and is calculated as the difference between gross rainfall and interception loss. Fig. 5.1 below shows examples of relative throughfall (including stemflow), calculated according to the Gash et al (1995) interception model, under different situations (see function hydrology_rainInterception()).

Examples of canopy interception with different \(S_{canopy}\) (canopy water storage capacity), \(ER_{ratio}\) (ratio between evaporation and rainfall rates) and \(p\) (throughfall coefficient; \(p = 1 - C_{canopy}\)), using the Gash (1995) model.

Figure 5.1: Examples of canopy interception with different \(S_{canopy}\) (canopy water storage capacity), \(ER_{ratio}\) (ratio between evaporation and rainfall rates) and \(p\) (throughfall coefficient; \(p = 1 - C_{canopy}\)), using the Gash (1995) model.

5.3.2 Liu (2001) model

In the single-storm form of the Liu (2001) model is similar to the Gash et al. (1995) model in terms of its inputs. Rainfall interception during a single storm can be predicted by gross rainfall amount \(P_G\), the ratio of evaporation rate to rainfall rate during the rainfall event, \(ER_{ratio}\), canopy gap fraction (\(C_{canopy}\)) and storage capacity (\(S_{canopy}\)):

\[\begin{equation} In = S_{canopy} \cdot \left[ 1.0 - \exp\left(-P_G \cdot \left( \frac{C_{canopy}} {S_{canopy}} \right) \right) \right] \cdot (1.0 - (ER/C_{canopy})) + (ER\cdot P_G) \end{equation}\]

Fig. 5.2 below shows examples of relative throughfall (including stemflow), calculated according to the Liu et al (2001) interception model, under different situations (see function hydrology_rainInterception()).

Examples of canopy interception with different \(S_{canopy}\) (canopy water storage capacity), \(ER_{ratio}\) (ratio between evaporation and rainfall rates) and \(p\) (throughfall coefficient; \(p = 1 - C_{canopy}\)), using the Liu (2001) model.

Figure 5.2: Examples of canopy interception with different \(S_{canopy}\) (canopy water storage capacity), \(ER_{ratio}\) (ratio between evaporation and rainfall rates) and \(p\) (throughfall coefficient; \(p = 1 - C_{canopy}\)), using the Liu (2001) model.

5.4 Runon, runoff, infiltration and percolation

Water from rainfall events is divided between runoff and infiltration, and infiltrated water is distributed among soil layers assuming a circulation through large macropores. In contrast, snow melt (\(Sm\)) does not participate from infiltration/runoff partitioning and percolation, because melting rates are assumed to be slow enough to infiltrate completely in the topmost layer.

The amount of rainfall event water that reaches the soil is the sum of net rainfall (\(Pr_{net}\)) and runon (\(Ro\), in mm), if any. The amount of rainfall event water infiltrating into the soil is \(If = Pr_{net} + Ro - Ru\), where \(Ru\) (in mm) is the rainfall water lost by surface runoff. Infiltration can be estimated using either Boughton (1989) or the Green & Ampt (1911) model.

5.4.1 Boughton (1989) model

Surface runoff \(Ru\) can be estimated using the USDA SCS curve number method, following Boughton (1989): \[\begin{equation} Ru=\frac{(Pr_{net} + Ro - 0.2 \cdot V_{fc, soil})^2}{(Pr_{net} + Ro - 0.8 \cdot V_{fc, soil})} \end{equation}\] where \(V_{fc, soil}\) (in mm) is the overall soil water retention capacity. The Boughton (1989) model is available in function hydrology_infiltrationBoughton() and examples of its behavior are given in the figure below.

Examples of infiltration/runoff calculation using the Boughton (1989) model for different values of net rainfall and overall retention capacity, \(V_{soil}\), calculated from different soil depths (topsoil+subsoil) and assuming that soil texture is 15% clay and 25% sand. Rock fragment content was 25% and 40% for the topsoil and subsoil, respectively.

Figure 5.3: Examples of infiltration/runoff calculation using the Boughton (1989) model for different values of net rainfall and overall retention capacity, \(V_{soil}\), calculated from different soil depths (topsoil+subsoil) and assuming that soil texture is 15% clay and 25% sand. Rock fragment content was 25% and 40% for the topsoil and subsoil, respectively.

5.4.2 Green-Ampt (1911) model

The description of the Green & Ampt (1911) model has been adapted from chapter 9 of Bonan (2019). The Green-Ampt equation uses the concept of a horizontal wetting front at some distance \(L_w\) into the soil. The soil above \(L_w\) is saturated (i.e. \(\theta_{sat}\)), whreas soil below the wetting front has a uniform content \(\theta_{dry}\) corresponding to conditions prior to the rainfall event. Darcy’s law provides the flux of water \(i_c\) infiltrating into the soil from the surface to the wetting front at depth \(z = -L_w\):

\[\begin{equation} i_c = K_{sat} \left[\frac{|\Psi_w| }{L_w} + 1 \right] \end{equation}\] where \(\Psi_w\) is the matric potential at the wetting front. The total cumulative infiltrated water \(I_c\) at time \(t\) is the time integral of \(i_c\), and is also equal to the \(L_w\) multiplied by the change in water, i.e. \(\Delta \theta = \theta_{sat} - \theta_{dry}\): \[\begin{equation} I_c = \Delta \theta L_w \end{equation}\] The infiltration rate can also be expressed as the time derivative of \(I_c\), i.e.: \[\begin{equation} i_c = \frac{dI_c}{dt} = \Delta \theta \frac{dL_w}{dt} \end{equation}\] Combining the two equations for \(i_c\) we obtain: \[\begin{equation} \Delta \theta \frac{dL_w}{dt} = K_{sat} \left[\frac{|\Psi_w| }{L_w} + 1 \right] \end{equation}\] and substituting \(I_c/\Delta \theta\) for \(L_w\) gives: \[\begin{equation} \frac{dI_c}{dt} = K_{sat} \left[\frac{|\Psi_w| \Delta \theta}{I_c} + 1 \right] \end{equation}\] By rearranging terms and integrating we obtain an expression for \(I_c(t)\): \[\begin{equation} I_c(t) = K_{sat}t + |\Psi_w| \Delta \theta \ln\left[1 + \frac{I_c(t)}{|\Psi_w| \Delta \theta} \right] \end{equation}\] This equation defines the cumulative infiltration to time \(t\) as a function of \(K_{sat}\), \(\Psi_{w}\), \(\theta_{sat}\) and \(\theta_{dry}\). Because \(I_c\) appears on both sides of the equation, it must be found using numerical root-finding (the Newton-Raphson method in our case). To estimate \(\Psi_{w}\) the following equation is used: \[\begin{equation} \Psi_{w} = \left(\frac{2b + 3}{2b + 6} \right)\Psi_{sat} \end{equation}\]

Parameter \(\theta_{dry}\) is the volumetric content of the topmost soil layer (antecendent soil conditions), and parameters \(\theta_{sat}\), \(\Psi_{sat}\), \(K_{sat}\) and \(b\) are estimated using Campbell’s water retention curves from the texture of the topmost layer. The infiltration time \(t\) (in \(h\)) is estimated from the sum of net rainfall (\(Pr_{net}\)) and runon (\(Ro\)) (both in \(mm\)) using rainfall intensity of the rainfall event (\(R_{int}\); in \(mm/h\)). The Green & Ampt (1911) model is available in function hydrology_infiltrationGreenAmpt() and examples of its behavior are given in the figure below.

Examples of infiltration using the Green-Ampt (1911) model for different values of net rainfall. The left panel shows the infiltration for different soil textural types, assuming a \(\theta_{dry} = 0.3\). The right panel shows the infiltration for different values of \(\theta_{dry}\), assuming a silt loam soil. In both cases rainfall intensity is assumed \(R_{int} = 6\,mm/h\).

Figure 5.4: Examples of infiltration using the Green-Ampt (1911) model for different values of net rainfall. The left panel shows the infiltration for different soil textural types, assuming a \(\theta_{dry} = 0.3\). The right panel shows the infiltration for different values of \(\theta_{dry}\), assuming a silt loam soil. In both cases rainfall intensity is assumed \(R_{int} = 6\,mm/h\).

5.4.3 Percolation

Rainfall infiltrating into the soil is divided among soil layers. Following Granier et al. (1999), part of the water reaching one soil layer percolates quickly through the macropores whereas the other fills micropores (and may contribute to the vertical, unsaturated, flow). The amount of water reaching the micropores of each soil layer \(s\), i.e. \(If_s\), is modelled using an extinction function that depends on macroporosity (see function hydrology_infiltrationRepartition).

5.5 Bare soil evaporation

Evaporation from the soil surface is the last component of the soil water balance to be calculated before calculating plant transpiration. Bare soil evaporation cannot happen if there is snow over the soil surface (i.e., if \(S_{snow}>0\)).

Potential evaporation from the soil (\(PE_{soil}\); in \(mm \cdot day^{-1}\)) is defined as the product between \(PET\) and \(L^{SWR}_{ground}\), the proportion of SWR absorbed by the ground: \[\begin{equation} PE_{soil} = PET \cdot L^{SWR}_{ground} \end{equation}\]

Actual evaporation from the soil surface is modeled as in Mouillot et al. (2001), who in turn followed Ritchie (1972). First, the model determines \(t_{soil}\), the time needed to evaporate the current water deficit (difference between field capacity and current moisture) in the surface soil layer: \[\begin{equation} t_{soil} = \left \{ \frac{V_{fc,1}\cdot(1- W_1)}{\gamma_{soil}} \right \} \end{equation}\] where \(V_{fc,1}\) is the water retention capacity of layer 1, \(W_1\) is the proportion of moisture in relation to field capacity of layer 1 and \(\gamma_{soil}\) is the maximum daily evaporation (\(mm \cdot day^{-1}\)). The calculated time is used to determine the ‘supplied’ evaporation, \(SE_{soil}\): \[\begin{equation} SE_{soil} = \gamma_{soil} \cdot (\sqrt{t_{soil}+1}-\sqrt{1}) \end{equation}\] The amount of water actually evaporated from the soil, \(Es\), is then calculated as the minimum between supply, \(SE_{soil}\), and demand (Federer 1982), i.e. \(PE_{soil}\) (see function hydrology_soilEvaporationAmount):

\[\begin{equation} Es = \min(PE_{soil}, SE_{soil}) \end{equation}\]

Figure 5.5 shows the cumulative evaporation from soils for different values of maximum evaporation rate.

Cumulative bare soil evaporation for different values of maximum evaporation rate (\(\gamma_{soil}\)). Topsoil layer (0 – 30 cm) is initialized at field capacity (\(V_1 = 50 mm\)). \(PE_{soil}\) was assumed not to be limiting.

Figure 5.5: Cumulative bare soil evaporation for different values of maximum evaporation rate (\(\gamma_{soil}\)). Topsoil layer (0 – 30 cm) is initialized at field capacity (\(V_1 = 50 mm\)). \(PE_{soil}\) was assumed not to be limiting.

5.6 Transpiration of the herbaceous layer

Transpiration of herbaceous layer is modelled analogously to the transpiration of woody plants. First, potential evapotranspiration from the herbaceous layer (\(PET_{herb}\); in \(mm \cdot day^{-1}\)) is defined as the product between \(PET\) and \(L^{SWR}_{herb}\), the proportion of SWR reaching the herbaceous layer: \[\begin{equation} PET_{herb} = PET \cdot L^{SWR}_{herb} \end{equation}\] Thus, the more dense the woody plant canopy, the lower herbaceous transpiration will be. Maximum herbaceous transpiration \(Tr_{herb, \max}\) depends on both \(PET_{herb}\) and the amount of transpiring surface, i.e. \(LAI_{herb}\). To estimate \(Tr_{\max}\) the model uses the empirical equation of Granier et al. (1999), as in section 6.1.1: \[\begin{equation} \frac{Tr_{herb,\max}}{PET}= -0.006\cdot (LAI_{herb})^2+0.134\cdot LAI_{herb} \end{equation}\]

Finally, actual herbaceous transpiration for a given soil layer \(s\) (\(Tr_{herb,s}\)) is reduced according to the soil water potential (\(\Psi_s\)), assuming that transpiration reduction follows Weibull function with fixed coefficients:

\[\begin{equation} Tr_{herb, s}= FRP_{herb,s} \cdot Tr_{herb,\max} \cdot \exp \left \{\ln{(0.5)}\cdot \left[ \frac{\Psi_i}{-1.5} \right] ^2 \right \} \tag{5.1} \end{equation}\] where \(FRP_{herb,s}\) is the proportion of fine roots in layer \(s\), defined using \(Z_{50} = 50 mm\) and \(Z_{95} = 500 mm\) (see 2.4.7).

5.7 Vertical water movement within the soil

In this chapter we describe water movement within the soil, whose design was taken from chapter 8 of Bonan (2019), and much of the present text has been summarized and adapted from this source.

5.7.1 Darcian flow and Richards equation

Water flows in the soil from high to low potential as described by Darcy’s law. The total potential is the sum of the matric potential (\(\Psi_{s}\), here expressed in \(m\), see 2.3.3) and vertical depth \(z\) (also in \(m\)). The vertical depth \(z\) is taken as positive in the upward direction, so that \(z = 0\) is the ground surface. In the vertical dimension the flux \(Q\), i.e. rate of water movement (in \(m s^{-1}\)), according to Darcy’s law is: \[\begin{equation} Q = -K(\theta)\cdot \frac{\delta (\Psi + z)}{\delta z} = -K(\theta)\cdot \frac{\delta \Psi}{\delta z} - K(\theta) \end{equation}\] where \(K(\theta)\) is the hydraulic conductivity (in \(m s^{-1}\)) that describes the ease with which water can move through the porous media space, and depends on soil moisture (\(\theta\)).

Consider a volume of soil with horizontal area \(\Delta x \Delta y\), thickness \(\Delta z\) and with non-rock proportion \(\lambda = 1 - P_{rocks}\). The mass flux entering the parcel from above is \(\rho_{w} Q_{in} \Delta x \Delta y\) and the flux out of it is \(\rho_{w} Q_{out} \Delta x \Delta y\). Conservation requires that the difference between fluxes and source/sinks equals the rate of change in water storage. The mass of water in the soil volume is \(\rho_{w} \theta \Delta x \Delta y \Delta z \lambda\) so the change in water over a time interval \(\Delta t\) is:

\[\begin{equation} \rho_{w} \frac{\Delta \theta}{\Delta t} \Delta x \Delta y \Delta z \lambda = - \rho_{w} (Q_{in} - Q_{out}) \Delta x \Delta y \end{equation}\]

and gives the continuity equation:

\[\begin{equation} \frac{\Delta \theta}{\Delta t} \lambda = - \frac{(Q_{in} - Q_{out})}{\Delta z} \end{equation}\]

which expressed in the notation of calculus is:

\[\begin{equation} \frac{\delta \theta}{\delta t} \lambda = - \frac{\delta Q}{\delta z} \end{equation}\]

Substituting Darcy’s law for \(Q\) gives the Richards equation:

\[\begin{equation} \frac{\delta \theta}{\delta t} \lambda = \frac{\delta}{\delta z} \left[ K(\theta) \frac{\delta \Psi}{\delta z} + K(\theta) \right] = K(\theta) \frac{\delta^2 \Psi}{\delta^2 z} + \frac{\delta K}{\delta z} \frac{\delta \Psi}{\delta z} + \frac{\delta K}{\delta z} \end{equation}\]

The former is known as the mixed-form Richards equation, because it includes \(\theta\) on the left hand and the vertical gradient in matric potential \(\Psi\) on the right hand. The head-based or \(\Psi\)-based form expresses the left hand term using \(\Psi\) as dependent variable:

\[\begin{equation} C(\Psi) \frac{\delta \Psi}{\delta t} \lambda = K(\Psi) \frac{\delta^2 \Psi}{\delta^2 z} + \frac{\delta K}{\delta z} \frac{\delta \Psi}{\delta z} + \frac{\delta K}{\delta z} \end{equation}\]

where \(C(\Psi) = \delta \theta / \delta \Psi\) is known as the specific moisture capacity, in \(m^{-1}\).

5.7.2 Finite difference approximation

The finite difference approximation for the Richards equation represents the soil as a network of \(N\) discrete nodal points that vary in space and time. The network is presented in Fig. 5.6.

Schematic representation of network of soil nodes

Figure 5.6: Schematic representation of network of soil nodes

For reasons of numerical stability, Richards equation is solved using an implicit time discretization win which the spatial derivatives \(\delta K / \delta z\), \(\delta \Psi / \delta z\) and \(\delta^2 \Psi / \delta ^2 z\) are written numerically using a central difference approximation at time \(n+1\), and the time derivative uses a backward difference approximation at \(n+1\).

Soil layer \(i\) has a thickness \(\Delta z_i\). Water content \(\theta_i\), matric potential \(\Psi_i\) and hydraulic conductivity \(K_i\) are defined at the center of the layer at depth \(z_i\) and are uniform over the layer. Here \(K_i\) values have been corrected by \(\lambda_i\), to account for the fact that the amount of stones in the soil reduce the conductivity, which is expressed per soil area unit (so, in fact, all terms in the mass balance equations are multiplied by \(\lambda_i\), except the source/sinks). Water inputs/outputs from other processes (i.e. infiltration, snow melt, bare soil evaporation, herbaceous transpiration and woody plant uptake/redistribution) are represented as a source/sink term \(S_i\) for each node \(i\).

For a given intermediate layer \(1 < i < N\), the outward flux \(Q_i\) between layers \(i\) and \(i+1\) is defined using:

\[\begin{equation} Q_i = -\frac{K_{i+1/2}}{\Delta z_{i+1/2}} (\Psi_i - \Psi_{i+1}) - K_{i+1/2} \end{equation}\]

where \(z_{i+1/2} = z_i - z_{i+1}\) and \(K_{i+1/2} = 0.5(K_i + K_{i+1})\) replaces the vertically varying hydraulic conductivity within the two layers with an effective conductivity for an equivalent homogeneous medium. The inward flux \(Q_{i-1}\) between layers \(i-1\) and \(i\) is defined similarly:

\[\begin{equation} Q_{i-1} = -\frac{K_{i-1/2}}{\Delta z_{i-1/2}} (\Psi_{i-1} - \Psi_{i}) - K_{i-1/2} \end{equation}\]

where \(z_{i-1/2} = z_{i-1} - z_{i}\) and \(K_{i-1/2} = 0.5(K_{i-1} + K_{i})\). If we apply the fluxes to the water balance equation, with fluxes and capacitance expressed for time \(n+1\), and considering the \(S_i\) source/sink term (in \(m s^{-1}\)), we obtain the following \(\Psi\)-based form:

\[\begin{equation} \frac{\Delta z_i}{\Delta t} C_i^{n+1} (\Psi_i^{n+1} - \Psi_i^n) \lambda_i = \frac{K_{i-1/2}^{n+1}}{\Delta z_{i-1/2}} (\Psi_{i-1}^{n+1} - \Psi_i^{n+1}) - \frac{K_{i+1/2}^{n+1}}{\Delta z_{i+1/2}} (\Psi_i^{n+1} - \Psi_{i+1}^{n+1}) \\ + K_{i-1/2}^{n+1} - K_{i+1/2}^{n+1} + S^{n+1}_i \end{equation}\]

For the first layer (\(i=1\)) the boundary condition is defined so that there is no inward flux (but a source/sink term is added), which yields the following \(\Psi\)-based form:

\[\begin{equation} \frac{\Delta z_1}{\Delta t} C_1^{n+1} (\Psi_1^{n+1} - \Psi_1^n) \lambda_1 = - \frac{K_{1+1/2}^{n+1}}{\Delta z_{1+1/2}} (\Psi_1^{n+1} - \Psi_{2}^{n+1}) - K_{1+1/2}^{n+1} + S^{n+1}_1 \end{equation}\]

Finally, for the last layer (\(i=N\)) a free drainage condition is specified as boundary condition, where \(Q_N = K_N\), which yields the following \(\Psi\)-based form:

\[\begin{equation} \frac{\Delta z_N}{\Delta t} C_N^{n+1} (\Psi_N^{n+1} - \Psi_N^n) \lambda_N = \frac{K_{N-1/2}^{n+1}}{\Delta z_{N-1/2}} (\Psi_{N-1}^{n+1} - \Psi_{N}^{n+1}) + K_{N-1/2}^{n+1} - K_{N}^{n+1} + S^{n+1}_N \end{equation}\]

The numerical form of the \(\Psi\)-based Richards equation for a soil with \(N\) layers is a tridiagonal system of \(N\) equations with \(N\) unknown values \(\Psi\) at time \(n+1\). An additional complexity is that \(\Psi\) at \(n+1\) also appears through \(K\) and \(C\). The solution requires linearizing with respect to \(\Psi\) through various methods. Here the predictor-corrector method is employed, a two-step solution that solves the equation twice (Haverkamp et al. 1977). The predictor step uses explicit linearization to solve for \(\Psi\) over one-half of a full time step (\(\Delta t/2\)) at time \(n +1/2\) with \(K\) and \(C\) from time \(n\). The resulting values for \(\Psi\) are used to evaluate \(K\) and \(C\) at \(n+1/2\). The corrector step uses these to obtain \(\Psi\) over the full time step (\(\Delta t\)) at \(n+1\). An implicit solution is used for the predictor step and the Crank-Nicolson method is used for the corrector step.

Source/sink terms are defined as follows for the first layer (\(i=1\)): \[\begin{equation} S_1 = If_1 + Sm - Es - Tr_{herb,1} - Ex_{woody,1} \end{equation}\] and for \(1 < i \leq N\): \[\begin{equation} S_i = If_i - Tr_{herb,i} - Ex_{woody,i} \end{equation}\] where \(Sm\) is snow melt, \(Es\) is bare soil evaporation, \(If_i\) is infiltration amount into layer \(i\), \(Tr_{herb, i}\) is herbaceous uptake from layer \(i\) and \(Ex_{woody}\) is woody plant uptake from layer \(i\). All right-hand terms are converted from \(mm/day\) to instantaneous rates (\(m/s\)) for proper units in \(S\).

Function hydrology_soilFlows() implements the vertical water movement flows, given a soil and source/sink vector for each layer. 24 time steps of 1 hour are used for vertical water simulation of each day.

Bibliography

Bonan, G. (2019). Climate change and terrestrial ecosystem modeling. Cambridge University Press, Cambridge, UK.
Boughton, W. (1989). A review of the USDA SCS curve number method. Australian Journal of Soil Research, 27, 511–523.
Deguchi, A., Hattori, S. & Park, H.-T. (2006). The influence of seasonal changes in canopy structure on interception loss: Application of the revised Gash model. Journal of Hydrology, 318, 80–102.
Federer, C. (1982). Transpirational supply and demand: Plant, soil, and atmospheric effects evaluated by simulation. Water Resources Research, 18, 355–362.
Gash, J., Lloyd, C. & Lachaud, G. (1995). Estimating sparse forest rainfall interception with an analytical model. Journal of Hydrology, 170.
Granier, A., Bréda, N., Biron, P. & Villette, S. (1999). A lumped water balance model to evaluate duration and intensity of drought constraints in forest stands. Ecological Modelling, 116, 269–283.
Green, H.W. & Ampt, G.A. (1911). Studies on Soil Physics. The Journal of Agricultural Science, 4, 1–24.
Haverkamp, R., Vauclin, M., Touma, J., Wierenga, P.J. & Vachaud, G. (1977). A Comparison of Numerical Simulation Models For OneDimensional Infiltration. Soil Science Society of America Journal, 41, 285–294.
Kergoat, L. (1998). A model for hydrological equilibrium of leaf area index on a global scale. Journal of Hydrology, 212-213, 268–286.
Liu, S. (2001). Evaluation of the Liu model for predicting rainfall interception in forests world-wide. Hydrological Processes, 15, 2341–2360.
Mouillot, F., Rambal, S. & Lavorel, S. (2001). A generic process-based SImulator for meditERRanean landscApes (SIERRA): Design and validation exercises. Forest Ecology and Management, 147, 75–97.
Ritchie, J. (1972). Model for predicting evaporation from a row crop with incomplete cover. Water resources research, 8, 1204–1213.
Watanabe, T. & Mizutani, K. (1996). Model study on micrometeorological aspects of rainfall interception over an evergreen broad-leaved. Agricultural and Forest Meteorology, 80, 195–214.