Modelling Flow and Transport in Unsaturated Zone

Abstract

Unsaturated zone transport models are indispensable tools for analyzing complex environmental pollution problems, and for developing practical management strategies. A quantitative study of water flow and contaminant transport in the unsaturated zone is a key factor in the improvement and protection of the quality of groundwater supplies. This is the region bounded above by the land surface and below by the groundwater table. It is the region through which water derived from precipitation and irrigation infiltrates and transports contaminants to reach the groundwater. This article presents an overview of the modelling process for water flow and contaminant transport in the unsaturated zone, input data requirements and related software packages.

Introduction

The unsaturated zone is the region through which water, together with pollutant carried by the water, must pass to reach the groundwater. Therefore various processes occurring within the unsaturated zone play a major role in determining both the quality and quantity of water recharging into the groundwater. A quantitative study of water flow and contaminant transport in the unsaturated zone is a key factor in the improvement and protection of the quality of groundwater supplies.

Numerous simulation models for water flow and solute transportation in the unsaturated region are now being used increasingly for numerous functions in both research and administration. Modelling techniques vary from straightforward analytical or even semi-analytical methods, to sophisticated numerical codes. Although analytical and even semi-analytical methods remain widely used for certain functions, the growing power of PCs along with the progression of more precise and numerically solid outcomes have inspired significantly broader usage of numerical codes in recent years. The extensive utilization of numerical models is additionally greatly improved by their accessibility in both the commercial and public areas, and by the advancement of innovative graphics-based interfaces which could significantly simplify their use.

Mathematical Equations of Water and Transport in Unsaturated Soils

Analytical, semi-analytical, and numerical models are usually based on the following three governing equations for water flow, solute transport, and heat movement, respectively:

Suitable simplifications (mostly for analytical approaches) or extensions thereof (e.g. for two- and three-dimensional systems) are also employed. In equation (1), frequently known as the Richards equation, h is the pressure head, z is the perpendicular coordinate positive upwards, θ is the water content, t is time, S is a sink term representing root water uptake or some other sources or sinks, and the hydraulic conductivity function (unsaturated) K(h) is, often given as the product of the relative hydraulic conductivity, Kr, and the saturated hydraulic conductivity, Ks. In equation (2), called the convection-dispersion equation (CDE), c is the concentration of the solution, D is the dispersion coefficient comprising both hydrodynamic dispersion and molecular diffusion, the retardation variable that represents adsorption (R), the volumetric fluid flux density (q), and Φ is a sink/source term that accounts for various zero- and first order or other reactions. In equation (3), T is temperature, λ is the apparent thermal conductivity, and C and Cw are the volumetric heat capacities of the soil and the liquid phase, respectively.

Solutions of the Richards equation (1) require knowledge of the unsaturated soil hydraulic functions, that is, the soil water retention curve, θ(h), describing the relationship between the water content θ and the pressure head h, and the hydraulic conductivity function (unsaturated),

K(h), defining the hydraulic conductivity K as a function of h or θ. While under certain conditions (i.e. for linear sorption, a concentration-independent sink term Φ, and a steady flow field) equations (2) and (3) are linear equations, equation (1) is often very nonlinear due to the nonlinearity of the soil's hydraulic attributes. Consequently, numerous analytical methods were derived previously for equations (2) and (3) and these kinds of analytical methods today are popular for evaluating solute and heat transport under steady-state conditions. Although a large number of analytical solutions of (1) exist, they can generally be applied only to drastically simplified problems. Most of the functions for water flow inside the vadose zone demands a numerical approach of the Richards formula.

Required Input Data

Simulation of water qualities in the unsaturated areas needs input info regarding the model parameters, the geometry of the system, the boundary conditions and, when simulating transient flow, initial conditions. With geometry parameters, the dimensions of the problem domain are defined. With the physical parameters, the physical properties of the system under consideration are described. With respect to the unsaturated zone, it concerns h(θ) (characteristics of soil water), and K(θ) (hydraulic conductivity).

For an appropriate explanation of the unsaturated flow, a correct description of the two hydraulic functions, K(θ) and h(θ), is important. K(θ), the hydraulic conductivity, declines significantly when the moisture content θ, declines from saturation. The experimental approach to calculating K(θ) at different moisture contents is fairly complicated but not too trustworthy. Additional methods were recommended to determine the K( θ ) function from more conveniently measurable characterizing qualities of the soil. In many studies, the hydraulic conductivity of the unsaturated soil is defined as product of a non-linear function of the effective water saturation, together with hydraulic conductivity at saturation. The relationship is shown by

… (4)

where,

The value of n is found to be 3.5 for coarse textured soils. n will vary with soil type. In literature, established empirical correlation between n and soil characteristic is available.

The relationship between the soil water pressure head h(θ) and moisture content θ, usually termed as the water retention curve or the soil moisture characteristic, is basically determined by the textural and the structural composition of the soil. Also the organic matter content may have an influence on the relationship. A characteristic feature of the water retention curve is that suction head (-h) decreases fairly rapidly with increasing moisture content. Hysteresis results might emerge, and, rather than being a single valued association, the h-θ connection includes a group of curves. The actual curve will have to be determined from the history of wetting and drying.

When root water uptake is additionally modelled, the criteria describing the association between root water uptake together with soil water state has to be supplied, along with plant details. If a functional flux-head association is employed as a lower limit stipulation, the criteria identifying the relationship between surface water and groundwater and, if required, the vertical resilience of inadequately porous layers must be provided.

The number and kind of criteria necessary for modeling flow as well as transport procedures in soils are dependent on the kind of model selected. These parameters can be categorized as control parameters (controlling the operation of the computer code), discretization data (grid and time stepping), and material parameters. The material criteria could be gathered in 6 units (Jury and Valentine, 1986) – static soil properties, water transport and retention attributes, time-dependent criteria, soil adsorption variables,standard chemical characteristics and tortuosity functions. Table 1 lists many of the relevant material model parameters.

Modelling of Unsaturated Flow

Analytical solutions, when available, provide better insight into the physics behind the transport phenomena and are computationally efficient and simple to use. However, analytical approaches are for the most part limited to situations of simple geometry domains, linear governing equations, and homogeneous systems. Along with efficient numerical methods and rapidly updated computer hardware, a large number of numerical models have been developed. However, the numerical technique cannot replace the analytical approaches completely, since numerical methods themselves need verification against analytical solutions because of discretization errors and convergence and stability problems that may be especially troublesome for advection-dominated and nonlinear adsorption processes.

Analytical approaches to the Richards formula for unsaturated flow under certain boundary and preliminary conditions are difficult to obtain because of the nonlinearity in soil hydraulic parameters. This difficulty is exaggerated in the case where the soil is heterogeneous. Usually, you have to depend on numerical methods for forecasting water motion in unsaturated soils, even for soils that are homogeneous. However, numerical approaches often suffer from convergence and mass balance problems.

Analytical solutions are only applicable to highly simplified systems, and are not well suited for situations normally encountered in the field. Initially, finite difference methods were mainly developed to predict unsaturated flow only; but finite element solution techniques also became available later. The nonlinearity of Richards equation is usually solved using an iterative procedure such as Newton or Picard methods. Perhaps the most important advantage of finite element techniques over standard finite difference methods is the ability to describe irregular system boundaries in simulations more accurately, as well as easily including nonhomogeneous medium properties. For one-dimensional simulations, finite difference methods are just as good as finite element schemes. However, several authors suggest that finite element methods lead to more stable and accurate solutions, thus permitting larger time steps and/or coarser grid systems, and hence leading to computationally more efficient numerical schemes.

To numerically solve coupled systems of equations, the solution process requires some manipulation at each time step so that the dependence of one equation on the solution of the other is dealt with accurately. One way to overcome this is to use a fully implicit approach to solve the equations simultaneously. Any nonlinearity of the generated system can be handled by Newton’s method. The implicit nature of this scheme allows for larger time steps in the simulation to find stable solutions as compared to the time steps for explicit schemes. Another method of using the totally implicit scheme is to apply the blended implicit-explicit approach. Yet somehow, the specific portion of the scheme implies that this algorithm is currently dependent upon a stability constraint which significantly limits the measurement of the time step and opens up numerical artifacts.

For heterogeneous soils that contain macropores, a different modelling approach is needed, as the presence of macropores in these soils may form a separate network for water flow. The common approach is to introduce two-region flow domain, one for macropores and the other for the soil matrix. In each flow domain, hydraulic conductivity is given independently. One more strategy would be to think of the heterogeneous soil as a group of stream tubes. It is assumed that there is no exchange of water between these tubes, and that within each tube, the hydraulic conductivity is defined, but varies between tubes.

Modelling of Solute Transport

Movement of dissolved solutes in soils is often defined by the advection-dispersion formula. Analytical answers are produced for various perimeter and preliminary situations. Although these solutions are obtained for narrowly specified circumstances, there are numerous applications like the confirmation of computer codes, forecast of solute motion for large times or distances where the use of numerical models become impractical, and the analysis of transport criteria from soil column research. The majority of analytical solutions pertain to semi-infinite and infinite media. Solutions are obtained by a variety of mathematical techniques including Green’s functions, separation of variables, characteristics method, Laplace transforms and Fourier transforms.

Forecast of solute migration under field situations requires the concurrent solution of the solute transport and unsaturated flow equations. First approximations involve or assume steady flow and constant water contents. Because of the natural complexity of unsaturated flow, strategies for forecasting solute transport have depended mostly on finite difference or finite element estimates of the governing formulas. Since advection-diffusion transport equations normally do not have closed form analytical outcomes, it is important to have accurate numerical approximations. Any time diffusion governs the physical procedure, regular finite difference and element methods function effectively in solving these equations. But yet, if advection is the principal procedure, these equations could display numerous numerical problems. In fact, any standard finite difference or element method will produce a solution, which exhibits non-physical oscillations.

One of the distinctive features of the porous media on the field scale is the spatial heterogeneity of transport properties. These features have a distinct effect on the spatial distribution of contaminant concentration, as has been observed in field experiments and demonstrated by simulation of contaminant transport in unsaturated, heterogeneous soil. Description of the mixing process due to spatial variability of the unsaturated hydraulic conductivity has been advanced with the development of numerical solutions, which assume spatially variable soil properties; stochastic models; and stochastic stream tube models, which decompose the field into a set of independent vertical soil columns.

Unsaturated Zone Modelling Software

Most of the early models developed for studying processes in the near-surface environment focused mainly on variably saturated water flow. They were used primarily in agricultural research for the purpose of optimizing moisture conditions to increase crop production. This focus has increasingly shifted to environmental research, with the primary concern now being the subsurface fate and transport of various agricultural and other contaminants. While the earlier models solved the governing equations (1) through (3) for relatively simplified system-independent boundary conditions (i.e. specified pressure heads or fluxes, and free drainage), models developed recently can cope with much more complex system-dependent boundary conditions evaluating surface flow and energy balances and accounting for the simultaneous movement of water, vapor, and heat. There are also composite models which simulate the processes both in unsaturated and saturated zones and other components of hydrological cycle. Few of the widely used unsaturated flow and composite models have been listed in Table 2.

Concluding Remarks

Forecasting water movement and contaminant transfer on a field-scale according to the existing monitoring and modelling methods is a challenging undertaking. The accuracy of the obtained predictions depend to a large extent upon the accuracy of available model input parameters and upon proper conceptualization of soil heterogeneity and other system complexities, such as the possible presence of non-equilibrium flow and transport, including preferential flow. Processes are often described and their parameters measured on a much smaller scale than those for which the model predictions are being sought. Consequently, many model parameters often need to be calibrated so that they reflect the bulk behaviour of the heterogeneous system, in which case they can be used for larger scale predictions.

Models can help guide field observations by identifying which parameters and processes control system behavior. Even so, you will find sizable unknowns in forecasts mostly because of our incapacity to demonstrate in depth spatial distributions of soil hydraulic attributes on the field-scale. Due to the high costs of data acquisition, few field measurements are usually available for the characterization of flow and contaminant transportation, although the spatial circulation of a contaminant plume could be very uneven. New measuring techniques that provide model parameters on the scale at which predictions are made are badly needed for successful applications of unsaturated flow and transport models in a predictive mode at the larger scale. Also, more research associated with the flow of water along with contaminant transport in the unsaturated zone of aquifers containing fractures and karstic conduits is needed for future investigations.

One may expect that unsaturated zone flow and transport models will be used increasingly for integrating fundamental knowledge about the vadose zone to yield tools for developing cost-effective, yet technically sound strategies for resource management and pollution remediation and prevention.

Abstract

Unsaturated zone transport models are indispensable tools for analyzing complex environmental pollution problems, and for developing practical management strategies. A quantitative study of water flow and contaminant transport in the unsaturated zone is a key factor in the improvement and protection of the quality of groundwater supplies. This is the region bounded above by the land surface and below by the groundwater table. It is the region through which water derived from precipitation and irrigation infiltrates and transports contaminants to reach the groundwater. This article presents an overview of the modelling process for water flow and contaminant transport in the unsaturated zone, input data requirements and related software packages.

Introduction

The unsaturated zone is the region through which water, together with pollutant carried by the water, must pass to reach the groundwater. Therefore various processes occurring within the unsaturated zone play a major role in determining both the quality and quantity of water recharging into the groundwater. A quantitative study of water flow and contaminant transport in the unsaturated zone is a key factor in the improvement and protection of the quality of groundwater supplies.

Numerous simulation models for water flow and solute transportation in the unsaturated region are now being used increasingly for numerous functions in both research and administration. Modelling techniques vary from straightforward analytical or even semi-analytical methods, to sophisticated numerical codes. Although analytical and even semi-analytical methods remain widely used for certain functions, the growing power of PCs along with the progression of more precise and numerically solid outcomes have inspired significantly broader usage of numerical codes in recent years. The extensive utilization of numerical models is additionally greatly improved by their accessibility in both the commercial and public areas, and by the advancement of innovative graphics-based interfaces which could significantly simplify their use.

Mathematical Equations of Water and Transport in Unsaturated Soils

Analytical, semi-analytical, and numerical models are usually based on the following three governing equations for water flow, solute transport, and heat movement, respectively:

Suitable simplifications (mostly for analytical approaches) or extensions thereof (e.g. for two- and three-dimensional systems) are also employed. In equation (1), frequently known as the Richards equation, h is the pressure head, z is the perpendicular coordinate positive upwards, θ is the water content, t is time, S is a sink term representing root water uptake or some other sources or sinks, and the hydraulic conductivity function (unsaturated) K(h) is, often given as the product of the relative hydraulic conductivity, Kr, and the saturated hydraulic conductivity, Ks. In equation (2), called the convection-dispersion equation (CDE), c is the concentration of the solution, D is the dispersion coefficient comprising both hydrodynamic dispersion and molecular diffusion, the retardation variable that represents adsorption (R), the volumetric fluid flux density (q), and Φ is a sink/source term that accounts for various zero- and first order or other reactions. In equation (3), T is temperature, λ is the apparent thermal conductivity, and C and Cw are the volumetric heat capacities of the soil and the liquid phase, respectively.

Solutions of the Richards equation (1) require knowledge of the unsaturated soil hydraulic functions, that is, the soil water retention curve, θ(h), describing the relationship between the water content θ and the pressure head h, and the hydraulic conductivity function (unsaturated),

K(h), defining the hydraulic conductivity K as a function of h or θ. While under certain conditions (i.e. for linear sorption, a concentration-independent sink term Φ, and a steady flow field) equations (2) and (3) are linear equations, equation (1) is often very nonlinear due to the nonlinearity of the soil's hydraulic attributes. Consequently, numerous analytical methods were derived previously for equations (2) and (3) and these kinds of analytical methods today are popular for evaluating solute and heat transport under steady-state conditions. Although a large number of analytical solutions of (1) exist, they can generally be applied only to drastically simplified problems. Most of the functions for water flow inside the vadose zone demands a numerical approach of the Richards formula.

Required Input Data

Simulation of water qualities in the unsaturated areas needs input info regarding the model parameters, the geometry of the system, the boundary conditions and, when simulating transient flow, initial conditions. With geometry parameters, the dimensions of the problem domain are defined. With the physical parameters, the physical properties of the system under consideration are described. With respect to the unsaturated zone, it concerns h(θ) (characteristics of soil water), and K(θ) (hydraulic conductivity).

For an appropriate explanation of the unsaturated flow, a correct description of the two hydraulic functions, K(θ) and h(θ), is important. K(θ), the hydraulic conductivity, declines significantly when the moisture content θ, declines from saturation. The experimental approach to calculating K(θ) at different moisture contents is fairly complicated but not too trustworthy. Additional methods were recommended to determine the K( θ ) function from more conveniently measurable characterizing qualities of the soil. In many studies, the hydraulic conductivity of the unsaturated soil is defined as product of a non-linear function of the effective water saturation, together with hydraulic conductivity at saturation. The relationship is shown by

… (4)

where,

The value of n is found to be 3.5 for coarse textured soils. n will vary with soil type. In literature, established empirical correlation between n and soil characteristic is available.

The relationship between the soil water pressure head h(θ) and moisture content θ, usually termed as the water retention curve or the soil moisture characteristic, is basically determined by the textural and the structural composition of the soil. Also the organic matter content may have an influence on the relationship. A characteristic feature of the water retention curve is that suction head (-h) decreases fairly rapidly with increasing moisture content. Hysteresis results might emerge, and, rather than being a single valued association, the h-θ connection includes a group of curves. The actual curve will have to be determined from the history of wetting and drying.

When root water uptake is additionally modelled, the criteria describing the association between root water uptake together with soil water state has to be supplied, along with plant details. If a functional flux-head association is employed as a lower limit stipulation, the criteria identifying the relationship between surface water and groundwater and, if required, the vertical resilience of inadequately porous layers must be provided.

The number and kind of criteria necessary for modeling flow as well as transport procedures in soils are dependent on the kind of model selected. These parameters can be categorized as control parameters (controlling the operation of the computer code), discretization data (grid and time stepping), and material parameters. The material criteria could be gathered in 6 units (Jury and Valentine, 1986) – static soil properties, water transport and retention attributes, time-dependent criteria, soil adsorption variables,standard chemical characteristics and tortuosity functions. Table 1 lists many of the relevant material model parameters.

Modelling of Unsaturated Flow

Analytical solutions, when available, provide better insight into the physics behind the transport phenomena and are computationally efficient and simple to use. However, analytical approaches are for the most part limited to situations of simple geometry domains, linear governing equations, and homogeneous systems. Along with efficient numerical methods and rapidly updated computer hardware, a large number of numerical models have been developed. However, the numerical technique cannot replace the analytical approaches completely, since numerical methods themselves need verification against analytical solutions because of discretization errors and convergence and stability problems that may be especially troublesome for advection-dominated and nonlinear adsorption processes.

Analytical approaches to the Richards formula for unsaturated flow under certain boundary and preliminary conditions are difficult to obtain because of the nonlinearity in soil hydraulic parameters. This difficulty is exaggerated in the case where the soil is heterogeneous. Usually, you have to depend on numerical methods for forecasting water motion in unsaturated soils, even for soils that are homogeneous. However, numerical approaches often suffer from convergence and mass balance problems.

Analytical solutions are only applicable to highly simplified systems, and are not well suited for situations normally encountered in the field. Initially, finite difference methods were mainly developed to predict unsaturated flow only; but finite element solution techniques also became available later. The nonlinearity of Richards equation is usually solved using an iterative procedure such as Newton or Picard methods. Perhaps the most important advantage of finite element techniques over standard finite difference methods is the ability to describe irregular system boundaries in simulations more accurately, as well as easily including nonhomogeneous medium properties. For one-dimensional simulations, finite difference methods are just as good as finite element schemes. However, several authors suggest that finite element methods lead to more stable and accurate solutions, thus permitting larger time steps and/or coarser grid systems, and hence leading to computationally more efficient numerical schemes.

To numerically solve coupled systems of equations, the solution process requires some manipulation at each time step so that the dependence of one equation on the solution of the other is dealt with accurately. One way to overcome this is to use a fully implicit approach to solve the equations simultaneously. Any nonlinearity of the generated system can be handled by Newton’s method. The implicit nature of this scheme allows for larger time steps in the simulation to find stable solutions as compared to the time steps for explicit schemes. Another method of using the totally implicit scheme is to apply the blended implicit-explicit approach. Yet somehow, the specific portion of the scheme implies that this algorithm is currently dependent upon a stability constraint which significantly limits the measurement of the time step and opens up numerical artifacts.

For heterogeneous soils that contain macropores, a different modelling approach is needed, as the presence of macropores in these soils may form a separate network for water flow. The common approach is to introduce two-region flow domain, one for macropores and the other for the soil matrix. In each flow domain, hydraulic conductivity is given independently. One more strategy would be to think of the heterogeneous soil as a group of stream tubes. It is assumed that there is no exchange of water between these tubes, and that within each tube, the hydraulic conductivity is defined, but varies between tubes.

Modelling of Solute Transport

Movement of dissolved solutes in soils is often defined by the advection-dispersion formula. Analytical answers are produced for various perimeter and preliminary situations. Although these solutions are obtained for narrowly specified circumstances, there are numerous applications like the confirmation of computer codes, forecast of solute motion for large times or distances where the use of numerical models become impractical, and the analysis of transport criteria from soil column research. The majority of analytical solutions pertain to semi-infinite and infinite media. Solutions are obtained by a variety of mathematical techniques including Green’s functions, separation of variables, characteristics method, Laplace transforms and Fourier transforms.

Forecast of solute migration under field situations requires the concurrent solution of the solute transport and unsaturated flow equations. First approximations involve or assume steady flow and constant water contents. Because of the natural complexity of unsaturated flow, strategies for forecasting solute transport have depended mostly on finite difference or finite element estimates of the governing formulas. Since advection-diffusion transport equations normally do not have closed form analytical outcomes, it is important to have accurate numerical approximations. Any time diffusion governs the physical procedure, regular finite difference and element methods function effectively in solving these equations. But yet, if advection is the principal procedure, these equations could display numerous numerical problems. In fact, any standard finite difference or element method will produce a solution, which exhibits non-physical oscillations.

One of the distinctive features of the porous media on the field scale is the spatial heterogeneity of transport properties. These features have a distinct effect on the spatial distribution of contaminant concentration, as has been observed in field experiments and demonstrated by simulation of contaminant transport in unsaturated, heterogeneous soil. Description of the mixing process due to spatial variability of the unsaturated hydraulic conductivity has been advanced with the development of numerical solutions, which assume spatially variable soil properties; stochastic models; and stochastic stream tube models, which decompose the field into a set of independent vertical soil columns.

Unsaturated Zone Modelling Software

Most of the early models developed for studying processes in the near-surface environment focused mainly on variably saturated water flow. They were used primarily in agricultural research for the purpose of optimizing moisture conditions to increase crop production. This focus has increasingly shifted to environmental research, with the primary concern now being the subsurface fate and transport of various agricultural and other contaminants. While the earlier models solved the governing equations (1) through (3) for relatively simplified system-independent boundary conditions (i.e. specified pressure heads or fluxes, and free drainage), models developed recently can cope with much more complex system-dependent boundary conditions evaluating surface flow and energy balances and accounting for the simultaneous movement of water, vapor, and heat. There are also composite models which simulate the processes both in unsaturated and saturated zones and other components of hydrological cycle. Few of the widely used unsaturated flow and composite models have been listed in Table 2.

Concluding Remarks

Forecasting water movement and contaminant transfer on a field-scale according to the existing monitoring and modelling methods is a challenging undertaking. The accuracy of the obtained predictions depend to a large extent upon the accuracy of available model input parameters and upon proper conceptualization of soil heterogeneity and other system complexities, such as the possible presence of non-equilibrium flow and transport, including preferential flow. Processes are often described and their parameters measured on a much smaller scale than those for which the model predictions are being sought. Consequently, many model parameters often need to be calibrated so that they reflect the bulk behaviour of the heterogeneous system, in which case they can be used for larger scale predictions.

Models can help guide field observations by identifying which parameters and processes control system behavior. Even so, you will find sizable unknowns in forecasts mostly because of our incapacity to demonstrate in depth spatial distributions of soil hydraulic attributes on the field-scale. Due to the high costs of data acquisition, few field measurements are usually available for the characterization of flow and contaminant transportation, although the spatial circulation of a contaminant plume could be very uneven. New measuring techniques that provide model parameters on the scale at which predictions are made are badly needed for successful applications of unsaturated flow and transport models in a predictive mode at the larger scale. Also, more research associated with the flow of water along with contaminant transport in the unsaturated zone of aquifers containing fractures and karstic conduits is needed for future investigations.

One may expect that unsaturated zone flow and transport models will be used increasingly for integrating fundamental knowledge about the vadose zone to yield tools for developing cost-effective, yet technically sound strategies for resource management and pollution remediation and prevention.