Modelling Water Dynamics, Transport Processes and Biogeochemical Reactions in Soil Vadose Zone

Abstract

Large numbers of numerical models are nowadays available for the description of physical and chemical processes affecting water flow and solute transport in soil vadose zone. This chapter explains basic principles of water flow and solute transport modelling in soil vadose (variably saturated) zone and some of the most important processes present in it. First part deals with water dynamics in the soil, that is, soil water content, pressure head, soil porosity, and water flow. Also, some of the measurement techniques used to estimate water dynamics in soil are explained. Water retention curve and soil hydraulic properties needed for modelling are briefly discussed with the explanation of basic (i.e. most commonly used) hydraulic relationship in soil (van Genuchten equation) and water flow (Richards equation) approaches. Second part includes solute transport description in vadose zone, including processes such as advection, diffusion, dispersion, and adsorption. Basic advection‐dispersion equation is explained and also the implementation of boundary and initial conditions in the numerical model. Preferential flow is shortly discussed with the basic principles behind its occurrence and modelling in the soil vadose zone. One real case one‐dimensional (1D) example of modelling with HYDRUS software is presented in which water flow and nitrate transport is simulated on the lysimeter study. Short overview of the most widely used numerical models for simulating vadose zone processes is also presented, whereas the final part is focused on chemical speciation modelling in relatively homogeneous soil solutions using visual MINTEQ interface.

1. Introduction

In the past few decades, more attention is devoted to environmental protection and pollution control in the soil‐plant‐atmosphere continuum. There are growing numbers of potentially harmful substances that are introduced in the environment, especially by agricultural practices and industry in arable and urban areas. Such substances can accumulate in the soil (and potentially get into food chain), surface water, or can be transported in a deeper soil zone and eventually reach groundwater. Intensive agriculture uses large amounts of fertilizer (organic and inorganic), plant protection products (pesticides), animal hormones, and include various other substances that may accumulate in environment (e.g. pathogens, bacteria, and trace metals). Although most of the substances applied in agricultural production are useful in the surface layer of the soil, due to their leaching into the deeper soil layers and groundwater and their bioaccumulation potential, they can also cause serious pollution and degrade natural resources.

Water flow and solute (pollutant) transport models can be used as tools for describing and predicting specific processes in variably saturated soil zone or vadose zone. For example, different models can be used for predicting and/or management of irrigation and drainage systems, crop growth, fertilizer application, and pesticide leaching to protect soil and water resources. Models are also equally necessary for the design of waste disposal sites (industrial, municipal) or long‐term management of various harmful substances (e.g. radioactive waste). A large number of models were developed to simulate the numerous simple hydraulic or complex biogeochemical processes and may be used for different purposes. Vadose zone is in focus of many research topics due to its complex nature and also the possibility of elimination and remediation of present/introduced contaminants before leaching to underground water resources.

Transfer of solutes in the soil is closely linked with the flow of water through vadose zone, which largely affects the concentration and biochemical reactions of various substances. Solute transport in soil vadose zone is one of the most demanding problems that occur in numerical modelling. It includes the transport of water and solutes, chemical reactions and microbial transformations. With the development of new numerical models, it is possible to describe more complex processes that are occurring in the soil‐plant‐atmosphere continuum. In this chapter, basic soil physical concepts and numerical modelling procedures (with example) are explained. Trace metal behaviour in the ecosystems, their mobility in soil, and metal chemical forms (species) in the soil solution are shortly discussed with a chemical speciation modelling example presented.

2. Water dynamics

2.1. Soil water content

Soil system can be defined as a three‐phase system that can be divided into solid, gas, and liquid phase. The fractions of liquid and gas are located in the voids between soil particles. These voids are defined as pores and their quantification can be defined as ratio of pore volume to the total (bulk) volume of a soil. Soil porosity can be estimated using the following expression:

$P=\frac{V_p}{V_s}$E1

where V_p represents volume of pores [L³] and V_s volume of the undisturbed soil sample [L³].

Soil water content represents the quantity of water contained in soil and is expressed as a ratio, which can range from completely dry to the soil porosity value at the point of saturation. Volumetric soil water content is defined as θ_v [L³L^-3]:

${\theta }_v=\frac{V_w}{V_s}$E2

where V_w is the volume of water in soil pores [L³], and V_s is the volume of the undisturbed soil sample [L³]. Its expression is unitless; however, usually, it is expressed as cm³/cm^-3 (or m³/m^-3) to emphasize its volumetric origin. Soil water content can be also defined as a mass, thus named gravimetric water content θ_g [MM^‐1]:

${\theta }_g=\frac{m_w}{m_s}$E3

where m_w is the mass of water and m_s is the mass of soil [M]. It is also unitless but often expressed as gram per gram (or kg per kg), following the same rule as mentioned earlier. Both values can be multiplied by 100 to express it as a percentage.

2.1.1. Soil water content measurements

Measurement of the soil water content can be direct or indirect depending on the used method. Direct measurements include the estimation of water quantity by removal from soil through evaporation, leaching or by chemical procedure. They include destructive soil sampling, and hence, additional soil samples need to be taken in order to achieve more reliable results. Therefore, small undisturbed cores are usually taken to determine water content and bulk density. Indirect methods rely on monitoring soil properties that are directly affected by soil water content (e.g. electrical conductivity). These methods require instrumentation placed in the soil or sensors placed over the soil. For this approach, the calibration is needed in addition to precise installation. The advantage of the in situ water content measurements are repeated measuring at the same location during given time period (e.g. years) without disruption of the soil system. There are various methods used for soil water content determination; here, the most common methods are shortly explained. More in depth explanation of a given methods could be found in Refs. [1, 2].

2.1.1.1. Direct measurement

The gravimetric method for the water content measurement represents simplest way to gather accurate soil water content. A soil sample of a known mass is placed in the container, dried in an oven, removed from the container, and allowed to cool in desiccator, then reweighted. The drying procedure is done by placing the sample in convection oven at 105°C for 24 h. From the measurement, gravimetric water content can be calculated and also converted into volumetric water content using soil bulk density.

2.1.1.2. Indirect measurement

Neutron probe uses radioactive material for measuring soil water content. A neutron meter is placed at the soil surface above the access tube in which probe is lowered into the soil to the desired measurement depth. The probe contains an americium 241/beryllium pellet that emits fast high‐energy neutrons. These high‐speed neutrons pass through the accesses tube and collide with hydrogen atoms in the surrounding soil and water. When neutrons hit H nuclei, they slow down and some are reflected back to the source tube and counted by the neutron detector. Because soil water is the primary source of hydrogen atoms, the count is directly related to the soil water content. This method can produce reliable results and can be used for measurements at multiple depths in few minutes. Its application might be questionable for shallow measurements (>15 cm), since the neutrons might escape from the soil instead of being detected. However, the main disadvantage is the use of radiative source that can be a potential health hazard for the device operator.

Time domain reflectometry (TDR) is a technique that involves measuring the travel time of an electromagnetic wave along a wave guide. The bulk soil properties affecting electromagnetic wave are electrical conductivity and dielectric permittivity. Electrical conductivity is a measure of the free electrons flow when exposed to an electric field. The dielectric permittivity is a measure of the displacement of constrained charges when exposed to an electric field. The speed of the electromagnetic wave in soil is dependent on the dielectric permittivity of soil matrix. The fact that water (80) has much larger dielectric constant than air (1), soil (3–7), or its organic components (2–5) is used to determine the volumetric water content of the soil. The TDR instrument consists of 2–3 parallel rods that are inserted into the soil and act as a waveguides. Electronics in the TDR instrument generate and sense the return of high‐energy signal that travels through the soil along the waveguide, that is, stainless steel rods. The high‐frequency signal is then converted to the volumetric water content. Readings can be affected by high clay content, high organic matter, or high soil salinity. The rods need to be inserted fully into the soil and have good contact with surrounding soil particles. The TDR probes are widely used due to its simplicity, accurate readings, and minimal soil disturbance.

Capacitance devices are used to determine the resonance frequency of a given soil. The capacitor is connected to an oscillator to form an electrical circuit; changes in soil moisture can be detected by changes in the circuit operating frequency. Probes usually consist of two or more electrodes (e.g. plates, rods, or metal rings) that are inserted into the soil. With the ring configuration, the probe is introduced into an access tube installed in the field. When parallel rod configuration is used, the probes are buried at the required depth into the soil and the soil represents a medium between capacitor electrodes. Compared to TDR, frequency domain sensors are relatively inexpensive and have a faster response time. However, because of the complex electrical field around the probe, the sensor needs to be calibrated specifically for different soil types.

2.2. Soil water potential

While the knowledge of soil water content is important, it is also important to know its energy potential since water can be held by the force fields of soils particles differently depending on soil type, (soil texture). For example, two soils with identical water content may have different soil water potential, such that one will easily allow plant water uptake (e.g. sandy soil), while in the other type, soil water will be extracted much harder (e.g. clay soil). The total energy state of soil water is defined by its equivalent potential energy, as determined by the various forces acting on the water per unit quantity. In most cases the kinetic energy of water can be neglected, since the flow rates in soil are very slow. Therefore, the energy state of soil water is defined by its equivalent potential energy, which is derived from its position in a force field. Water in the soil will move from area with high soil water potential energy to area of lower potential energy. Driving force for the flow is the change in potential energy with distance (soil water potential gradient).

These driving forces determine the following:

• Direction and magnitude of water flow

• Plant water uptake

• Drainage amount

• Capillary rise

• Soil temperature changes

• Solute transport.

The potential energy of water in the soil is defined relative to its reference or standard state. Standard state water is pure (no solutes), free (no external forces other than gravity) water at a reference pressure (atmospheric), reference temperature, and reference elevation [3]. Soil water potential is defined as the difference in potential energy per unit of volume, mass, or weight of water compared to the standard (reference state). Depending on the choice of unit for quantity, three different systems can be used (Table 1.)

Name	Symbol	Definition	Dimensions	SI units
Chemical potential	μt	Energy per mass	L^2T^-2	J kg^-1
Soil water potential	ψt	Energy per volume	ML^-1T^-2	N m^-2
Soil water potential head	H	Energy per weight	L	m

Table 1.

Different unit systems for expressing soil water potential (adapted with permission from [4]).

Total soil water potential is defined as the amount of work per unit quantity of pure water that must be completed by external forces to transfer reversibly and isothermally an infinitesimal amount of water from the standard state to the soil at the point under consideration [3]. The transformation of water from the reference states can be divided into the components caused by each force field acting on soil water. These components are forces caused by gravity, hydrostatic pressure, capillarity, solute, air pressure, and swelling [5].

Following are presented the definitions of the most important components of total soil water potential (ψ) which is represented by the sum of its active components:

$\psi ={\psi }_z+ {\psi }_p+{\psi }_s++{\psi }_m+ {\psi }_a$E4

Gravitational potential (ψ_z) is defined as the difference in energy per unit volume or weight between standard water and soil water due to gravity. This component quantifies the effect of the gravitational force field on the energy of soil water.

Hydrostatic pressure potential (ψ_p) is defined as the difference in energy per unit volume or weight between standard water and soil water due to the pressure exerted by overlying free water. This component quantifies the pressure effect from overlying water on the energy of water.

Osmotic (solute) potential (ψ_s) is defined as the difference in energy per unit of volume or weight between standard water and soil water due to the presence of solutes. This component quantifies the effect of solutes on the energy of soil water.

Matric potential (ψ_m) is defined as the difference in energy per unit volume or weight between standard water and soil water due to capillarity and adsorption. This component quantifies the effect of the capillarity and adsorption on the energy of soil water.

Air potential (ψ_a) is defined as the difference in energy per unit volume or weight between standard water and soil water due to effect of soil air pressure. This component quantifies the effect of the air pressure in soil porous system on the energy of soil water.

Some of the components of water potential can be neglected like osmotic pressure and also the effect of air pressure in most of the cases due to its low effect (and estimation difficulty) on the global soil water potential. Following these assumptions, total soil water potential head or hydraulic head is when expressed per unit weight:

$H=h+z$E5

Thus, hydrostatic pressure (h) and gravity (z) dominate the potential energy of water under unsaturated condition.

2.2.1. Measuring soil water potential components

Tensiometer is a measuring device used to determine matric water potential (ψ_m) in the vadose zone. The device consists of an airtight glass or plastic tube filled with water and connected to a porous cup at the bottom. Tensiometers are placed in the soil and the water inside the tube comes into equilibrium with the soil solution (i.e. it is at the same pressure potential as the water held in the soil matrix). Then, the reading is collected from the pressure gauge (water or mercury) at the top of the device. Typically, the measurement range is 0–80 centibars. These devices are easy to use and inexpensive. They require a close contact with surrounding soil around porous cup that might sometimes be hard to achieve, for example, on swelling or coarse soil types. More details about this method can be found in [6].

Piezometers are used to measure hydrostatic potential (ψ_p) and positive pressure head (h), since they are used below water table. A piezometer is a hollow plastic tube installed in the soil, which is open to the atmosphere at the top and located in the saturated soil at the bottom. The bottom of the tube has perforated screened section which allows water to enter the piezometer. Water rises in the tube until the hydrostatic pressure of the water inside the tube is the same as hydrostatic pressure of surrounding soil water. Hydrostatic pressure can be calculated from the water level readings in the piezometer.

Thermocouple psychrometers are used to measure soil matric potential based on the relative humidity of the water vapour in the soil [7]. At equilibrium, the vapour water potential is equal to the liquid water potential. Since the vapour and liquid are at the same elevation, the components of the soil water potential measured by the psychrometer are the sum of the matric and osmotic potential, assuming atmospheric air pressure. The unit consists of a measuring device and soil sensor, which is buried into the soil (ceramic cup) and connected via cable to a measuring device. Electrical circuit (thermocouple) is used to measure temperature. The output of the psychrometer is expressed in voltage, which represents the difference in temperature measured in the ceramic cup and the reference (constant) temperature. If the soil is dry, the output voltage will be greater. A calibration equation is used to convert readings to water potential.

Electrical resistance sensors measure the electrical resistance of a porous block that is in contact with surrounding soil. Electrical resistance between electrodes embedded in a porous medium (block) is proportional to its water content, which is related to the soil water matric potential of the surrounding soil. Electrical resistance decreases as the soil and the block lose water. The most common sensor is a gypsum block. The block is part of a simple DC circuit and buried in the soil. Cable is connecting the block and a voltage measuring device. These sensors are not very accurate [8] and are best suited to manage irrigation systems if precise measurements are not needed.

Soil water constants are used to describe water content across different water potential range in soil and are related to the energy required to extract water from soil (Figure 1).

Maximal water capacity (SWMAX) is the maximal water content of soil, that is, at (or near) saturation. The potential energy gradient is downward through the soil profile, mainly due to gravity forces and through macropores.

Field capacity (FC) is the amount of water that remains 2–3 days after the saturation of a soil with water after gravity movement of water has largely ceased. The water is held in the soil at tension of -0.33 bar (pF 2.0) by matric forces (in micro‐ and mesopores). Water held between saturation and field capacity is subject to free drainage over short time periods, and it is generally considered unavailable to plants.

Permanent wilting point (PWP) is the lowest amount of water in the soil at which matric forces hold water too tight for plant extraction (-15 bar or pF 4.2).

Available field capacity (aFC) or plant available water (PAW) is the difference between field capacity (FC) and wilting point of a soil. It is considered as water available for plants to extract from the soil moisture zone. Plant available water is mostly located in soil micro‐ and mesopores.

Figure 1.

Scheme of the mayor soil water constants: maximal water capacity, field capacity, permanent wilting point and plant available water depending on the soil water potential range.

2.3. Soil water retention curve

Soil water retention curve represents the relationship between the water content (θ), and the soil water potential (h). In the literature, different names could be found such as soil water characteristic curve, capillary pressure saturation relationship, and/or pF curve. The water retention curve provides information on how tight water is held in soil porous system and how much energy would need to extract it from the different pores. The main characteristics of soil water retention curves are visible from Figure 2 where the x‐axis shows the relative water content in the soil (θ), while the y‐axis shows pressure head (h). If the pressure head values are close to 0, the soil is almost completely saturated, as θ decreases, the binding force is getting stronger (more energy is needed to extract water from the soil). At the low‐pressure head values (close to the border wilting point pF 4.2 or -15,000 cm), the water that is retained in the soil is located in the smallest pores.

Figure 2.

Water retention curve example for loam, sand, and clay texture soil based on soil hydraulic parameters taken from [9].

Coarser textured soils (sandy) lose water more quickly than fine textured soils (clay) as a direct reflection of the size distribution of pores in the soil. As most of the pores in the coarser soils have greater diameter, water will percolate during small negative soil water potential, while in the finer textured soils (clay, loam, silty loam), water drainage occurs at very high values of negative soil water potential.

Hysteresis in soil is defined as the difference in the relationship between the water content of the soil and the corresponding water potential obtained under wetting and drying process. Soil water retention curve is usually developed by going from high to low water content producing drying curve. However, if the measurements started from saturated conditions (e.g. pressure plate), this would produce wetting curve. It means that water content in the drying (or drainage) curve of water potential is larger than water content in the wetting curve for the same value of water potential.

Most common methods for soil water retention curve estimation are pressure plate and pressure cells, although there are other methods to determine its shape like hanging water column, suction tables, or soil freezing.

Pressure plate has a range for measuring water retention curve (h=0 to -15,330 cm) and usually can provide very accurate measurements in the wet range. The pressure plate was introduced in the 1930s by Richards [10]. Device is used to make indirect measurements of soil pressure head (matric potential) by imposing a known pressure potential on saturated soil samples until free water no longer flows from the system. When the sample comes to equilibrium, its water potential will be equivalent to the applied pressure.

Pressure chambers are usually used to estimate the dryer part of soil water retention curve. Pressure chambers that hold a single intact soil sample are called pressure (tempe) cells and are usually used for h range from 0 to -1000 or even -3000 cm. The cell consists of plastic housing, a porous ceramic plate at the bottom, a metal ring to secure the soil sample and a rubber sealing between the housing and the ring. The positive known pressure is applied to the sample to extract water and the sample is weighted at specific pressure values (when equilibrium is reached).

The shape of water retention curve can be explained using a wide range of mathematical expressions [11–13]. However, most common expression used is the van Genuchten‐Mualem one [14]:

$\theta \left(h\right)= {\theta }_r+ \frac{{\theta }_s- {\theta }_r}{{(1+ {\left|\alpha h\right|}^n)}^m}\text{for h}< 0$E6

$\theta \left(h\right)= {\theta }_s\text{ for h }\ge 0$E7

$K\left(h\right)= K_s{S}_e^l{(1-{(1-S_e^{\frac{1}{m}})}^m)}^2$E8

$S_e= \frac{\theta -{\theta }_r}{{\theta }_{s- }{\theta }_r}$E9

$m=1- \frac{1}{n}; n\text{>1}$E10

where θ_r and θ_s denote residual and saturated volumetric water content [L³L^-3], respectively, K_s is the saturated hydraulic conductivity [LT^-1], S_e is the effective saturation, α [L^-1] and n [-] are the shape parameters, and l [-] is a pore connectivity parameter. The pore connectivity parameter value is usually taken from an average for many soils (l=0.5) [15].

The hydraulic conductivity characterizes the ability of a soil to transmit water and as such is inversely related to the resistance to water flow [16]. The hydraulic conductivity decreases as soil becomes unsaturated and smaller quantity of pore space is filled with water. The unsaturated hydraulic conductivity function shows the dependency of the hydraulic conductivity on the various ranges of water content or pressure head. Hydraulic conductivity of saturated soil is much higher in coarser texture soils (sand) compared with clay or loam texture soils.

The fitting of soil water retention curve data can be obtained by optimizing some of the parameters used in the equation. This requires non‐linear optimization method. The RETC (RETention Curve) program [17] may be used to predict the hydraulic conductivity from observed soil water retention data assuming that one observed conductivity value (not necessarily at saturation) is available. The program also permits one to fit analytical functions simultaneously to observed water retention and hydraulic conductivity data.

2.4. Soil hydraulic properties estimation

Modelling of water flow and solute transport in soils is predefined with hydraulic properties of studied soil system. Soil hydraulic properties can be obtained directly by conducting laboratory or field measurements and experiments. However, these techniques can be costly and time‐consuming which have resulted in different method for the estimation of soil hydraulic properties. It is possible to estimate these properties from soil information, which are widely available or easy to measure. Mathematical functions to estimate soil hydraulic parameters from basic soil information, such as soil texture, soil organic matter content, bulk density, etc. are called pedotransfer functions (PTFs). Large databases, including soil hydraulic properties and corresponding textures, bulk densities, organic matter content, field capacity, wilting point etc., are available for different soil types. One example is program ROSETTA [18] that predicts parameters of van Genuchten‐Mualem functions. Soil hydraulic parameters are used in agricultural, environmental, and hydrological modelling in which the number of parameters needed depends on simulated processes and used model. Some models require knowledge of the shape of soil water retention curve (SWAP, HYDRUS). Other models, for example, the SWAT model, require water retention values at given matric potentials. Direct point predictions for given matric potentials can lead to more accurate estimations than if water retention values of these matric potentials are derived from predicted SWRC (parameter estimation). Therefore, it is important to have both point predictions and parameter estimations. It is also well known that the performance of prediction models is highly dependent on the quality of the data set (number and type of measured properties, sample size and its heterogeneity) used for their development.

3. Vadose zone modelling

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‐dimensional (1D) water flow can be written as:

$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}\left(K(h)\frac{\partial h}{\partial z}\right)+\frac{\partial K(h)}{\partial z}-S(h)$E11

where θ represents volumetric water content [L³ 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³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, K_s. 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\left(h\right)=\alpha (h)S_p$E12

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 h₁). For a pressure head less than h₄ (the willing point pressure head), water uptake is also assumed to be zero. Water uptake is optimal between pressure heads of h₂ and h₃. The pressure head h₃ may be adjusted depending on the transpiration rate so that it is more negative when transpiration rates are low (optimal root water uptake occurs over a wider range in pressure head at lower transpiration rates).

3.2. Solute transport modelling

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 simultaneously 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=qc$E13

where q_a 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 = D_n + D_m where D is the hydrodynamic dispersion coefficient, D_n is mechanical dispersion, and D_m is molecular diffusion. This can be expressed by Fick's law:

$q_d=-\theta D\frac{\partial c}{\partial z}$E14

where q_d is solute flux density due to hydrodynamic dispersion [ML^-2T^-1], D is dispersion coefficient [L²T^-1], c is solute concentration [ML^-3], θ is soil water content [L³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 (longitudinal) 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).

$\theta D=D_L\left|q\right|+\theta D_m\tau$E15

where D_L is longitudinal dispersivity [L], q is water flux [LT^-1], D_m is molecular diffusion [L²T^-1], θ is soil water content [L³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}$E16

where θ is soil water content [L³L^-3], c is solute concentration [ML^-3], q_d is solute flux density due to hydrodynamic dispersion [ML^-2T^-1], q_a 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 c)}{\partial t}+\frac{\partial (\rho s)}{\partial t}=\frac{\partial }{\partial z}\left(\theta D\frac{\partial c}{\partial z}-qc\right)-\varnothing$E17

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

3.2.1. Adsorption (linear, nonlinear)

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:

$s=K_{d}c$E18

where K_d is the distribution coefficient [L³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:

$s=K_f{c}^{\beta }$E19

$s= \frac{K_{d}c}{1+\eta c}$E20

where K_f [M^-βL^{-3 β}] and β [-] are coefficients used in the Freundlich isotherm, and η [L³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 distribution) 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 (q0=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

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‐permeability 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 permeability models (adapted with the permission from [16]).

3.5. Available vadose zone models

A large number of programs are available for the description of water flow and solute transport in soil vadose zone. Below is the list and brief description of the most widely used software for modelling water flow and solute transport, which are applicable to unsaturated and partially saturated conditions.

HYDRUS‐1D software includes one‐dimensional finite‐element model for simulating the movement of water, heat, and multiple solutes in variably saturated media [37].

HYDRUS (2D/3D) is a software package for simulating water, heat, and solute movement in two‐ and three‐dimensional variably saturated media [38].

MACRO is a one‐dimensional, process‐oriented, dual‐permeability model for water flow and reactive solute transport in soil [30].

TOUGH (‘transport of unsaturated groundwater and heat') suite of software codes are multidimensional numerical models for simulating the coupled transport of water, vapour, non‐condensable gas, and heat in porous and fractured media [39].

SWAP (soil, water, atmosphere, and plant) simulates transport of water, solutes and heat in unsaturated/saturated soils at field scale level, during growing seasons and for long‐term time series [40].

The RZWQM is an integrated physical, biological, and chemical process model that simulates plant growth and movement of water, nutrients, and pesticides in run‐off and percolate within agricultural management systems [41].

Animo is a detailed process‐oriented simulation model for the evaluation of nitrate leaching to groundwater, N and P loads on surface waters and greenhouse gas emission [42].

LEACHM (the leaching estimation and chemistry model) refers to a suite of simulation models describing the water and chemical regime in unsaturated or partially saturated soil profiles [43].

Daisy is a dynamic model for the simulation of water and nitrogen dynamics and crop growth in agro‐ecosystems [44].

4. Water and solute transport using HYDRUS software: case study in Croatia

The example will show the procedure of modelling water flow and nitrogen species transport in soil vadose zone. In this example, urea and NPK fertilizers were added to soil profile. The study site is located in the eastern Croatia in the intensive agricultural production area. The more information about the study can be found in [45]. The project main goal was to evaluate the influence of high fertilizer load, mostly nitrogen based, on the groundwater resources. Here, the results of modelling study on one selected site will be presented from year 2014. HYDRUS‐1D [37] was used to simulate water flow and nitrate transport and its ability to reproduce observed water and nitrate outflows (collected by lysimeters) was assessed. Simulations are carried out in one‐dimensional domain using measured soil hydraulic parameters, climatic daily data, crop growth parameters, and fertilizer application rates as the input for model. Evapotranspiration rates were calculated using CROPWAT model [46] (based on the Penman‐Monteith approach). The soil type at the site was classified as Haplic Gleysol Calcaric Eutric Siltic (Horizons: Ap‐Bg‐Cr‐Cg), from which saturated hydraulic conductivity, retention curve (at different pressure head values), and saturated water content were measured. Optimization of necessary remaining hydraulic parameters (θr, α, and n) was performed using RETC software [17] by fitting the measured data. Figure 7 shows the graph of the water retention curve data (circle) and the fitted water retention equation (line) for two upper soil layers. The fitting of the van Genuchten equation to data was confirmed by high R² values (>0.97).

Figure 7.

Soil water retention curves for (a) 0–25 soil layer, and (b) 25–50 cm soil layer derived from RETC software.

After acquiring all the necessary input data, we can proceed to modelling with HYDRUS. In the main processes window (Figure 8 left), water flow and solute transport are selected with additional root water uptake and root growth option since Barley (Hordeum vulgare L.) was grown at the site. After selecting geometry (length of the profile i.e. 50 cm) and time information (duration of the simulation, that is, 365 days) hydraulic parameters (θ_r, θ_s, α, n, K_s, l) for selected model are inserted (Figure 8 right).

Figure 8.

Main modelling process selection in HYDRUS‐1D (left) and input of soil hydraulic parameters (right).

In the soil hydraulic model window (Figure 9 left), there are multiple option to choose from, varying from single porosity to dual porosity/permeability model with different approaches. In this example, van Genuchten‐Mualem equation was used without considering the effect of soil hysteresis. Initial and boundary conditions were chosen. Atmospheric boundary condition (which includes precipitation and evapotranspiration) was selected at the top and seepage face at the bottom boundary (to mimic lysimeter plate). Root water uptake was simulated using Feddes [20] model. As for water flow, it is necessary to select solute transport parameters (Figure 9 right). The number of solute are set to three since the simulation consider nitrification chain (urea > ammonium > nitrate). The first‐order reaction term representing nitrification of urea to ammonium was 0.38 per day [47]. The first‐order reaction term representing nitrification of ammonium to nitrate was 0.2 per day [47]. The first‐order reaction term for the volatilization of ammonium to ammonia was 0.0552 per day [48]. The distribution coefficient for ammonium (K_d) is assumed to be 3.5 cm³ g^-1 [47]. The unrestricted passive root uptake of urea, ammonium and nitrate was assumed. For solute boundary condition, concentration flux and zero concentration gradient were selected for upper and lower boundary condition, respectively.

Figure 9.

Snapshot of soil hydraulic (left) and solute transport (right) selection window in HYDRUS‐1D.

After inserting all the necessary data (derived from field observation and laboratory measurements of addition optimization) and running the model, pre‐processing (input parameters, Figure 10 left) and post‐processing (results, Figure 10 right) window are displayed in the HYDRUS.

Figure 10.

Snapshot of pre‐processing (input data) and post‐processing (results) window in HYDRUS‐1D.

Cumulative water outflow, which was measured in lysimeter at site 4 during 2014 and simulated ones using HYDRUS‐1D are presented on the Figure 11. The amount of water outflow from lysimeter was mainly the result of high precipitation events. The largest outflows were collected during wet part of the year, while during main crop vegetation season the outflows were very small due to large crop water demand (large transpiration intensity). Very high value of R² (0.96) between the measured and simulated values indicates that the HYDRUS model was capable to reproduce field data with high efficiency. Simulated values of nitrate outflow reflect the water flow pattern (Figure 12) which was measured in the same lysimeter (site 4) during 2014. Although the model derived larger cumulative nitrate outflows compared to the measured ones, the R² value of 0.72 indicate good model ability to simulate nitrogen transformation (from urea and ammonium to nitrate) and transport. The larger simulated nitrate values could be due to denitrification process that might occur in soil which was not considered in the modelling. From Figure 12, it can be seen that the main nitrate leaching occurs at the end of the simulation period after barley harvesting (day 164) and during wet period (autumn/winter). Our results indicate that numerical models can be very helpful for estimating water flow and nitrate dynamics under field conditions. From this example, it can be seen the possible negative influence of the nitrogen fertilizer application and their potential of leaching below root zone. One of the possible numerical models usages is their application in crop water demand (irrigation) and fertilization optimization or pesticide management in agriculture which can eventually lead to the protection of environment.

Figure 11.

Observed cumulative water outflows from lysimeter (site 4) in 2014 and simulated ones using HYDRUS‐1D.

Figure 12.

Observed cumulative nitrate outflows from lysimeter (site 4) in 2014 and simulated ones using HYDRUS‐1D.

5. Trace metals mobility in soil‐plant system

In some recent reports by European Commission [49], contamination by trace metal elements (Cd, Cu, Pb, Zn, Ni, Cr) and some nutrients (N, P) was defined as one of the main pressures to environmental resources in Europe. Contamination of arable soil resources by trace metal elements due to different anthropogenic activities is increasing rapidly and continuously in the last decades at the global level and often with detrimental ecological scenarios [50]. Trace metal elements represent relatively wide group of essential and some of the most toxic elements for biota which occur generally at very low levels in the ‘clean’ pedospheres. For instance, the total contents of trace elements in non‐contaminated mineral top soils range from 1 to 100 ppm [51], but in contaminated soils, their concentrations may be higher by up to several orders of magnitude.

Trace metal behaviour in the ecosystems, their chemical forms (species), mobility and risk of inclusion into the food chain, greatly depend on the environmental conditions. For example, the contamination of soil by trace metals is followed by a cascade of reactions with soil surfaces and the concentration of metals in the soil solution is controlled by a number of inter‐related processes: oxidation/reduction, precipitation/dissolution, adsorption/desorption, inorganic and organic complex formation [52]. Which processes will predominate depends on the physical, chemical, and biological properties of the (non)contaminated soil, as well as the important environmental factors (e.g. moisture, temperature, aeration). The most dominant parameters controlling soil trace metals chemistry and availability to plants are pH and soil OM. At lower soil pH values, hydrogen ions are adsorbed to soil particles, which increases the positive charge on inorganic and organic soil components, resulting in weaker adsorption of metal cations and their increased mobility in soil. SOM has a two‐sided role in metal mobility in soil: (i) particulate SOM is retaining trace metals through the formation of metal‐SOM complexes, thus decreasing their mobility [53], and (ii) dissolved organic matter (DOM) can increase trace metal mobility due to formation of metal‐DOM complexes [54] that are substantially less bounded to soil particles than a free metal ions.

Particular metal forms in soil fractions (soluble, solid, liquid) are possible to detect by adequate analytical technique [55] or predict by some of computational models (e.g. next section). However, inside the cultivated soil (vs. soil without plants), even at very small scales (up to several mm), the physically‐chemical and/or biological characteristics may differ substantially, and thus consequently influence trace metal biogeochemistry. It is especially pronounced at the soil‐root interface, or so‐called rhizosphere microarea. From the biogeochemical perspective, the rhizosphere is very dynamic and heterogenic microarea dominantly controlled by plant roots and released different organic metabolites.

Inside of relatively numerous group of trace metals, during the last decades, copper (Cu) has been intensively studied from different scientific perspectives given on several next facts: it is essential phytonutrient at small concentrations but easily becomes phytotoxic at higher levels; its biochemistry is pH‐dependent; it is one of the most reactive trace metal elements with soil OM and over the huge reactive (mostly negatively charged surface) interface of OM, Cu strongly competes with some other positively charged metals (e.g. Cd, Zn) and thus impacts on its bioavailability. In the next section, one of the most advanced biogeochemical modelling approaches elaborates the influence of particular soil OM substances on Cu chemical speciation in relatively homogenous and realistic rhizosphere conditions.

5.1. Visual MinteQ example: copper chemical speciation

Chemical speciation modelling was recognized as a very useful tool for studying various types of elements (nutrients, trace metals) and their mobility in natural systems such as rhizosphere. In the next, modelling in Visual MINTEQ chemical equilibrium program [56]) was performed for three levels of Cu (non‐contaminated and contaminated system), at three pH levels (covering the most naturally‐occurring pH reactions in different rhizosphere environments) and at three levels of soil dissolved OM (DOC; corresponds for mineral to organic soil types) (Table 2).

		DOC10			DOC20			DOC30
	Species^1^,^2	Cu40	Cu250	Cu500	Cu40	Cu250	Cu500	Cu40	Cu250	Cu500
		% of total concentration			% of total concentration			% of total concentration
pH 5	Cu^2+	4	18	28	1	6	11	0	3	5
	CuNO3^+	0	1	2	0	0	1	0	0	0
	CuCO3 (aq)	0	0	0	0	0	0	0	0	0
	Cu(OH)2 (aq)	0	0	0	0	0	0	0	0	0
	Cu(CO3)2^2-	0	0	0	0	0	0	0	0	0
	HA1‐Cu (6) (aq)	60	61	55	55	67	68	50	66	70
	HA2‐Cu (6) (aq)	36	20	14	44	27	21	49	31	25
pH 7	Cu^2+	0	1	2	0	0	0	0	0	0
	CuNO3^+	0	0	0	0	0	0	0	0	0
	CuCO3 (aq)	0	1	4	0	0	1	0	0	0
	Cu(OH)2 (aq)	0	0	0	0	0	0	0	0	0
	Cu(CO3)2^2-	0	0	0	0	0	0	0	0	0
	HA1‐Cu (6) (aq)	5	13	18	3	9	13	3	7	11
	HA2‐Cu (6) (aq)	95	85	75	97	91	85	97	93	89
pH 9	Cu^2+	0	0	0	0	0	0	0	0	0
	CuNO3^+	0	0	0	0	0	0	0	0	0
	CuCO3 (aq)	0	1	4	0	0	0	0	0	0
	Cu(OH)2 (aq)	0	0	1	0	0	0	0	0	0
	Cu(CO3)2^2-	0	0	2	0	0	0	0	0	0
	HA1‐Cu (6) (aq)	0	1	2	0	1	1	0	0	1
	HA2‐Cu (6) (aq)	100	97	91	100	99	98	100	99	99

Table 2.

Distribution (%) of Cu species in tested soil solution, estimated by Visual MINTEQ chemical equilibrium software (NICA‐Donnan model) as affected by soil pH (5, 7, and 9), dissolved organic carbon (DOC; 10, 20, and 30 mgL^-1) and different soil Cu total concentration (40, 250, and 500 mg kg^-1).

¹HA‐Cu: humic acid‐complexed Cu via 1‐carboxylic and 2‐phenolic functional groups.

²Species with the < 0.5% of total concentration are not shown.

Even though the total soil Cu content by itself is not an adequate measure to determine Cu mobility and phytoavailability, a strong positive correlation between total element concentration and its bioavailable fraction is often reported, especially in contaminated soils or in trials on soils spiked by metals [57]. In this model, chosen soil total Cu concentrations were (in mg kg^‐-1); 40 (correspond to most of non‐contaminated soil conditions), 250 and 500 (corresponds from medium to highly contaminated soil conditions). The mobility of trace metals in soil depends ultimately on their chemical speciation, which is actually a function of pH and the presence of inorganic and organic ligands in the soil solution. The non‐ideal competitive adsorption (NICA)‐Donnan (sub)model for the adsorption of cations onto dissolved organic matter (DOC) and software database on equilibrium constants [58] was applied in this model. Cu concentrations used for Cu speciation here were average values obtained from experimental data from CaCl₂ extracts, and for other elements in soil solution from saturated soil water extracts [53]. Ionic strength was set to be calculated and temperature of 25°C for all the calculations. Other software default settings were not modified.

From the data presented in Table 2, it can be seen that the majority of Cu in the modelled solutions was founded as a Cu‐DOC complex, suggesting that DOC might be the main factor affecting Cu ultimate fate after its release from the soil solid phase into the rhizosphere solution. Increase in total Cu concentration led to an increase in free Cu²⁺ ion in the modelled solutions, but only at lower concentrations of DOC and at acid pH (Table 2). Even though in every investigated scenario, the majority of Cu in the modelled solution was found to be complexed with DOC (Table 2), an increased soil pH caused the reduction of the percentage of free Cu ion in the solution. It is known that free metal ion is considered the most mobile species in the soil, thus data are implicating a decreased metal mobility in soil with an increase in soil pH. Furthermore, with an increase in soil pH Cu shifted from carboxylic groups in humic acids (HA1‐Cu (6)) to phenolic functional groups (HA2‐Cu (6)), which is in agreement with some previous models that phenolic groups might be more important for metal (e.g. Cd, Cu) binding under higher pH of surrounding media [55].

6. Conclusion

Located at the atmosphere, plant, soil, and water interface, vadose zone is represented by a variety of linked complex processes. For solving problems such as the transport of nutrients, pesticides, pharmaceuticals, colloids, bacteria, viruses, hormones, and toxic trace elements, carbon sequestration, and bioremediation of organic contaminants, a thorough understanding and coupling of multiple hydrogeological, geochemical, and microbiological processes is needed. Models should be considered as one of the most advanced and useful ‘tools’ which, when used properly, can predict different scenarios, with positive or negative outcome, that can occur in the natural systems. With development of numerical models, such complex problems can be solved more successfully using different mathematical expressions and approaches. The accuracy of model predictions rely largely upon a quality and quantity of input parameters required for a specific problems, mostly due to large heterogeneity of the soil. Thus, fundamental knowledge of basic soil physics is needed, which combined with new measurement techniques provides satisfactory foundation for performing modelling. In recent years, scientists have been mostly engaged in the coupling of different numerical models, since no single model is yet available for describing such complex system as soil vadose zone [59]. Development of coupled numerical models capable of describing unstable preferential flow in soils, as well as models coupled with sophisticated geochemical models capable of describing complex kinetic chemical and biological reactions will remain a focus of research in the near future.