HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Google Scholar Articles by Simunek, J. Articles by van Genuchten, M. Th. PubMed Articles by Simunek, J. Articles by van Genuchten, M. Th. GeoRef GeoRef Citation Agricola Articles by Simunek, J. Articles by van Genuchten, M. Th. Related Collections Dual Porosity/Permeability Models Preferential Flow Models Solute Transport Models

SPECIAL SECTION: VADOSE ZONE MODELING

Modeling Nonequilibrium Flow and Transport Processes Using HYDRUS

Jirka imnek^a,* and Martinus Th. van Genuchten^b

^a Dep. of Environmental Sciences, Univ. of California, Riverside, CA 92521
^b USDA-ARS, U.S. Salinity Lab., 450 West Big Springs Rd., Riverside, CA 92507
^* Corresponding author (Jiri.Simunek{at}ucr.edu).

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.

Received 17 April 2007.

	ABSTRACT

TOP ABSTRACT INTRODUCTION Overview of Conceptual Models Specific Models for Water... Specific Models for Solute... Numerical Implementation Implications for Formulation of... Applications Discussion Conclusions Appendix REFERENCES

 
Accurate process-based modeling of nonequilibrium water flow and solute transport remains a major challenge in vadose zone hydrology. Our objective here was to describe a wide range of nonequilibrium flow and transport modeling approaches available within the latest version of the HYDRUS-1D software package. The formulations range from classical models simulating uniform flow and transport, to relatively traditional mobile-immobile water physical and two-site chemical nonequilibrium models, to more complex dual-permeability models that consider both physical and chemical nonequilibrium. The models are divided into three groups: (i) physical nonequilibrium transport models, (ii) chemical nonequilibrium transport models, and (iii) physical and chemical nonequilibrium transport models. Physical nonequilibrium models include the Mobile-Immobile Water Model, Dual-Porosity Model, Dual-Permeability Model, and Dual-Permeability Model with Immobile Water. Chemical nonequilibrium models include the One Kinetic Site Model, the Two-Site Model, and the Two Kinetic Sites Model. Finally, physical and chemical nonequilibrium transport models include the Dual-Porosity Model with One Kinetic Site and the Dual-Permeability Model with Two-Site Sorption. Example calculations using the different types of nonequilibrium models are presented. Implications for the formulation of the inverse problem are also discussed. The many different models that have been developed over the years for nonequilibrium flow and transport reflect the multitude of often simultaneous processes that can govern nonequilibrium and preferential flow at the field scale.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION Overview of Conceptual Models Specific Models for Water... Specific Models for Solute... Numerical Implementation Implications for Formulation of... Applications Discussion Conclusions Appendix REFERENCES

 
THERE IS increasing evidence that flow and transport processes in soils often cannot be described using classical models that assume uniform flow and transport (e.g., Nkedi-Kizza et al., 1984; Hendrickx and Flury, 2001; Pot et al., 2005; Köhne et al., 2006). Many laboratory and field experiments have demonstrated the presence of nonequilibrium flow and transport conditions in soils. Nonequilibrium water flow and solute transport in the unsaturated zone can be simulated at present by means of a large number of models of various degrees of complexity and dimensionality. Modeling approaches range from relatively simple analytical solutions for solute transport (e.g., van Genuchten, 1981; Toride et al., 1993) to complex numerical codes (e.g., imnek et al., 2005; Jacques and imnek, 2005). While such programs as STANMOD (imnek et al., 1999) that implement analytical solutions undoubtedly will remain useful for simplified analyses of solute transport during steady-state flow (e.g., for analyzing solute breakthrough curves measured in the laboratory, or for initial or approximate analysis of field-scale transport problems), numerical models are now increasingly being used since they can be applied more readily than analytical models to realistic laboratory and field problems. The use of numerical models has been further popularized during the last 20 yr or so because of the appearance of appropriate software packages in both the public and commercial domains and the development of increasingly sophisticated graphics-based interfaces that can simplify their use tremendously.

Attempts to describe nonequilibrium transport have traditionally been developed along two lines: physical and chemical nonequilibrium models (van Genuchten and Cleary, 1979). While physical nonequilibrium models assume that nonequilibrium flow or transport is caused by physical factors (e.g., van Genuchten and Wierenga, 1976), chemical nonequilibrium models assume that chemical factors are the cause of nonequilibrium transport (e.g., Selim et al., 1976; van Genuchten and Wagenet, 1989). Only a few researchers (e.g., Brusseau et al., 1989; Selim et al., 1999; Pot et al., 2005) have combined the physical and chemical nonequilibrium approaches to account for both possible causes of nonequilibrium, thus improving the description of solute transport in soils.

Over the years, several publicly available numerical codes have been developed that consider a number of options for simulating nonequilibrium water flow and solute transport (e.g., Pruess, 1991; Jarvis, 1994; van Dam et al., 1997). Unique to the HYDRUS-1D software package (imnek et al., 2005, 2008) is the wide range of approaches that can be selected for simulating nonequilibrium processes. Our objective here is to describe the large number of physical and chemical nonequilibrium approaches available in the latest version of HYDRUS-1D (imnek et al., 2008). The models range from classical models simulating uniform flow and transport, to traditional dual-porosity physical and two-site chemical nonequilibrium models, to complex dual-permeability models that consider both physical and chemical nonequilibrium. Since our focus is mainly on solute transport and since the nonequilibrium models for water flow have been reviewed relatively recently (imnek et al., 2003), we only very briefly review here the governing water flow equations to define variables later used in the solute transport equations. We will present several examples calculated with the different nonequilibrium approaches to show the effect of various transport and reaction parameters, and to demonstrate the consequences of increased complexity in the models. Implications for the formulation of the inverse problem are also discussed.

	Overview of Conceptual Models

TOP ABSTRACT INTRODUCTION Overview of Conceptual Models Specific Models for Water... Specific Models for Solute... Numerical Implementation Implications for Formulation of... Applications Discussion Conclusions Appendix REFERENCES

 
A large number of alternative physical and chemical nonequilibrium models can be formulated. Figures 1 and 2 show schematics of a range of possible physical equilibrium and nonequilibrium models for water flow and solute transport. Figure 3 shows similar schematics of various chemical nonequilibrium models that have been incorporated into HYDRUS-1D. Traditional flow and transport models are based on the classical description of uniform flow and transport in soils (the Uniform Flow Model in Fig. 1a and 2a). In this model, the porous medium is viewed as a collection of impermeable soil particles (or of impermeable soil aggregates or rock fragments), separated by pores or fractures through which flow and transport takes place. Variably saturated water flow through such a porous system is usually described using the Richards equation and solute transport using the classical advection–dispersion equation. Definitions of various water contents and concentrations used in the different models are given in the Appendix.

View larger version (17K): [in this window] [in a new window]   FIG. 1. Conceptual physical nonequilibrium models for water flow and solute transport. In the plots, is the water content, mo and im in (b) and (c) are water contents of the mobile and immobile flow regions, respectively; M and F in (d) are water contents of the matrix and macropore (fracture) regions, respectively, and M,mo, M,im, and F in (e) are water contents of the mobile and immobile flow regions of the matrix domain, and in the macropore (fracture) domain, respectively; c are concentrations of corresponding regions, with subscripts having the same meaning as for water contents, while S is the total solute content of the liquid phase.

View larger version (28K): [in this window] [in a new window]   FIG. 2. Conceptual physical nonequilibrium models for water flow and solute transport.

The mobile–immobile water model can be further expanded by assuming that both water and solute can move into and out of the immobile domain (imnek et al., 2003), leading to the Dual-Porosity Model in Fig. 1c and 2c. While the water content inside of the soil particles or aggregates is assumed to be constant in the Mobile–Immobile Water Model, it can vary in the Dual-Porosity Model since the immobile domain is now allowed to dry out or rewet during drying and wetting processes. Water flow into and out of the immobile zone is usually described using a first-order rate process. Solute can move into the immobile domain of the Dual-Porosity Model by both molecular diffusion and advection with flowing (exchanging) water. Since water can move from the main pore system into the soil aggregates and vice versa, but not directly between the aggregates themselves, water in the aggregates can be considered immobile from a larger scale point of view.

The limitation of water not being allowed to move directly between aggregates is overcome in a Dual-Permeability Model (e.g., Gerke and van Genuchten, 1993a,b). Water and solutes in such models also move directly between soil aggregates, as shown in Fig. 1d and 2d. Dual-permeability models assume that the porous medium consists of two overlapping pore domains, with water flowing relatively fast in one domain (often called the macropore, fracture, or interporosity domain) when close to full saturation, and slow in the other domain (often referred to as the micropore, matrix, or intraporosity domain). Like the Dual-Porosity Model, the Dual-Permeability Model allows the transfer of both water and solutes between the two pore regions.

Finally, the Dual-Permeability Model can be further refined by assuming that inside of the matrix domain an additional immobile region exists into which solute can move by molecular diffusion (the Dual-Permeability Model with MIM in Fig. 1e and 2e).

Chemical nonequilibrium models implemented into HYDRUS-1D are schematically shown in Fig. 3 . The simplest chemical nonequilibrium model assumes that sorption is a kinetic process (the One Kinetic Site Model in Fig. 3a), usually described by means of a first-order rate equation. This model can be expanded into a Two-Site Sorption model by assuming that the sorption sites can be divided into two fractions (e.g., Selim et al., 1976; van Genuchten and Wagenet, 1989). The simplest two-site sorption model arises when sorption on one fraction of the sorption sites is assumed to be instantaneous, while kinetic sorption occurs on the second fraction (Two-Site Model in Fig. 3b). This model can be further expanded by assuming that sorption on both fractions is kinetic and proceeds at different rates (the Two Kinetic Sites Model in Fig. 3c). The Two Kinetic Sites Model reduces to the Two-Site Model when one rate is so high that it can be considered instantaneous, to the One Kinetic Site Model when both rates are the same, or to the chemical equilibrium model when both rates are so high that they can be considered instantaneous.

View larger version (17K): [in this window] [in a new window]   FIG. 3. Conceptual chemical nonequilibrium models for reactive solute transport. In the plots, is the water content, mo and im in (d) are water contents of the mobile and immobile flow regions, respectively; M and F in (e) are water contents of the matrix and macropore (fracture) regions, respectively; c are concentrations of the corresponding regions, se are sorbed concentrations in equilibrium with the liquid concentrations of the corresponding regions, and sk are kinetically sorbed solute concentrations of the corresponding regions.

The resulting combined physical and chemical nonequilibrium approach may be simulated with HYDRUS-1D using the Dual-Porosity Model with One Kinetic Site (Fig. 3d). This model considers water flow and solute transport in a dual-porosity system (or a medium with mobile–immobile water) while assuming that sorption in the immobile zone is instantaneous. Following the two-site kinetic sorption concept, however, the sorption sites in contact with the mobile zone are now divided into two fractions, subject to either instantaneous or kinetic sorption. Since the residence time of solutes in the immobile domain is relatively large, equilibrium probably exists between the solution and the sorption complex here, in which case there is no need to consider kinetic sorption in the immobile domain. The model, on the other hand, assumes the presence of kinetic sorption sites in contact with the mobile zone since water can move relatively fast in the macropore domain and thus prevent chemical equilibrium.

Finally, chemical nonequilibrium can also be combined with the Dual-Permeability Model. This last nonequilibrium option implemented into HYDRUS-1D (the Dual-Permeability Model with Two-Site Sorption in Fig. 3e) assumes that equilibrium and kinetic sites exist in both the macropore (fracture) and micropore (matrix) domains.

A complete list of the different models summarized here, and described in more detail below, is given in Table 1 , including the specific equations used for the water flow and solute transport models. The table additionally lists the various parameters, and their total number, that will appear in any particular solute transport model when used for steady-state flow conditions.

	Specific Models for Water Flow

TOP ABSTRACT INTRODUCTION Overview of Conceptual Models Specific Models for Water... Specific Models for Solute... Numerical Implementation Implications for Formulation of... Applications Discussion Conclusions Appendix REFERENCES

Uniform Flow Model
Numerical models for water flow in soils (Fig. 1a and 2a) are usually based on the following equation:

[1]

or its extensions (e.g., for two- and three-dimensional systems). In Eq. [1], often referred to as the Richards equation, z is the vertical coordinate positive upward [L], t is time [T], h is the pressure head [L], is the water content [L³ L^–3], S is a sink term representing root water uptake or some other source or sink [T^–1], and K(h) is the unsaturated hydraulic conductivity function, often given as the product of the relative hydraulic conductivity, K_r (dimensionless), and the saturated hydraulic conductivity, K_s [L T^–1]. Solutions of the Richards Eq. [1] require knowledge of the unsaturated soil hydraulic functions made up of the soil water retention curve, (h), which describes the relationship between the water content and the pressure head h, and the unsaturated hydraulic conductivity function, K(h), which defines the hydraulic conductivity K as a function of h or . HYDRUS-1D considers both relatively traditional models (Brooks and Corey, 1964; van Genuchten, 1980) for the hydraulic functions, as well as more recent alternative single- (e.g., Kosugi, 1996) and dual-porosity (Durner, 1994) models.

Dual-Porosity Model
Dual-porosity models (Fig. 1c and 2c) assume that water flow is restricted to the macropores (or interaggregate pores and fractures), and that water in the matrix (intraaggregate pores or the rock matrix) does not move at all. This conceptualization leads to two-region type flow and transport models (van Genuchten and Wierenga, 1976) that partition the liquid phase into mobile (flowing, interaggregate), _mo, and immobile (stagnant, intraaggregate), _im, regions [L³ L₃^–]:

[2]

The dual-porosity formulation for water flow can be based on a mixed formulation of the Richards Eq. [1] to describe water flow in the macropores (the preferential flow pathways) and a mass balance equation to describe moisture dynamics in the matrix as follows (imnek et al., 2003):

[3]

where S_mo and S_im are sink terms for the mobile and immobile regions, respectively [T^–1], and _w is the transfer rate for water between the inter- and intraaggregate pore domains [T^–1]. imnek et al. (2003) and Köhne et al. (2004) discussed different formulations that can be used to evaluate the mass transfer rate _w.

Dual-Permeability Model
Different dual-permeability approaches (Fig. 1d, 1e, 2d, and 2e) may be used to describe flow and transport in structured media. While some models invoke similar equations for flow in the fracture and matrix regions, others use different formulations for the two regions. A typical example of the first approach, implemented in HYDRUS-1D, is the work of Gerke and van Genuchten (1993a,b, 1996), who applied the Richards equation to each of the two pore regions. The flow equations for the macropore (fracture) (subscript f) and matrix (subscript m) pore systems in their approach are given by

[4]

where w is the ratio of the volumes of the macropore or fracture domain and the total soil system (dimensionless). This approach is relatively complicated in that the model requires characterization of water retention and hydraulic conductivity functions (potentially of different form) for both pore regions, as well as a hydraulic conductivity function of the fracture–matrix interface. Note that the water contents, _f and _m in Eq. [4], have different meanings than in Eq. [3], where they represented water contents of the total pore space (i.e., = _mo + _im), while here they refer to water contents of the two separate (fracture or matrix) pore domains such that

[5]

Hence, lowercase subscripts in the dual-permeability model refer to the local (pore-region) scale, while uppercase subscripts refer to the global (total soil medium) scale.

	Specific Models for Solute Transport

TOP ABSTRACT INTRODUCTION Overview of Conceptual Models Specific Models for Water... Specific Models for Solute... Numerical Implementation Implications for Formulation of... Applications Discussion Conclusions Appendix REFERENCES

 
Uniform Transport
Solute transport in numerical models is usually described using the relatively standard advection–dispersion equation (Fig. 1a and 2a) of the form

[6]

or various extensions thereof (e.g., for two- and three-dimensional systems, or for multiple phases or components). In Eq. [6], c is the solution concentration [M L^–3], s is the sorbed concentration [M M^–1], D is the dispersion coefficient accounting for both molecular diffusion and hydrodynamic dispersion [L² T^–1], q is the volumetric fluid flux density [L T^–1] evaluated using the Darcy–Buckingham law, and is a sink–source term that accounts for various zero- and first-order or other reactions [M L^–3 T^–1].

While HYDRUS-1D considers a general nonlinear sorption equation that can be simplified into a Langmuir or Freundlich isotherm (imnek et al., 2005), for simplicity we assume here only linear adsorption of the form

[7]

where K_d is the distribution coefficient [L³ M^–1]. Linear sorption leads to the following definition of the retardation factor R (dimensionless):

[8]

In our examples below we will be using Eq. [7] for linear sorption and applying the resulting definition of R to all of the liquid domains involved (e.g., the total liquid phase, the mobile and immobile regions, or the matrix and macropore domains) with their appropriate parameters. Although considered in HYDRUS-1D, the effect of molecular diffusion will be neglected in the various examples. The dispersion coefficient D accounting only for hydrodynamic dispersion [L² T^–1] is thus defined as

[9]

where is the dispersivity [L] and v the average pore velocity [L T^–1]. The same definition of the dispersion coefficient will be used for all mobile phases. HYDRUS-1D additionally considers molecular diffusion for transport in the gaseous phase, which will not be discussed here either. Finally, in the examples to follow we will use the time T needed to reach one pore volume of effluent. This time is given by

[10]

where L is some distance (e.g., column length) within the transport domain being considered. For the physical nonequilibrium models discussed below, the pore volume for different domains is defined using corresponding water contents and fluxes.

Physical Nonequilibrium Transport Models
Mobile–Immobile Water and Dual-Porosity Models
The concept of two-region, dual-porosity type solute transport (Fig. 1b, 1c, 2b, and 2c) was implemented already in earlier versions (1.0 and 2.0) of HYDRUS-1D to permit consideration of physical nonequilibrium transport. While the physical nonequilibrium transport model in the earlier versions was combined only with uniform water flow Eq. [1], Version 3.0 of HYDRUS-1D was expanded to also consider the dual-porosity water flow model Eq. [3] with a transient immobile water content. In both implementations, the governing solute transport equations are as follows:

[11a]

[11b]

[11c]

in which solute exchange between the two liquid regions is modeled as the sum of an apparent first-order diffusion process and advective transport (where applicable). In Eq. [11], c_mo and c_im are concentrations of the mobile and immobile regions [M L^–3], respectively; s_mo and s_im are sorbed concentrations of the mobile and immobile regions [M M^–1], respectively; D_mo is the dispersion coefficient in the mobile region [L² T^–1], q_mo is the volumetric fluid flux density in the mobile region [L T^–1], _mo and _im are sink–source terms that account for various zero- and first-order or other reactions in both regions [M L^–3 T^–1]; f_mo is the fraction of sorption sites in contact with the mobile water content (dimensionless), _mim is the mass transfer coefficient [T^–1], and _s is the mass transfer term for solutes between the mobile and immobile regions [M L^–3 T^–1]. Equation [11a] describes solute transport in the mobile (macropore) zone, Eq. [11b] is a mass balance for the immobile (micropore) domain, while Eq. [11c] (_s) describes the rate of mass transfer between the mobile and immobile domains. The second (advective) term of _s in Eq. [11] is equal to zero for the Mobile–Immobile Model since that model does not consider water flow between the two regions. In the Dual-Porosity Model, c* is equal to c_mo for _w > 0 and c_im for _w < 0.

Dual-Permeability Model
Analogous to Eq. [4], the dual-permeability formulation for solute transport is based on advection–dispersion type equations for transport in both the fracture and matrix regions as follows (Gerke and van Genuchten, 1993a,b) (Fig. 1d and 2d):

[12a]

[12b]

[12c]

The variables in Eq. [12] have similar meanings as in Eq. [11], except that they refer now to two overlapping domains, i.e., the matrix (subscript m) and fracture (subscript f) domains. Equation [12a] describes solute transport in the fracture domain, Eq. [12b] transport in the matrix domain, and Eq. [12c] advective–dispersive mass transfer between the fracture and matrix domains. Equation [12] assumes complete advective–dispersive transport descriptions for both the fractures and the matrix. Van Genuchten and Dalton (1986) and Gerke and van Genuchten (1996), among others, discussed possible expressions for the first-order solute mass transfer coefficient, _dp [T^–1].

Dual-Permeability Model with Immobile Water
The Dual-Permeability Model with Immobile Water (Fig. 1e and 2e) assumes that the liquid phase of the matrix can be further partitioned into mobile (flowing), _m,m [L³ L^–3], and immobile (stagnant), _im,m [L³ L^–3], regions as follows:

[13]

where _m is the volumetric water content of the matrix pore system [L³ L^–3]. The governing advection–dispersion equation for transport in the matrix region (Eq. [12b]) is then replaced with the modified equations (Eq. [11]) (e.g., Pot et al., 2005) to yield

[14a]

[14b]

[14c]

[14d]

[14e]

where c_im,m and c_m,m are solute concentrations in the immobile and mobile zones of the matrix region [M L^–3], respectively; _m,m and _im,m represent various reactions in the mobile and immobile parts of the matrix [M L^–3 T^–1], respectively; f_m is again the fraction of sorption sites in contact with the mobile region of the matrix (dimensionless), _dpm is the mass transfer coefficient between mobile and immobile zones of the matrix region [T^–1], and _s* is the mass transfer term for solutes between the mobile and immobile regions of the matrix domain [M L^–3 T^–1]. Equation [14a] now describes solute transport in the fracture domain, Eq. [14b] transport in the mobile zone of the matrix domain, Eq. [14c] is a mass balance for the immobile zone of the matrix domain, Eq. [14d] describes mass transfer between the fracture and matrix domains, while Eq. [14e] describes mass transfer between the mobile and immobile zones within the matrix domain.

Chemical Nonequilibrium Transport Models
One Kinetic Site Model
When sorption in the Uniform Transport Model is considered (Fig. 3a) to be kinetic, Eq. [6] needs to be supplemented with an equation describing the kinetics of the sorption process. This is usually done by assuming a first-order process as follows:

[15]

where s_e^k is the sorbed concentration that would be reached at equilibrium with the liquid-phase concentration [M M^–1], s^k is the sorbed concentration of the kinetic sorption sites [M M^–1], _k is a first-order rate constant describing the kinetics of the sorption process [T^–1], and _k represents a sink–source term that accounts for various zero- and first-order or other reactions at the kinetic sorption sites [M L^–3 T^–1].

Two-Site Model
Similarly to the mobile–immobile water concept, the concept of two-site sorption (Selim et al., 1976; van Genuchten and Wagenet, 1989) (Fig. 3b) was implemented already in Versions 1.0 and 2.0 of HYDRUS-1D to permit consideration of nonequilibrium adsorption–desorption reactions. The two-site sorption concept assumes that the sorption sites can be divided into two fractions:

[16]

Sorption s^e [M M^–1], on one fraction of the sites (Type 1 sites) is assumed to be instantaneous, while sorption s^k [M M^–1], on the remaining (Type 2) sites is considered to be a first-order kinetic rate process. The system of equations describing the Two-Site Model is given by:

[17a]

[17b]

[17c]

[17d]

where f_e is the fraction of exchange sites assumed to be in equilibrium with the liquid phase (dimensionless), and _k is a first-order rate constant [T^–1]. Equation [17a] describes solute transport in the total system, Eq. [17b] equilibrium sorption onto the instantaneous sorption sites, Eq. [17c] is a mass balance of the kinetic sorption sites (van Genuchten and Wagenet, 1989), while Eq. [17d] represents the sorbed concentration of the kinetic sites when equilibrium would be reached with the liquid-phase concentration.

Two Kinetic Sites Model
To facilitate simulations of the transport of colloids or microorganisms (such as viruses and bacteria), Version 3.0 of HYDRUS-1D also implemented a Two Kinetic Sites Model (Fig. 3c) using the attachment–detachment approach:

[18]

where s₁^k and s₂^k are sorbed concentrations of the first and second fractions of kinetic sorption sites [M M^–1], respectively; k_a1 and k_a2 are attachment coefficients for the first and second fractions of kinetic sorption sites [T^–1], respectively; k_d1 and k_d2 are detachment coefficients for the first and second fractions of kinetic sorption sites [T^–1], respectively; and _k1 and _k2 represent sink–source terms for the first and second fractions of kinetic sorption sites [M L^–3 T^–1], respectively. Note that the Two Kinetic Sites Model can be used (and often is used) to describe different processes. While the first kinetic process could be chemical attachment, the second kinetic process could represent physical straining (e.g., Bradford et al., 2004; Gargiulo et al., 2007, 2008). Note that in Eq. [18] we do not give the nonlinear blocking coefficients accounting for, for example, Langmuirian blocking to attachment sites or depth-dependent straining that are considered in HYDRUS-1D (e.g., Bradford et al., 2004).

It is easily shown that the formulation based on attachment–detachment coefficients is mathematically identical to the formulation using first-order mass transfer coefficients. For example, by comparing Eq. [15] with [18] we have

[19]

Physical and Chemical Nonequilibrium Transport Models
Dual-Porosity Model with One Kinetic Site
This model (Fig. 3d) is similar to the Dual-Porosity Model (Eq. [2]) in that the porous medium is divided into mobile and immobile domains such that = _mo + _im. The current model, however, additionally divides the sorption sites in contact with the mobile zone, similarly to the Two-Site Model (Eq. [16]), into two fractions involving instantaneous and kinetic sorption such that the total sorption concentration at equilibrium is given by

[20]

where s_mo^e is the sorbed concentration in equilibrium with the liquid-phase concentration of the mobile region of the Dual-Porosity Model [M M^–1], s_mo,e^k is the sorbed concentration of the kinetic sites in contact with the mobile region of the Dual-Porosity Model when at equilibrium [M M^–1], f_mo is the fraction of sorption sites in contact with mobile water (the remainder is in contact with immobile water), and f_em is the fraction of sorption sites in equilibrium with the mobile liquid phase (the remaining kinetic sites are also in contact with the mobile liquid phase). The complete Dual-Porosity Model with One Kinetic Site is described using the following equations:

[21a]

[21b]

[21c]

[21d]

[21e]

[21f]

[21g]

where _ph and _ch are first-order rate constants [T^–1] accounting for physical and chemical rate processes, respectively; _s1 is the mass transfer term for solute exchange between the mobile and immobile regions [M L^–3 T^–1]; _s2 represents mass transfer to the kinetic sorption sites in the mobile region [M L^–3 T^–1]; and _mo, _im, and _mo,k represent sink–source terms for the equilibrium phases in the mobile zone, the immobile zone, and the kinetic sorption sites [M L^–3 T^–1], respectively. Equation [21a] describes transport in the mobile phase, Eq. [21b] is a mass balance for the immobile phase, and Eq. [21c] a mass balance for the kinetic sorption sites in contact with the mobile zone. Equations [21d] and [21e] describe mass transfer rates between the mobile and immobile zones and to the kinetic sorption sites, respectively, while Eq. [21f] and [21g] represent sorption onto the equilibrium and kinetic sorption sites in contact with the mobile zone, respectively.

Dual-Permeability Model with Two-Site Sorption
Finally, simultaneous physical and chemical nonequilibrium processes are implemented in HYDRUS-1D by assuming applicability of the Dual-Permeability Model (Gerke and van Genuchten, 1993a; imnek et al., 2003) and dividing the sorption sites of both the fracture and matrix domains into equilibrium and kinetic sites (Fig. 3e). This model leads to the following set of equations (Pot et al., 2005):

[22]

where s_m^k and s_f^k are sorbed concentrations of Type 2 (kinetic) sites in the matrix and fracture domains [M M^–1], respectively; f_m and f_f are fractions of the exchange sites assumed to be in equilibrium with the solution phases (dimensionless) of the matrix and fracture domains, respectively; _f, _m, _f,k, and _m,k represent reactions in the equilibrium phases of the fracture and matrix domains and at the kinetic sites of the fracture and matrix domains [M L^–3 T^–1], respectively; and _ch,m and _ch,f are again first-order rate constants for the matrix and fracture domains [T^–1], respectively. Note that the distribution coefficients can be different in the different regions (i.e., K_df K_dm).

	Numerical Implementation

TOP ABSTRACT INTRODUCTION Overview of Conceptual Models Specific Models for Water... Specific Models for Solute... Numerical Implementation Implications for Formulation of... Applications Discussion Conclusions Appendix REFERENCES

 
All of the models presented here were implemented in the HYDRUS-1D software package and as such are all solved with very similar numerical techniques. The Galerkin-type linear finite element method was used for spatial discretization of the governing partial differential equations, while finite difference methods were used to approximate temporal derivatives. A fully implicit finite difference scheme with Picard linearization was used to solve the Richards equation, while a Crank–Nicholson finite difference scheme was used for solution of the advection–dispersion equations. For the dual-permeability models, we always first solved the equations describing processes in the matrix, after which the equations describing processes in the fractures were solved. Complete details about the invoked numerical techniques are provided in the HYDRUS-1D technical manual (imnek et al., 2008).

	Implications for Formulation of the Inverse Problem

TOP ABSTRACT INTRODUCTION Overview of Conceptual Models Specific Models for Water... Specific Models for Solute... Numerical Implementation Implications for Formulation of... Applications Discussion Conclusions Appendix REFERENCES

 
The physical and chemical nonequilibrium models both involve relatively large numbers of parameters, many of which cannot be (or cannot easily be) measured independently. Some parameters need to be obtained by calibrating a particular model against laboratory or field measurements. This is usually done using a parameter estimation procedure (e.g., imnek and Hopmans, 2002; imnek et al., 2002) in which the sum of squared deviations between measurements and model predictions, organized in an objective function, is minimized. While the definition of the objective function for models that consider only uniform flow and transport may be relatively straightforward, the objective functions for calibration of the nonequilibrium models can become extremely complicated. For the equilibrium flow and transport models, one usually defines the objective function in terms of measured pressure heads, water contents, solution concentrations, or actual or cumulative water or solute fluxes. Note from Fig. 1a that the system is always described using unique values of the water content, pressure head, or solution concentration. In addition to the resident concentration, the objective function can also be defined using the flux concentration or the total solute mass at a specified location. The total solute mass must include not only the mass in the liquid phase but also the mass sorbed to either instantaneous or kinetic sorption sites (the latter if a kinetic sorption model is used).

By comparison, when physical nonequilibrium models are used, one can encounter different water contents such as the mobile water content, the immobile water content, the water content of the fracture domain, the water content of the matrix, or the total water content. The nonequilibrium models typically will also involve different types of concentrations, including concentrations of the mobile and immobile zones, sorbed concentrations associated with instantaneous and kinetic sorption sites, and total concentrations. Since different measurement techniques may lead to different types of concentrations (e.g., resident, flux, or total concentrations), a flexible inverse approach should also allow for the different modes of concentrations. Different types of water and solute fluxes can furthermore be defined in the nonequilibrium models, including total fluxes, fluxes in the mobile zone, or fluxes in the fracture and matrix domains. To provide users with flexibility during the model calibration process, these various water contents, concentrations, and water and solute fluxes must be considered in the objective function. Table 2 provides a list of variables in the different equilibrium and nonequilibrium models that can be incorporated in the objective function when HYDRUS-1D is used for parameter estimation purposes.

	Applications

TOP ABSTRACT INTRODUCTION Overview of Conceptual Models Specific Models for Water... Specific Models for Solute... Numerical Implementation Implications for Formulation of... Applications Discussion Conclusions Appendix REFERENCES

Uniform Transport
The uniform (equilibrium) solute transport model during steady-state water flow is fully defined in terms of four parameters (Table 1). Of these, two parameters (the water content, , and the fluid flux density, q) are related to water flow, and only two parameters (the dispersivity, , and the distribution coefficient, K_d) are directly related to solute transport. Figure 4 shows the well-known effects of the dispersivity (left) and the distribution coefficient (right) on calculated breakthrough curves. Notice that the red lines in Fig. 4 (left and right) represent the same simulation. While higher values of the dispersivity lead to earlier arrival of the solute and more pronounced tailing compared with lower values, larger values of the distribution coefficient delay the arrival of solute in the effluent. Note that the selected distribution coefficients of 0, 1, and 3 lead to retardation factors of 1, 4, and 10, respectively. Also note that one pore volume (Eq. [10]) for the selected parameters corresponds to 1 d.

View larger version (9K): [in this window] [in a new window]   FIG. 4. Breakthrough curves calculated using the Uniform Transport Model for a 10-cm-long soil column and the following parameters: q = 5 cm d–1, = 0.5, and b = 1.5 g cm–3; on the left, Kd = 1 cm3 g–1 and = 0.1, 1, and 10 cm; and on the right, = 1 cm and Kd = 0, 1, and 3 cm3 g–1.

View larger version (10K): [in this window] [in a new window]   FIG. 5. Breakthrough curves calculated using the Mobile–Immobile Water Model for a 10-cm-long soil column and the following parameters: q = 3 cm d–1, = 0.5, mo = 1 cm, Kd = 1 cm3 g–1, b = 1.5 g cm–3. On the left: fmo = 0.6, mo = 0.3, im = 0.2, and mim = 0.1, 0.5, 10 d–1. On the right: mim = 0.5 d–1 and (a) fmo = 0.4, mo = 0.2, im = 0.3; (b) fmo = 0.6, mo = 0.3, im = 0.2; and (c) fmo = 0.8, mo = 0.4, im = 0.1.

[23]

As shown in Fig. 6, the average outflow concentration quickly increases initially, similarly as for the fracture domain, but then stabilizes at about 0.5 since only about half of the outflow comes from the fracture domain. When solute arrives also in the outflow from the matrix domain, the average outflow concentration gradually increases again until it reaches unit concentration. The solute breakthrough curves start deviating from those calculated using the Uniform Flow Model individually for the matrix and macropore regions when the mass transfer coefficient is increased (red lines in Fig. 6) because of exchange of solute between the two regions. For relatively large values of the mass transfer coefficient, this causes rapid equilibration of the concentrations of the two domains. The breakthrough curves from the two regions will then start converging and will again resemble the results of the uniform transport model applied to the entire soil (blue lines in Fig. 6). For a fully equilibrated system, R equals 4 and one pore volume equals 0.88 d, causing the solute to arrive at approximately 3.5 d.

View larger version (21K): [in this window] [in a new window]   FIG. 6. Breakthrough curves calculated using the Dual-Permeability Model for a 10-cm-long soil column and the following parameters: qm = 3 cm d–1, qf = 30 cm d–1, = m = f = 0.5, w = 0.1, m = f = 1 cm, Kdm = Kdf = 1 cm3 g–1, b = 1.5 g cm–3, and dp = 0, 0.1, and 0.5 d–1. Matrix, fracture, and total breakthrough curves are represented by thin, medium, and thick lines, respectively.

View larger version (25K): [in this window] [in a new window]   FIG. 7. Breakthrough curves calculated using the Dual-Permeability Model with MIM for a 10-cm-long soil column and the following parameters: qm = 3 cm d–1, qf = 30 cm d–1, = m = f = 0.5, w = 0.1, m = f = 1 cm, Kdm = Kdf = 1 cm3 g–1, b = 1.5 g cm–3, dp = 0.1 d–1, dpm = 0.1 d–1, im,m (thim) = 0.0, 0.1, 0.3, and 0.4, and fm = 1, 0.8, 0.4, and 0.2, respectively. Matrix, fracture, and total breakthrough curves are represented by thin, medium, and thick lines, respectively.

View larger version (13K): [in this window] [in a new window]   FIG. 8. Breakthrough curves calculated using the One Kinetic Site Model for a 10-cm-long soil column and the following parameters: q = 5 cm d–1, = 0.5, = 1 cm, Kd = 1 cm3 g–1, b = 1.5 g cm–3, and k = 0.1, 0.5, and 1