**3.1 Water flow modelling**

Water flow modelling in the unsaturated zone is based on Richards equation [19]. It is a non‐ linear partial differential equation, which is often difficult to approximate since it does not have a closed‐form analytical solution. The expression of Richards equation for one‐dimen‐ sional (1D) water flow can be written as:

$$\frac{\partial \mathcal{O}}{\partial t} = \frac{\partial}{\partial z} \left( K(h) \frac{\partial h}{\partial z} \right) + \frac{\partial K(h)}{\partial z} - S(h) \tag{11}$$

where *θ* represents volumetric water content [L<sup>3</sup> L-3], *h* is pressure head [L], *z* is vertical coordinate (positive upwards) [L], *t* is time [T], *K* is the unsaturated hydraulic conductivity [LT-1], *S* represents a sink term (root water uptake) [L<sup>3</sup> L-3T-1]. Richards equation is based on continuity equation and Darcy's law. The continuity equation in general states that the change in the water content (storage) in a given volume is due to spatial changes in the water flux. Darcy's law is a constitutive equation that describes the flow of a fluid through a porous medium. The Darcy‐Buckingham equation is formally similar to the Darcy's equation, except that the proportionality constant (i.e. the unsaturated hydraulic conductivity) in the Darcy‐ Buckingham equation is a nonlinear function of the pressure head (or water saturation), while *K(h)* in Darcy's equation is a constant equal to the saturated hydraulic conductivity, *Ks*. Root water uptake can be also considered through a sink term in the Richards equation. Sink term is usually modelled using Feddes equation [20]:

$$S(h) = \alpha(h)S\_p \tag{12}$$

where *S(h)* is the root uptake removed from a unit of soil in time [T-1], α(*h*) is water stress response function, which varies between 0 and 1 [T-1], and *Sp* is potential root water uptake rate [T-1].

**Figure 3.** Plant water stress response function according to Feddes equation.

The stress response function is shown in **Figure 3**. Water uptake is assumed to be 0 close to saturation due to a lack of oxygen in the root zone (pressure head greater than h1). For a pressure head less than *h*<sup>4</sup> (the willing point pressure head), water uptake is also assumed to be zero. Water uptake is optimal between pressure heads of *h*<sup>2</sup> and *h*3. The pressure head *h*<sup>3</sup> may be adjusted depending on the transpiration rate so that it is more negative when transpi‐ ration rates are low (optimal root water uptake occurs over a wider range in pressure head at lower transpiration rates).

#### **3.2. Solute transport modelling**

have a closed‐form analytical solution. The expression of Richards equation for one‐dimen‐

( ) ( ) ( ) *h Kh K h S h*

where *θ* represents volumetric water content [L<sup>3</sup> L-3], *h* is pressure head [L], *z* is vertical coordinate (positive upwards) [L], *t* is time [T], *K* is the unsaturated hydraulic conductivity

continuity equation and Darcy's law. The continuity equation in general states that the change in the water content (storage) in a given volume is due to spatial changes in the water flux. Darcy's law is a constitutive equation that describes the flow of a fluid through a porous medium. The Darcy‐Buckingham equation is formally similar to the Darcy's equation, except that the proportionality constant (i.e. the unsaturated hydraulic conductivity) in the Darcy‐ Buckingham equation is a nonlinear function of the pressure head (or water saturation), while *K(h)* in Darcy's equation is a constant equal to the saturated hydraulic conductivity, *Ks*. Root water uptake can be also considered through a sink term in the Richards equation. Sink term

> ( ) ( ) *<sup>p</sup> S h hS* = a

where *S(h)* is the root uptake removed from a unit of soil in time [T-1], α(*h*) is water stress response function, which varies between 0 and 1 [T-1], and *Sp* is potential root water uptake

The stress response function is shown in **Figure 3**. Water uptake is assumed to be 0 close to saturation due to a lack of oxygen in the root zone (pressure head greater than h1). For a

æ ö = +- ç ÷ ¶¶ ¶ ¶ è ø (11)

L-3T-1]. Richards equation is based on

(12)

*tz z z* ¶¶ ¶¶

sional (1D) water flow can be written as:

144 Groundwater - Contaminant and Resource Management

q

[LT-1], *S* represents a sink term (root water uptake) [L<sup>3</sup>

is usually modelled using Feddes equation [20]:

**Figure 3.** Plant water stress response function according to Feddes equation.

rate [T-1].

Transport of various substances in soil is associated with water flow. Some of the substances present in a soil system dissolve in water and they are transported through the soil. On the other hand, some of the substances do not dissolve in water and they are transported simul‐ taneously with water. Substances either do not react with surrounding soil system and do not change during time or react with soil and change due to chemical reactions, microbiological transformations, etc. Given the solute transport is either conservative (transported solute mass is constant) or non‐conservative (transported solute mass is usually decreasing due to adsorption, nitrification, degradation, volatization, etc.). Some of the main processes that affect solute transport in soil are explained in this section.

*Advection* represents transport of solute mass at the average rate caused by the water flux.

$$q\_a = q\mathbf{c} \tag{13}$$

where *qa* is solute flux density due to advection [ML-2T-1], *c* is solute concentration [ML-3], and *q* is water flux density [LT-1]. If the solute transport in the soil would be controlled only by advection its transport velocity would be identical as the average water flow in soil.

*Hydrodynamic dispersion* includes solute spreading by mechanical dispersion and molecular diffusion in direction of the flow (longitudinal dispersion) and perpendicular to it (transverse dispersion). Mechanical dispersion occurs in the soil as a result of differences in the pore size, the difference in the flow path length and ongoing mixing between pores (due to arrangement of pores in the soil) and the difference in transport velocity within a pore [16]. This process happens in the micro (within the pores) and macro (preferential flow through cracks in the soil) scale. Molecular diffusion represents transfer of solutes due to concentration gradient. Hydrodynamic dispersion is the result of two processes *D = Dn + Dm* where *D* is the hydrody‐ namic dispersion coefficient, *Dn* is mechanical dispersion, and *Dm* is molecular diffusion. This can be expressed by Fick's law:

$$q\_d = -\theta D \frac{\partial c}{\partial z} \tag{14}$$

where *qd* is solute flux density due to hydrodynamic dispersion [ML-2T-1], *D* is dispersion coefficient [L2 T-1], *c* is solute concentration [ML-3], *θ* is soil water content [L3 L-3] and *z* represent spatial coordinate [L].

*The advection‐dispersion equation* (ADE) is most common mathematical description used for solute transport in the soil. For numerical or analytical solution of ADE, it is necessary to know the value of the dispersion coefficient. Hydrodynamic dispersion coefficient can be estimated by field or laboratory experiments. The most commonly used method for estimating dispersion coefficient is the application of tracer (a non‐reacting compound) to soil column and estimation of the dispersion coefficient from gathered data. From the experiment, a *breakthrough curve* is derived, which plots relative concentration of a given substance versus time, where relative concentration is defined as the ratio of the actual concentration to the source concentration.

The dispersion coefficient in one‐dimensional (1D) system is proportional to the average water flow velocity in porous system where proportionality constant is referred to as the (longitu‐ dinal) dispersion [21]. Dispersion can be derived from Newton's law of viscosity which states that velocities within a single capillary tube follow a parabolic distribution, with the largest velocity in the middle of the pore and zero velocities at the walls (this can be applied on soil porous system, **Figure 4a**). Solute transport due to dispersion is the result of the unequal distribution of water flow rate in the pores of different sizes. Since soils consist of pores of different radii, solute fluxes will be significantly different, with some solutes again traveling faster than others (**Figure 4b**).

$$
\rho \theta D = D\_L \left| q \right| + \theta D\_n \tau \tag{15}
$$

where *DL* is longitudinal dispersivity [L], *q* is water flux [LT-1], *Dm* is molecular diffusion [L2 T-1], *θ* is soil water content [L3 L-3], and *τ* is tortuosity factor [‐].

**Figure 4.** Distribution of single pore velocity (a), and distribution of velocities in a larger‐complex pore system (b, adapted with the permission from [16]).

Continuity equation is used to calculate the mass balance of solute in soil, that is, it states the changes of total solute concentration in time per volume of soil:

$$\frac{\partial(\theta c)}{\partial t} = -\frac{\partial(q\_d + q\_a)}{\partial z} \tag{16}$$

where *θ* is soil water content [L3 L-3], *c* is solute concentration [ML-3], *qd* is solute flux density due to hydrodynamic dispersion [ML-2T-1], *qa* is solute flux density due to advection [ML-2T-1], and *z* is coordinate [L].

Adsorption dispersion equation is based on the continuity, advection and dispersion equation. The ADE for one‐dimensional solute transport during transient water flow in a variably saturated medium:

$$\frac{\partial(\theta \mathbf{c})}{\partial t} + \frac{\partial(\rho \mathbf{s})}{\partial t} = \frac{\partial}{\partial z} \left( \theta D \frac{\partial \mathbf{c}}{\partial z} - q\mathbf{c} \right) - \mathcal{Q} \tag{17}$$

where *c* is solute concentration [ML-3], *s* is adsorbed concentration [MM-1], *θ* is soil water content [L3 L-3], *ρ* is soil bulk density [ML-3], *D* is dispersion coefficient [L2 T-1], *q* is volumetric flux [LT-1], and *Ø* is rate constant representing reactions [ML-3 T-1].

#### *3.2.1. Adsorption (linear, nonlinear)*

*The advection‐dispersion equation* (ADE) is most common mathematical description used for solute transport in the soil. For numerical or analytical solution of ADE, it is necessary to know the value of the dispersion coefficient. Hydrodynamic dispersion coefficient can be estimated by field or laboratory experiments. The most commonly used method for estimating dispersion coefficient is the application of tracer (a non‐reacting compound) to soil column and estimation of the dispersion coefficient from gathered data. From the experiment, a *breakthrough curve* is derived, which plots relative concentration of a given substance versus time, where relative concentration is defined as the ratio of the actual concentration to the source concentration.

The dispersion coefficient in one‐dimensional (1D) system is proportional to the average water flow velocity in porous system where proportionality constant is referred to as the (longitu‐ dinal) dispersion [21]. Dispersion can be derived from Newton's law of viscosity which states that velocities within a single capillary tube follow a parabolic distribution, with the largest velocity in the middle of the pore and zero velocities at the walls (this can be applied on soil porous system, **Figure 4a**). Solute transport due to dispersion is the result of the unequal distribution of water flow rate in the pores of different sizes. Since soils consist of pores of different radii, solute fluxes will be significantly different, with some solutes again traveling

> qt

where *DL* is longitudinal dispersivity [L], *q* is water flux [LT-1], *Dm* is molecular diffusion

**Figure 4.** Distribution of single pore velocity (a), and distribution of velocities in a larger‐complex pore system (b,

Continuity equation is used to calculate the mass balance of solute in soil, that is, it states the

() ( ) *d a c qq t z* ¶ ¶+

L-3], and *τ* is tortuosity factor [‐].

*D Dq D* = + *L m* (15)

= - ¶ ¶ (16)

q

changes of total solute concentration in time per volume of soil:

q

faster than others (**Figure 4b**).

146 Groundwater - Contaminant and Resource Management

T-1], *θ* is soil water content [L3

adapted with the permission from [16]).

[L2

For solving solute transport equation, additional information is needed to describe the relationships between various substances whose transport is modelled. Adsorption is a physical and chemical process in which one substance is bound to the surface of the other phase (in this case, the binding of substances occur mostly to the soil solid phase). The most widely used and simplest way of describing this process is to assume instantaneous sorption and to use adsorption isotherms. The most elementary form of the adsorption isotherm is the linear isotherm given by the following equation:

$$\mathbf{s} = \mathbf{K}\_d \mathbf{c} \tag{18}$$

where *Kd* is the distribution coefficient [L3 M-1]. This assumption simplifies the mathematical description of solute transport; however, adsorption is generally nonlinear and often depends on the presence of various competing substances in the soil solution. The most commonly applied models used to describe the nonlinear adsorption are presented by [22] Freundlich and [23] Langmuir:

$$\mathbf{s} = \mathbf{K}\_f \mathbf{c}^\emptyset \tag{19}$$

$$\mathbf{s} = \frac{K\_d \mathbf{c}}{\mathbf{l} + \eta \mathbf{c}} \tag{20}$$

where *Kf* [M-βL-3 β] and *β* [-] are coefficients used in the Freundlich isotherm, and *η* [L3 M-1] is coefficient used for description of the Langmuir isotherm. Linear adsorption represents a case where the Freundlich equation has *β* equal to 1. For cases when *β*< 1, less solute mass is adsorbed per unit increase in *c* at high concentrations compared with low concentrations (**Figure 5a**). This happens when the amount of solute added exceeds the ability of adsorption of a type of soil. In the Langmuir model with increasing *η* value, the isotherm becomes more non‐linear and is approaching the maximum sorbed concentration (**Figure 5b**).

**Figure 5.** Examples of nonlinear adsorption isotherm using the Freundlich (a), and Langmuir (b) equations assuming different *β* and *η* values.

#### **3.3. Initial and boundary conditions**

Before starting the simulation process, it is necessary to determine the initial conditions of water flow (the potential distribution at *t*=0) and solute transport (distribution of solute concentration at *t*=0). This means that it is necessary to describe the initial state of the simulated soil system in terms of the relative amounts of water (pressure head or water content distri‐ bution) and concentration of simulated solute or solutes in the soil profile at the moment of the beginning of the simulation. Boundary conditions are the conditions specified at the edges of the transport model area and they define how the site specific model interacts with its environment. The water flow and transport equations can be solved analytically or numerically after determination of the initial and boundary conditions. Complex interactions between the transport domain and surrounding unsaturated zone should be considered carefully for any problem having in mind that this interaction determines the dynamics of water flow and solute transport (velocity and quantity) on the domain boundaries. To solve the solute transport equation, most numerical models used three types of boundary conditions. Dirichlet (first) type of boundary conditions is used when the concentration of the solute at the boundary of domain is known. In some cases, for example, when a boundary is impermeable (*q*0=0) or when water flow is directed out of the region, Neumann (second) type of boundary conditions is used. Cauchy (third) type of boundary conditions can be used to describe solute transport on the domain boundaries (a combination of the first two types of boundary conditions). Since Cauchy boundary conditions define the solute flux across a boundary, the solute flux entering the transport domain will be known exactly. This specified solute flux is then divided into advective and dispersive components. On the other side, Dirichlet boundary condition controls only the concentration on the boundary, but not the solute flux because its advective and dispersive contributions will be larger than for the Cauchy boundary condition [4].

#### **3.4. Preferential flow modelling**

adsorbed per unit increase in *c* at high concentrations compared with low concentrations (**Figure 5a**). This happens when the amount of solute added exceeds the ability of adsorption of a type of soil. In the Langmuir model with increasing *η* value, the isotherm becomes more

**Figure 5.** Examples of nonlinear adsorption isotherm using the Freundlich (a), and Langmuir (b) equations assuming

Before starting the simulation process, it is necessary to determine the initial conditions of water flow (the potential distribution at *t*=0) and solute transport (distribution of solute concentration at *t*=0). This means that it is necessary to describe the initial state of the simulated soil system in terms of the relative amounts of water (pressure head or water content distri‐ bution) and concentration of simulated solute or solutes in the soil profile at the moment of the beginning of the simulation. Boundary conditions are the conditions specified at the edges of the transport model area and they define how the site specific model interacts with its environment. The water flow and transport equations can be solved analytically or numerically after determination of the initial and boundary conditions. Complex interactions between the transport domain and surrounding unsaturated zone should be considered carefully for any problem having in mind that this interaction determines the dynamics of water flow and solute transport (velocity and quantity) on the domain boundaries. To solve the solute transport equation, most numerical models used three types of boundary conditions. Dirichlet (first) type of boundary conditions is used when the concentration of the solute at the boundary of domain is known. In some cases, for example, when a boundary is impermeable (*q*0=0) or when water flow is directed out of the region, Neumann (second) type of boundary conditions is used. Cauchy (third) type of boundary conditions can be used to describe solute transport on the domain boundaries (a combination of the first two types of boundary conditions). Since Cauchy boundary conditions define the solute flux across a boundary, the solute flux entering the transport domain will be known exactly. This specified solute flux is then divided into advective and dispersive components. On the other side, Dirichlet boundary condition controls only the concentration on the boundary, but not the solute flux because its advective and dispersive contributions will be larger than for the Cauchy boundary condition [4].

different *β* and *η* values.

**3.3. Initial and boundary conditions**

148 Groundwater - Contaminant and Resource Management

non‐linear and is approaching the maximum sorbed concentration (**Figure 5b**).

The term preferential flow combines all transport where water and solutes move along certain pathways, while bypassing other volume fractions of the porous soil matrix [24, 25]. In heterogeneous structured soils, which contains large interconnected voids (e.g. root channels, fissures, earthworm pathways) water and transported solutes bypass soil matrix creating non‐ equilibrium conditions in pressure heads and solute concentrations between preferential flow paths and the soil matrix‐pore region. Large number of approaches has been developed to model preferential flow in soil vadose zone. Most of these models try to separately describe flow and solute transport in matrix and fracture pore regions, that is, dual porosity, dual permeability models and multiporosity or multipermeability models [26–29] (**Figure 6**). Dual‐ porosity and dual‐permeability models both assume that the porous medium consists of two interacting regions, one associated with the interaggregate, macropore, or fracture system, and one comprising micropores (or intra‐aggregate pores) inside soil matrix (soil aggregates). While dual‐porosity models assume that water in the soil matrix is stagnant (immobile), dual‐ permeability models allow for water movement in the matrix domain as well. Dual‐permea‐ bility models in which water can move in both the inter‐ and intra‐aggregate pore regions are now also becoming more popular [27, 30]. The main difference between available dual permeability is the implementation of water flow in and between two pore regions (i.e. fracture and matrix). Different approaches are used to estimate water flow in fracture and matrix domains like the Poiseuille's equation [31], the Green and Ampt or Philip infiltration models [32], the kinematic wave equation [30, 33, 34], and the Richards equation [27]. Multiporosity and/or multipermeability models are based on the same concept as dual‐porosity and dual‐ permeability models but include additional interacting pore regions [35, 36]. In these models, the transport of solute mass is determined with the transfer rate which describes the transport of solutes between the fracture and matrix domain by the sum of diffusive and convective fluxes. Straightforward descriptions of main preferential flow modelling approaches are given in reviews by [29] and [25].

**Figure 6.** Scheme of transport processes assumptions in (from left to right) single‐porosity, dual‐porosity and dual per‐ meability models (adapted with the permission from [16]).
