The DualPerm Modulefor HYDRUS (2D/3D)
Simulating Two-Dimensional
Water Movement and Solute Transportin Dual-Permeability Variably-Saturated Porous Media
Version 1

Jirka Šimůnek1, Miroslav Šejna2, and M. Th. van Genuchten3
January 2012
1

Department of Environmental Sciences
University of California Riverside
Riverside, CA, 92521, USA
2

PC Progress, Ltd.
Prague, Czech Republic
3

Department of Mechanical Engineering
Federal University of Rio de Janeiro
Rio de Janeiro, Brazil© 2012 J. Šimůnek and M. Šejna. All rights reserved.

2

Table of Contents
Table of Contents................................................................................................................................. 3
List of Figures ...................................................................................................................................... 5
List of Tables ........................................................................................................................................ 7
Abstract ................................................................................................................................................ 9
1. Introduction ............................................................................................................................11
2. Dual-Permeability Water Flow ................................................................................................ 13
2.1.

Governing Flow Equations ..........................................................................................13

2.2.

Mass Water Transfer....................................................................................................14

2.3.

Soil Hydraulic Properties ............................................................................................15

3. Dual-Permeability Water Flow .............................................................................................17
3.1.

Governing Solute Transport Equations .......................................................................17

3.2.

Solute Mass Transfer ...................................................................................................17

4. Preprocessing and Input Data ..............................................................................................19
4.1.

Soil Hydraulic Model ...................................................................................................19

4.2.

Water Flow Parameters ................................................................................................19

4.3.

Solute Transport Parameters .......................................................................................21

4.4.

Solute Reaction Parameters .........................................................................................22

4.5.

Initial Conditions .........................................................................................................24

4.6.

Anisotropy ....................................................................................................................24

5. Post-Processing.......................................................................................................................25
5.1.

Results – Graphical Display ........................................................................................25

5.2.

The v_meanF.out Output File ......................................................................................26

6. Verification and Test Examples ............................................................................................27
6.1.

Example 1 (Test1d): Surface Flux Into Fracture. ..........................................................28

6.2.

Example 2 (Test1a3): Surface Flux Into Matrix. ...........................................................29

6.3.

Examples 3 (Test3) - Tension Disc Infiltration Experiment ...........................................30

References .....................................................................................................................................313

4

List of Figures
Figure 1.

Conceptual physical nonequilibrium models for water flow and solute transport. .11

Figure 2.

The Soil Hydraulic Model dialog window. .............................................................19

Figure 3.

The Water Flow Parameters dialog window. ..........................................................20

Figure 4.

The Solute Transport Parameters dialog window. ..................................................21

Figure 5.

Parts of the "Reaction Parameters for Solute - 1" dialog window. .........................23

Figure 6.

The "Domain Properties" part of the Data Tab of the Navigator Bar. ....................24

Figure 7.

The "Results - Graphical Display" part of the Data Tab of the Navigator Bar. ......25

Figure 8.

The "Results" part of the View Tab of the Navigator Bar. .....................................25

Figure 9.

The "Dual-Permeability Options" switch at the Edit (Tools) Bar. ..........................26

Figure 10.

Pressure heads in the matrix (a) and fracture (b) domains, water contents in thematrix (c) and fracture (d) domains, water mass exchange rates (e), solute mass
exchange rates (f), and concentrations in the matrix (h) and fracture (g) domains
calculated with the dual-permeability model of Gerke and van Genuchten[1993a,b]..................................................................................................................28

Figure 11.

Infiltration ad mass exchange fluxes (a), water contents in the matrix (b) andfracture (c) domains, water mass exchange rates (d) and concentrations in thematrix (e) and fracture (f) domains calculated with the dual-permeability model ofGerke and van Genuchten [1993a,b].......................................................................29

Figure 12.

Pressure head profiles for the matrix, isotropic fracture, and fracture with KxA/KzA=10,
and fracture with KxA/KzA=0.1 after 19 minutes.. .......................................................30

5

6

List of Tables
Table 1.

Water flow parameters for the dual-permeability module ...........................................20

Table 2.

Soil specific solute transport parameters... ..................................................................22

Table 3.

Solute specific solute transport parameters..................................................................22

Table 4.

Solute reaction parameters... ........................................................................................22

Table 5.

Duplicated outputs from the DualPerm module for the matrix and fracture domains...25

Table 6.

The v_meanF.out output file... .....................................................................................26

Table 7.

List of various test examples used to verify various new options of the DulPerm
module..........................................................................................................................277

8

AbstractŠimůnek, J., M. Šejna, and M. Th. van Genuchten, The DualPerm Module for HYDRUS (2D/3D)
Simulating Two-Dimensional Water Movement and Solute Transport in Dual-Permeability PorousMedia, Version 1.0, PC Progress, Prague, Czech Republic, 32 pp., 2012.
This report documents the DualPerm module [Šimůnek et al., 2003; Šimůnek and van Genuchten,
2008] for the HYDRUS (2D/3D) software package simulating two-dimensional variablysaturated water movement and solute transport in dual-permeability porous media. The dualpermeability flow and transport model assumes that soil consists of two mobile pore regions, one
representing the matrix and one the macropores. The dual-permeability formulation for waterflow is based on a mixed form of the Richards equation, describing water flow in both thefractures (macropores) and the matrix (micropores) domains [Gerke and van Genuchten, 1993a].
The dual-permeability formulation for solute transport is based on a convection-dispersion
equation, describing solute transport in both the fractures (macropores) and the matrix(micropores) domains [Gerke and van Genuchten, 1993a]. The mass transfer of water between
the two domains is driven by the gradient of pressure heads. The mass transfer for soluteincludes both convective mass transfer with water mass transfer, as well as diffusive masstransfer driven by the concentration gradient. Applications of the DualPerm module aredemonstrated on several examples. The DualPerm module is not intended for large transportdomains and should be used preferable for smaller domains.
This report serves as both a Technical and User Manual and is a reference document of theGraphical User Interface of the DualPerm related parts of the HYDRUS software package.

9

DISCLAIMER
This report documents the DualPerm module of HYDRUS (2D/3D). The DualPerm module wasdeveloped as a supplemental module of the HYDRUS software package, to model the twodimensional variably-saturated water movement and solute transport in dual-permeability porous
media. The software has been verified against selected test cases. However, no warranty is giventhat the program is completely error-free. If you do encounter problems with the code, find errors,or have suggestions for improvement, please contact one of the authors at
Jirka Šimůnek
Tel/Fax: 1-951-827-7854
Email: jiri.simunek@ucr.edu

10

1. Introduction
There is increasing evidence that flow and transport processes in soils often cannot be describedusing 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 insoils. Nonequilibrium water flow and solute transport in the unsaturated zone can be simulated atpresent 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., Šimůnek et al.,2008; Jacques and Šimůnek, 2005]. Numerical models are now increasingly being used since
they can be applied more readily to realistic laboratory and field problems as compared toanalytical models. The use of numerical models has been further popularized during the last 20
years or so because of the appearance of appropriate software packages in both the public andcommercial domains, and the development of increasingly sophisticated graphics-based
interfaces that can simplify their use tremendously.
Preferential flow in structured media (both macroporous soils and fractured rocks) can be
described using a variety of dual-porosity, dual-permeability, multi-porosity, and/or multipermeability models [Šimůnek et al., 2003; Šimůnek and van Genuchten, 2008] (Fig. 1). Dualporosity and dual-permeability models both assume that the porous medium consists of two
interacting regions, one associated with the inter-aggregate, macropore, or fracture system, andone comprising micropores (or intra-aggregate pores) inside soil aggregates or the rock matrix.
While dual-porosity models assume that water in the matrix is stagnant, dual-permeability
models allow for water flow in the matrix as well. While the uniform flow model and dualporosity models are available in the standard version of HYDRUS (2D/3D) (Fig. 1abc), dualpermeability models are not (Fig. 1d).

Figure 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 andimmobile flow regions, respectively; θM and θF in (d) are water contents of the matrix andmacropore (fracture) regions, respectively.

11

Dual-porosity models have long been applied to solute transport studies. Especially popular
early on were dual-porosity models in which distinct mobile and immobile flow regions are
assumed to be present. Dual-permeability models, in which water can move in both the interand intra-aggregate pore regions are now also becoming more popular. Available dualpermeability models differ mainly in how they implement water flow in and between the twopore regions, especially with respect to the degree of simplification and empiricism.
Approaches to calculating water flow in macropores or inter-aggregate pores range from those
invoking Poiseuille’s equation, the Green and Ampt or Philip infiltration models, the kinematicwave equation, and the Richards equation [Gerke and van Genuchten, 1993a].
We have implemented the dual-permeability module DualPerm (Fig. 1d) as a supplemental
module for HYDRUS (2D/3D) based on the approach suggested by Gerke and van Genuchten[1993a]. This report serves as both a Technical and User Manual and is a reference document of
the Graphical User Interface of the DualPerm-related parts of the HYDRUS software package
simulating two-dimensional variably-saturated water movement and solute transport in dualpermeability porous media. The numerical solution of the two-dimensional variably-saturated
water flow, and heat and solute movement are described in the HYDRUS (2D/3D) documentation,
and will not be repeated here. The DualPerm module is fully supported by the HYDRUS (2D/3D)
graphical user interface [Šejna et al., 2011].
The DualPerm module may be used to analyze water and solute movement in unsaturated, partiallysaturated, or fully saturated porous media. DualPerm can handle flow domains delineated byirregular boundaries. The flow region itself may be composed of nonuniform soils having an
arbitrary degree of local anisotropy. Flow and transport can occur in the vertical plane, thehorizontal plane, or in a three-dimensional region exhibiting radial symmetry about a vertical axis.
The water flow part of the model considers prescribed head and flux boundaries, as well asboundaries controlled by atmospheric conditions.
The governing flow and transport equations are solved numerically using standard Galerkin-typelinear finite element schemes. Applications of the DualPerm module are demonstrated on several
examples.

122. Dual-Permeability Water Flow
2.1. Governing Flow Equations
Different dual-permeability approaches (Fig. 1) may be used to describe flow and transport instructured media. While some models invoke similar equations for flow in the fracture andmatrix regions, others use different formulations for the two regions. A typical example of thefirst approach, implemented also in HYDRUS-1D, is the work of Gerke and van Genuchten[1993ab, 1996] who applied the Richard equation to each of the two pore regions. The flow
equations for the macropore (fracture) (subscript f) and matrix (subscript m) pore systems in theirapproach are given by∂θ f (h f )∂t= A ∂h fΓ∂ + K izA, f   - S f (h f ) - w
 K f (h f )  K ij , fw∂xi 
∂x j  A ∂h∂θ m (hm ) ∂ ΓA= K m (hm )  K ij ,m m + K iz ,m   - S m (hm ) + w∂t∂xi 
∂x j
1- w

 

(1)

where subscripts f and m refer to the fracture and matrix domains, w is the ratio of the volumes ofthe macropore or fracture domain and the total soil system [-], θ is the volumetric water content[L3L-3], h is the pressure head [L], S is a sink term [T-1], xi (i=1,2) are the spatial coordinates [L], t is
time [T], Γw is the transfer rate for water from the inter- to the intra-aggregate pores, KijA are
components of a dimensionless anisotropy tensor KA, and K is the unsaturated hydraulicconductivity function [LT-1] given by
K ( h, x , y , z ) = K s ( x , y , z ) K r ( h, x , y , z )(2)where Kr is the relative hydraulic conductivity and Ks the saturated hydraulic conductivity [LT-1].
The anisotropy tensor KijA in (1), which can be different for the two domains, is used to account foran anisotropic medium. The diagonal entries of KijA equal one and the off-diagonal entries zero for
an isotropic medium. If (1) is applied to planar flow in a vertical cross-section, x1=x is the horizontal
coordinate and x2=z is the vertical coordinate, the latter taken to be positive upward. Einstein'ssummation convention is used in (1) and throughout this report. Hence, when an index appears
twice in an algebraic term, this particular term must be summed over all possible values of the index.
Solutions of the Richards equation (1) require knowledge of the unsaturated soil hydraulicfunctions made up of the soil water retention curve, θ(h), which describes the relationshipbetween 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 θ. The dualpermeability approach, as developed by Gerke and van Genuchten [1993a], is relatively
complicated in that the model requires characterization of water retention and hydraulicconductivity functions (potentially of different form) for both pore regions, as well as a hydraulic13

conductivity function of the fracture-matrix interface. Note that the water contents θf and θm in(1) have different meanings than in the dual-porosity model where they represented watercontents of the total pore space (i.e., θ = θmo + θim), while here they refer to water contents of thetwo separate (fracture or matrix) pore domains such thatθ = θ F + θ M = wθ f + (1 − w)θ m(3)

Hence, lower case subscripts in the dual-permeability model refer to the local (pore-region) scale,
while upper case subscripts refer to the global (total soil medium) scale.
2.2. Water Mass Transfer
The rate of exchange of water between the fracture and matrix regions, Γw, is usually assumed to
be proportional to the difference in pressure heads between the two pore regions [Gerke and van
Genuchten, 1993a]:Γ w = α w ( h f - hm )(4)in which αw is a first-order mass transfer coefficient [L-1T-1]. Since pressure heads are now
needed for both regions, this approach requires retention curves for both pore regions. For porous
media with well-defined geometries, the first-order mass transfer coefficient, αw, can be defined
as follows [Gerke and van Genuchten, 1993b]:

αw =

β
d

2

K aγ w(5)

where d is an effective ‘diffusion’ pathlength (i.e. half the aggregate width or half the fracture
spacing) [L], β is a shape factor that depends on the geometry [-], and γw is a scaling factor(=0.4) obtained by matching the results of the first-order approach at the half-time level of the
cumulative infiltration curve to the numerical solution of the horizontal infiltration equation[Gerke and van Genuchten, 1993b]. The value of β varies with the geometry assumed for theaggregates [van Genuchten and Dalton, 1986]; it is 3 for rectangular slabs. Gerke and van
Genuchten [1996] suggested a method for deriving values of the shape factor β in (5) for nonidealized geometries or mixtures of shapes in order to obtain a more general 'macroscopic'
approach. Gerke and van Genuchten [1996] evaluated the effective hydraulic conductivity Ka[LT-1] of the fracture-matrix interface using a simple arithmetic average involving both hf and hm
as follows=
K a ( h ) 0.5  K a ( h f ) + K a ( hm ) 

(6)

The use of (6) implies that the medium contains geometrically well-defined rectangular or othertypes of macropores or fractures (e.g. Edwards et al. [1979], van Genuchten and Dalton [1986],
and Gerke and van Genuchten [1996]). While geometrically based models are conceptuallyattractive, they may be too difficult to use for field applications, partly because structured soils

14and rocks usually contain mixtures of aggregates and matrix blocks of various sizes and shapes,
but also because the parameters in (6) may not be identifiable. Hence, rather than using (6)directly, one could also lump β, d, and γw into one effective hydraulic conductivity Ka* of thefracture-matrix interface to give

α w = K a* ( h )(7)in which case Ka* can be used as a calibration parameter.

2.3. Soil Hydraulic Properties
The unsaturated soil hydraulic properties, θ(h) and K(h), in (1) are in general highly nonlinear
functions of the pressure head. The DualPerm module of HYDRUS implements the soil-hydraulicfunctions of van Genuchten [1980] who used the statistical pore-size distribution model of Mualem[1976] to obtain a predictive equation for the unsaturated hydraulic conductivity function in terms ofsoil water retention parameters. The expressions of van Genuchten [1980] are given byθs - θr θr +
θ ( h) = [1+ |α h| n] m
 θ
s


h<0

(8)

h≥0
2

K ( h ) = K s Sel [1 - (1 - Se1 / m ) m ](9)m = 1 - 1/n , n > 1

(10)

where

The above equations contain six independent parameters: θr, θs, α, n, Ks, and l. The poreconnectivity parameter l in the hydraulic conductivity function was estimated [Mualem, 1976] to beabout 0.5 as an average for many soils. These soil hydraulic functions are used for both the matrix
and fracture domains.

15

16

3. Dual-Permeability Solute Transport
3.1. Governing Solute Transport Equation
Analogous to Eq. (1), the dual-permeability formulation for solute transport is based on advectiondispersion type equations for transport in both the fracture and matrix regions as follows [Gerke andvan Genuchten, 1993a] (Fig. 1):∂θ f c f
∂t

+ρ

∂s f
∂t

=

∂c f  ∂q f ,i c fΓ∂ f-φf − s θ f Dij
 w
∂xi 
∂x j  ∂xi ∂qm ,i cm∂θ m cm
∂s
Γ
∂ 
m ∂c
+ρ m =- φm + s θ m Dij m  ∂t∂t∂xi ∂x j  ∂xi
1− w

(11)

where subscripts f and m again refer to the fracture and matrix domains, c and s are soluteconcentrations in the liquid [ML-3] and solid [MM-1] phases, respectively; q is the volumetricflux density [LT-1], ρ is the soil bulk density [M L-3], Dij is the dispersion coefficient tensor [L2T1] for the liquid phase, Γs is the mass transfer term for solute exchange between the fracture andmatrix domains [ML-3T-1], andφ represents the sink/source term, accounting for first-orderdegradation and zero-order production processes, as well as root solute uptake. The first equationof (11) describes solute transport in the fracture domain, the second equation transport in thematrix domain. Equation (11) assume complete advective-dispersive transport descriptions forboth the fractures and the matrix.
Similarly as in the standard solute transport model of HYDRUS (2D/3D), the solute transportmodel in the dual-permeability module can be linked with the chemical nonequilibrium two-site
sorption model or the physical nonequilibrium dual-porosity (mobile-immobile water) model.
For governing equations of these two nonequilibrium models, see the technical manual of
HYDRUS (2D/3D) [Šimůnek et al., 2011].
3.2. Solute Mass Transfer
The advective-dispersive mass transfer between the fracture and matrix domains is traditionallydescribed using the following form:Γ s = α s (1 − w)θ m (c f − cm ) + Γ w c *(12)where αs is the first-order solute mass transfer coefficient [T-1] accounting for physical rate
process. c* in (12) is equal to cf for Γw>0 and cm for Γw<0. van Genuchten and Dalton [1986]
and Gerke and van Genuchten [1996], among others, discussed possible expressions for the firstorder solute mass transfer coefficient, αs [T-1], which can be defined as:

17

αs =

β
d2

Da

(28)

in which Da is an effective diffusion coefficient [L2T-1], which represents the diffusion propertiesof the fracture-matrix interface as well as other parameters.

18

4. Preprocessing and Input Data
4.1. Soil Hydraulic Model
The Dual-Permeability module is activated when the "Dual-Permeability" check box is selectedin the Soil Hydraulic Model window (Fig. 2). Once this check box is selected, the edit box"Fraction of Boundary Flux into Fracture" appears. The value of this parameter can be
specified by a user and should vary between 0 and 1. When 0 is specified, all water at boundaries
with a specified boundary flux will enter into the matrix. When 1 is specified, all water atboundaries with a specified boundary flux will enter into the fracture domain.

Figure 2. The Soil Hydraulic Model dialog window.

4.2. Water Flow Parameters
The soil hydraulic parameters for the dual-permeability module are specified in the Water Flow
Parameters window (Fig. 3). Six soil hydraulic parameters (θr, θs, α, n, Ks, and l) are needed forboth the matrix and fracture domains to describe the retention curve and hydraulic conductivity
function. Additional 5 parameters are needed to characterize the matrix-fracture interface (w, β,γw, d, and Kas). Detailed description of all soil hydraulic parameters is given in Table 1.
19

Figure 3. The Water Flow Parameters dialog window (for three materials).

Table 1. Water flow parameters for the dual-permeability module.
Parameter Notation
Description
Qr
Residual soil water content of the matrix [L3/L3] (see eq. (8)).
θr
Qs
Saturated soil water content of the matrix [L3/L3].
θs
Alpha
Parameter α in the soil water retention function of the matrix [L-1].
αnnParameter n in the soil water retention function of the matrix [-].
Ks
Ks
Saturated hydraulic conductivity of the matrix [LT-1].
llTortuosity parameter in the conductivity function of the matrix [-].
QrFr
Residual soil water content of the fracture domain [L3/L3].
θr
QsFr
Saturated soil water content of the fracture domain [L3/L3].
θsAlphaFr
Parameter α in the soil water retention function of the fracture domain
α

20

nFr

n

KsFr
lFr

Ks
l

w

w

Beta

β

Gamma
a

γw

Ksa

Kas

d

[L-1].
Parameter n in the soil water retention function of the fracture domain[-].
Saturated hydraulic conductivity of the fracture domain [LT-1].
Tortuosity parameter in the conductivity function of the fracturedomain [-].
Ratio of the volumes of the macropore or fracture domain and the totalsoil system [-] (see eq. (5)). Must be the same for all soil layers.
Shape factor that depends on the geometry [-] (see eq. (5)); it is 3 forrectangular slabs.
Scaling factor (=0.4) [-] (see eq. (5)).
Effective ‘diffusion’ pathlength (i.e. half the aggregate width or half thefracture spacing) [L] (see eq. (5)).
Saturated hydraulic conductivity Kas [LT-1] of the fracture-matrix
interface (see eqs. (5) and (6)).

4.3. Solute Transport Parameters
The solute transport parameters for the dual-permeability module are specified in the Solute
Transport Parameters window (Fig. 4). Detailed description of these parameters is given in
Tables 2 and 3.

Figure 4. The Solute Transport Parameters dialog window.
a) Soil Specific Parameters
These parameters are common for all three components, but different for different soil layers.

21

Table 2. Soil specific solute transport parameters.
Parameter
Bulk. D.
DispL.M.
Disp.T.
Fract.M

Notation

ThImob.M.

θim

ρ
λL
λT

Description
3

Bulk density [M/L ]Longitudinal dispersivity in the matrix domain [L]Transverse dispersivity [L]Dimensionless fraction of adsorption sites classified as type-1 sites,
i.e., sites with instantaneous sorption when the chemical
nonequilibrium option is considered [-] in the matrix domain. Set thisparameter equal to 1 when equilibrium transport is considered. Frac
becomes the dimensionless fraction of adsorption sites in contact withmobile water when the physical nonequilibrium option is considered[-]. In that case Frac should be set equal to 1 when all sorption sites
are in contact with mobile water.
Immobile water content in the matrix domain. Set equal to 0 when thephysical nonequilibrium option is not considered.

fb) Solute Specific Parameters
These parameters are common for different soil layers, but may be different for different
components.
Table 3. Solute specific solute transport parameters.
Parameter
Diffus. W.
Mass.Tr.

Notation
Dw
Da

Description
Molecular diffusion in the liquid phase [LT-2]
Effective diffusion coefficient [L2T-1], which represents the diffusion
properties of the fracture-matrix interface as well as other parameters.

4.4. Solute Reaction Parameters
The solute transport parameters for the dual-permeability module are specified in the Solute
Reaction Parameters window (Fig. 5). Detailed description of these parameters is given in
Tables 4.
Table 4. Solute Reaction parameters.
Parameter
KdM
Nu

Notation
Kd

Beta

β

AlphaMωη

Description
Adsorption isotherm (distribution) coefficient for the matrix [M-1L3].
Adsorption isotherm (Langmuir) coefficient [L3/M] (used for bothmatrix and fracture domains).
Adsorption isotherm (Freundlich) exponent [-] (used for both matrix
and fracture domains).
First-order rate coefficient for nonequilibrium adsorption in the matrix

22

SinkL1µw

SinkS1µs

Disp.F.
SinkL1'λLµw'

SinkS1'µs'

Frac.F.

f

SinkL0
SinkS0
KdF

γw
γs
Kd

AlphaFω[T-1].
First-order degradation rate constant for dissolved phase [T-1] (usedfor both matrix and fracture domains).
First-order degradation rate constant for solid phase [T-1] (used forboth matrix and fracture domains).
Longitudinal dispersivity in the fracture domain [L]First-order degradation rate constant (representing the chain reaction)
for dissolved phase [T-1] (used for both matrix and fracture domains).
First-order degradation rate constant (representing the chain reaction)
for solid phase [T-1] (used for both matrix and fracture domains).
Dimensionless fraction of adsorption sites classified as type-1 sites,
i.e., sites with instantaneous sorption when the chemical
nonequilibrium option is considered [-] in the fracture domain. Set
this parameter equal to 1 when equilibrium transport is considered.
Frac becomes the dimensionless fraction of adsorption sites in contact
with mobile water when the physical nonequilibrium option isconsidered [-]. In that case Frac should be set equal to 1 when all
sorption sites are in contact with mobile water.
Zero-order (production) rate constant for dissolved phase [ML-3T-1].
Zero-order (production) rate constant for solid phase [T-1].
Adsorption isotherm (distribution) coefficient for the fracture domain[M-1L3].
First-order rate coefficient for nonequilibrium adsorption in thefracture domain [T-1].

Figure 5. Parts of the "Reaction Parameters for Solute -1" dialog window (for three materials).

234.5. Initial Conditions
Initial conditions for water flow can be specified either in terms of pressure heads or watercontents. When the pressure head is used to specify the initial conditions, the same pressure headis used in both matrix and fracture domains (i.e., it is assumed that the two domains are in
equilibrium). When the water content is used to specify the initial conditions, the specified watercontent is assigned to the matrix domain, and the corresponding pressure head (calculated fromthe retention curve of the matrix domain) is assigned to both the matrix and fracture domains.
The same initial concentrations are assumed for both the matrix and transport domains.

4.6. Anisotropy
In the standard HYDRUS code, the angle of anisotropy, and the 1st and 2nd component of theanisotropy tensor KijA are specified in GUI. In the DualPerm module, this input is reinterpreted forthe two (matrix and fracture) domains differently. The angle of anisotropy (Anisotropy - Angle;
Fig. 6) must be the same in both transport domains. The input "Anisotropy - 1st Component" (Fig.
6) is interpreted as a ratio of the horizontal and vertical components (KxA/KzA) of the anisotropy
tensor of the matrix domain, while the input "Anisotropy - 2st Component" (Fig. 6) is interpreted
as a ratio of the horizontal and vertical components (KxA/KzA) of the anisotropy tensor of the fracturedomain.

Figure 6. The "Domain Properties" part of the Data Tab of the Navigator Bar.

24

5. Post-processing
5.1. Results – Graphical Display
Multiple outputs from the DualPerm module are duplicated for the matrix and fracture domains(Table 5). Each of these outputs can be displayed as done in the standard HYDRUS version byclicking on corresponding commands in the "Results - Graphical Display" (Fig. 7) and "Results Other Information" parts of the Navigator bar. Matrix or fracture information is displayed
depending on the "Dual-Permeability Options" switch (Display Matrix or Display Fractures)that can be selected either in the "Results" part of the View Tab of the Navigator Bar (Fig. 8) oron the Edit (Tools) Bar (Fig. 9). In addition to the standard output, water and solute mass transfer
between the matrix and fracture domains can be also displayed (Fig. 7).
Table 5. Duplicated outputs from the DualPerm module for the matrix and fracture domains.
Fileh.out and hF.out
th.out and thF.out
v.out and vF.out
Conc1.out and
ConcF1.out
ObsNod.out and
ObsNodF.out
Solute1.out andSoluteF1.out

Description
Nodal values of pressure heads in the matrix and fracture domains.
Nodal values of water contents in the matrix and fracture domains.
Nodal values of fluxes in the matrix and fracture domains.
Nodal values of concentrations in the matrix and fracture domains.
Observation node values of pressure heads, water contents, andconcentrations in the matrix and fracture domains.
Boundary and reaction solute fluxes for the matrix and fracture domains.

Figure 7. The "Results - Graphical Display" part of the Data Tab of the Navigator Bar.

Figure 8. The "Results" part of the View Tab of the Navigator Bar.

25

Figure 9. The "Dual-Permeability Options" switch at the Edit (Tools) Bar.

5.2. The v_meanF.out Output File
An additional output file v_meanF.out is created, which is not recognized by GUI and which
provides information about various actual and cumulative boundary fluxes and average pressure
heads for the fracture domain (Table 6).
Table 6. The v_meanF.out output file.
Variablet
vFrac
CumFracvFracSvCumFracs
vMeanF(iKod)

CumQF(iKod)
hMeaF(iKod)

Description
Time
Water mass transfer between the fracture and matrix domain [L/T]
Cumulative water mass transfer between the fracture and matrix domain [L]Solute mass transfer between the fracture and matrix domain [M/L2/T]
Cumulative solute mass transfer between the fracture and matrix domain[M/L2]
Water flux across the boundary with Kode=iKod [L2/T] (iKod=1: constanthead or flux boundary; =2: seepage face; =3: variable head or flux; =4:
atmospheric; =5: drains; =6: free or deep drainage.
Cumulative water flux across the boundary with Kode=iKod [L2]
Average pressure head at the boundary with Kode=iKod [L]26

6. Example Problems
The DualPerm module supplements the standard HYDRUS program, which includes variablysaturated water flow and solute transport model, and thus the water flow and solute transport partsof the model have been widely tested earlier. Therefore, verification and test examples presented inthis section concentrate on the demonstration of the new features of the model.
The new module and most of its new features have been extensively tested on test problems ofincreasing complexity (Table 7). The first set of examples described in Sections 6.1 and 6.2 (Test1athrough Test1f; see Table 7) considered surface flux either into the fracture (Section 6.1) or matrix(Section 6.2) domain and different sizes of the aggregates (parameter d). The second test of
examples (Test2a through Test2f) considered surface ponding. The third test of examples (Test3
through Test3b) simulated tension disc infiltration experiment, considering either isotropic or
anisotropic fracture domain. The final two examples (DPerm_2SM and DPerm_MIM) considered,in addition to the dual-permeability medium, also the nonequilibrium problem in one (immobile
water in the matrix) or both (two-site sorption) domains.
Table 7. List of various test examples used to verify various new options of the DulPerm module.
Name
Test1aTest1a1
Test1a2
Test1a3
Test1a4
Test1b
Test1cTest1dTest1d1
Test1eTest1fTest2aTest2b
Test2cTest2dTest2eTest2f
Test3
Test3aSolute
No
No
No
Yes
No
No
No
Yes
Yes
Yes
Yes
No
No
No
Yes
Yes
Yes
Yes
Yes

Test3bYes

DPerm_2SM
DPerm_MIM

Yes
Yes

Description
Flux upper BC. Flux into fracture, d=1.0 cm.
Flux upper BC. Flux into matrix, d=1.0 cm.
Flux upper BC. Flux into matrix, d=1.0 cm. Higher flux.
Flux upper BC. Flux into matrix, d=1.0 cm. Solute.
Atmospheric upper BC. Root uptake, d=1.0 cm.
Flux upper BC. Flux into fracture, d=0.1 cm.
Flux upper BC. Flux into fracture, d=0.0316 cm.
Flux upper BC. Flux into fracture, d=1.0 cm. Solute.
Flux upper BC. Flux into fracture, d=1.0 cm. Solute inflow.
Flux upper BC. Flux into fracture, d=0.1 cm. Solute.
Flux upper BC. Flux into fracture, d=0.0316 cm. Solute.
Surface ponding, d=1.0 cm.
Surface ponding, d=0.2 cm.
Surface ponding, d=3.0 cm.
Surface ponding, d=1.0 cm, solute.
Surface ponding, d=0.2 cm, solute.
Surface ponding, d=3.0 cm, solute.
Ring infiltration into the isotropic dual-permeability system.
Ring infiltration into dual-permeability system; different anisotropy infracture (10:1) and matrix (1:1) domains.
Ring infiltration into dual-permeability system; different anisotropy in
Fracture (1:10) and Matrix (1:1) domains.
Dual-permeability model, two-site sorption model.
Dual-permeability model, dual-porosity (MIM) in matrix.

27

6.1. Example 1 (Test1d): Surface Flux Into FractureThis example demonstrates the application of the dual-permeability model of Gerke and van
Genuchten [1993a] to infiltration into a macroporous soil profile 40 cm deep. The following soil
hydraulic parameters for the matrix (θr=0.105, θs=0.50, α1=0.005 cm-1, n=1.5, l=0.5, Ks=1.05 cmd-1), the fracture (θr=0.0, θs=0.50, α1=0.1 cm-1, n=2.0, l=0.5, Ks=2000.0 cm d-1), and the matrixfracture interface (w=0.05, β=3, γ=0.4, a=1.0, and Kas=0.01 cm d-1, see Eq. (5) above) were used.
Water mass transfer was described using Eq. (4). Only convective solute mass transfer between
the two pore regions, with a dispersivity of 2 cm, was considered. The initial pressure head forboth regions was equal to –1000 cm, while water was applied only to the soil fracture.

Figure 10. Pressure heads in the matrix (a) and fracture (b) domains, water contents in the matrix(c) and fracture (d) domains, water mass exchange rates (e), solute mass exchange rates (f), andconcentrations in the matrix (h) and fracture (g) domains calculated with the dual-permeability
model of Gerke and van Genuchten [1993a,b].

28

6.2. Example 2 (Test1a3): Surface Flux Into Matrix
This example was taken from Šimůnek et al. [2003]. Figure 11 shows an application of the dualpermeability model of Gerke and van Genuchten [1993a] to infiltration into a macroporous soil
profile 40 cm deep. The same system and the same parameters were used as above (example 1)
except that water was applied only to the soil matrix until ponding occurred. The soil matrix at
the surface quickly became saturated since the applied water flux (50.0 cm d-1) greatly exceeded
the saturated hydraulic conductivity of the matrix (Ks = 1.05 cm d-1) (Fig. 11a). Water started toflow into the fracture domain (Fig. 10a) where it moved vertically downwards (Fig. 11c) andinfiltrated into the matrix domain (Fig. 10d). The two distinct moisture fronts in Figure 11b (andconcentration fronts in Figure 11e) are due to lateral transfer from the fracture domain (theleading edge), and flow in the matrix domain (the second front).

Figure 11. Infiltration and mass exchange fluxes (a), water contents in the matrix (b) and fracture(c) domains, water mass exchange rates (d) and concentrations in the matrix (e) and fracture (f)domains calculated with the dual-permeability model of Gerke and van Genuchten [1993a,b].

29

6.3. Examples 3 (Test3) - Tension Disc Infiltration Experiment
The tension disc infiltration experiment (with a disc radius of 10 cm) is simulated in this test
example. The axisymmetrical transport domain is 50 cm wide and 150 cm deep. The soil profile isassumed to consist of three soil layers (with each upper layers 20 cm thick). The soil hydraulic andsolute transport and reaction parameters are shown above in Figures 3, 4, and 5. The solute free,equilibrium initial profile with groundwater table at the bottom of the soil profile is considered.
The tension applied at the disc is -0.4 cm. The soil matrix is considered to be isotropic, while thefracture domain is assumed to be either isotropic, or to have the ratio KxA/KzA equal to 10 or 0.1.
Pressure head profiles for these three conditions with different anisotropies of the fracture domain
are shown in Figure 12. Note that while we shoe in Figure 12a the pressure head profile in thematrix only for the first scenario (i.e., both domains are isotropic), the matrix pressure headprofiles for the other scenarios are not very different.

Figure 12. Pressure head profiles for the matrix, isotropic fracture, and fracture with KxA/KzA=10,
and fracture with KxA/KzA=0.1 after 19 minutes.

30

References
Edwards, W. M., R. R., van der Ploeg, and W. Ehlers, A numerical study of the effects ofnoncapillary-sized pores upon infiltration, Soil Sci. Soc. Am. J., 43, 851-856, 1979.
Gärdenäs, A., J. Šimůnek, N. Jarvis, and M. Th. van Genuchten, Two-dimensional modelling ofpreferential water flow and pesticide transport from a tile-drained field, J. Hydrology, 329,
647-660, 2006.
Gerke, H. H., and M. Th. van Genuchten, A dual-porosity model for simulating the preferential
movement of water and solutes in structured porous media, Water Resour. Res., 29, 305-319,
1993a.
Gerke, H. H., and M. Th. van Genuchten, Evaluation of a first-order water transfer term forvariably saturated dual-porosity flow models, Water Resour. Res., 29, 1225-1238, 1993b.
Gerke, H. H., and M. Th. van Genuchten, Macroscopic representation of structural geometry for
simulating water and solute movement in dual-porosity media, Adv. Water Resour., 19, 343357.
Hendrickx, J. M. H., and M. Flury, Uniform and preferential flow, Mechanisms in the vadose zone.
In: Conceptual Models of Flow and Transport in the Fractured Vadose zone. National ResearchCouncil, National Academy Press, Washington, D.C., pp. 149-187, 2001.
Jacques, D., and J. Šimůnek, User Manual of the Multicomponent Variably-Saturated Flow andTransport Model HP1, Description, Verification and Examples. Version 1.0., SCK•CENBLG-998, Waste and Disposal, SCK•CEN, Mol, Belgium, 79 pp., 2005.
Köhne, J. M., S. Köhne, and J. Šimůnek, Multi-process herbicide transport in structured soilcolumns: Experiment and model analysis, J. Contam. Hydrology, 85, 1-32, 2006.
Mualem, Y., A new model for predicting the hydraulic conductivity of unsaturated porous media,
Water Resour. Res., 12, 513-522, 1976.
Nkedi-Kizza, P., J. W. Biggar, H. M. Selim, M. Th. van Genuchten, and P. J. Wierenga,
Equivalence of two conceptual models for describing Ion exchange during transport throughan aggregated Oxisol, Water Resour. Res., 20(8), 1123-1130, 1984.
Pot, V., J. Šimůnek, P. Benoit, Y. Coquet, A. Yra and M.-J. Martínez-Cordón, Impact of rainfall
intensity on the transport of two herbicides in undisturbed grassed filter strip soil cores, J. of
Contaminant Hydrology, 81, 63-88, 2005.
Šejna, M., J. Šimůnek, and M. Th. van Genuchten, The HYDRUS Software Package for Simulating
Two- and Three-Dimensional Movement of Water, Heat, and Multiple Solutes in VariablySaturated Media, User Manual, Version 2.0, PC Progress, Prague, Czech Republic, pp. 280,2011.
Šimůnek, J., N. J. Jarvis, M. Th. van Genuchten, and A. Gärdenäs, Review and comparison ofmodels for describing non-equilibrium and preferential flow and transport in the vadose zone,
Journal of Hydrology, 272, 14-35, 2003.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, Modeling Subsurface Water Flow and Solute
Transport with HYDRUS and Related Numerical Software Packages, In: Garcia-Navarro &Playán (eds.), Numerical Modelling of Hydrodynamics for Water Resources, AnInternational Workshop, Centro Politecnico Superior, University of Zaragoza Spain, June 1821 2007. Taylor & Francis Group, London, ISBN 978-0-415-44056-1, 95-114, 2007.

31Šimůnek, J., and M. Th. van Genuchten, Modeling nonequilibrium flow and transport with
HYDRUS, Vadose Zone Journal, doi:10.2136/VZJ2007.0074, Special Issue “Vadose Zone
Modeling”, 7(2), 782-797, 2008.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, Development and applications of the
HYDRUS and STANMOD software packages and related codes, Vadose Zone J 7(2), 587600, 2008.
Šimůnek, J., M. Th. van Genuchten, and M. Šejna, The HYDRUS Software Package for Simulating
Two- and Three-Dimensional Movement of Water, Heat, and Multiple Solutes in VariablySaturated Media, Technical Manual, Version 2.0, PC Progress, Prague, Czech Republic, pp. 258,
2011.
Toride, N., F. J. Leij, and M. Th. van Genuchten, A comprehensive set of analytical solutions fornonequilibrium solute transport with first-order decay and zero-order production, Water Resour.
Res., 29(7), 2167-2182, 1993.
van Genuchten, M. Th., A closed-form equation for predicting the hydraulic conductivity ofunsaturated soils, Soil Sci. Soc. Am. J., 44, 892-898, 1980.
van Genuchten, M. Th., Non-equilibrium transport parameters from miscible displacement
experiments, Research Report No. 119, U.S. Salinity Laboratory, Riverside, CA, 1981.
van Genuchten, M. Th., and F. N. Dalton, Models for simulating salt movement in aggregatedfield soils, Geoderma, 38, 165-183, 1986.
van Genuchten, M. Th., and J. Šimůnek, Integrated modeling of vadose zone flow and transport
processes, Proc. Unsaturated Zone Modelling: Progress, Challenges and Applications, Eds. R.
A. Feddes, G. H. de Rooij, and J. C. van Dam, Wageningen UR Frontis Series, Vol. 6,
Chapter 2, pp. 37- 69, x-xi, Kluwer Academic Publishers, Dordrecht, The Netherlands, 2004.

32