Published online September 13, 2005

Modeling Variably Saturated Water Flow and Multicomponent
Reactive Transport in Constructed WetlandsGünter Langergraber* and Jirka Šimůnek

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

ABSTRACT

important factors in controlling subsurface flow and/or
contaminant transport, and/or evaluating regional andglobal water cycles (Scanlon et al., 2002).
Many or most vadose zone models consider the transport of only one solute and assume that the fate andtransport of this solute is independent of all other species that may be present in the soil solution. In reality,
the soil solution is always a mixture of many chemical
species and microorganisms that may mutually interact,
create complexed species, precipitate or dissolve, affecteach others degradation and decay, and/or compete witheach other for sorption sites (van Genuchten and Šimůnek, 2004). Many important environmental problems
hence require a simultaneous analysis of prevailing flow
and transport processes with multiple chemical and biological reactions and transformations. Examples involving mostly geochemical reactions as listed by Šimůnek
and Valocchi (2002) include acid mine drainage (Walter
et al., 1994; Lichtner, 1996), radionuclide transport (Viswanathan et al., 1998), and reactive permeable barriers
for aquifer remediation (Fryar and Schwartz, 1994). Examples that additionally involve biochemical processes
include degradation of NTA (Tebes-Stevens et al., 1998),
fate and transport of metal–organic mixed wastes (Rittmann and VanBriesen, 1996; VanBriesen, 1998), theformation of redox zones in organic-contaminated aquifers (Abrams et al., 1998; Essaid et al., 1995), C and N
cycles in soils (Parton et al., 1988), bacteria-inducedchanges in the soil hydraulic properties (Rockhold et al.,2002), and complex biogeochemical processes in CWs(Langergraber, 2001, 2003).
Constructed wetlands are engineering structures usedworldwide to improve water quality (Kadlec et al., 2000;
Langergraber and Haberl, 2001; Haberl et al., 2003).
They involve a complex mixture of water, substrate,
plants, litter, and a variety of microorganisms to produce
the optimal conditions for removing OM, N, P, andfor decreasing the concentrations of toxic trace metals,
organic chemicals, and pathogens. Constructed wetlands
include surface flow and subsurface flow CWs (Kadlec andKnight, 1996). Surface-flow CWs are generally densely
vegetated and typically have water depths of ⬍0.4 m.
In subsurface-flow CWs no free water level is visible.
Subsurface vertical flow CWs with intermittent feeding
are now used widely due to their efficiency in removing
ammonia N (e.g., Langergraber and Haberl, 2001). Al-

Constructed wetlands (CWs) are becoming increasingly popularworldwide for removing organic matter (OM), nutrients, trace elements, pathogens, or other pollutants from wastewater and/or runoff
water. We present a multicomponent reactive transport model CW2D
(i.e., Constructed Wetlands 2D), as an extension of the HYDRUS-2D
variably saturated water flow and solute transport software package.
CW2D was developed to model the biochemical transformation anddegradation processes in subsurface-flow CWs. Such wetlands involvea complex mixture of water, substrate, plants, litter, and a variety of
microorganisms to provide optimal conditions for improving water
quality. The water flow regime in subsurface-flow CWs can be highly
dynamic and requires the use of a transient variably saturated flow
model. The biochemical components defined in CW2D include dissolved oxygen (DO), three fractions of OM (readily and slowly biodegradable, and inert), four N compounds (ammonium, nitrite, nitrate,and dinitrogen), inorganic P, and heterotrophic and autotrophic microorganisms. Organic N and organic P were modeled as part of the
OM. The biochemical degradation and transformation processes werebased on Monod-type rate expressions. All process rates and diffusion
coefficients were assumed to be temperature dependent. Heterotrophic bacteria were assumed to be responsible for hydrolysis, mineralization of OM (aerobic growth), and denitrification (anoxic growth).
Autotrophic bacteria were assumed to be responsible for nitrification,
which was modeled as a two-step process. Lysis was considered tobe the sum of all decay and sink processes. We demonstrate theperformance of the model for one- and two-stage subsurface verticalflow CWs. Model simulations of water flow, tracer transport, andselected biochemical compounds are compared with experimental
observations. Limitations of the model are discussed, and needs formodel improvements are summarized.

O

ur ability to model flow and transport processes
in the vadose zone between the soil surface andthe groundwater table has increased enormously during
the past several decades (van Genuchten and Šimůnek,2004). Vadose zone models are now routinely used fora large number of applications, both in research of subsurface processes and management of our soil and water
resources. Typical applications are evaluation of theperformance and effectiveness of engineered covers tominimize infiltration into underlying waste, quantifying
groundwater recharge, evaluating the impact of climate
and land-use changes on subsurface flow, identifying

G. Langergraber, Institute of Sanitary Engineering and Water Pollution Control, BOKU Univ. of Natural Resources and Applied Life
Sciences, Vienna, Muthgasse 18, A-1190 Vienna, Austria; J. Šimůnek,
Dep. of Environmental Sciences, Univ. of California, Riverside, CA
92521, USA. Received 25 Nov. 2004. *Corresponding author (guenter.
langergraber@boku.ac.at).

Abbreviations: ASM, Activated Sludge Model; COD, chemical oxygen demand; CW, constructed wetland; CW2D, Constructed Wetlands
2D; DO, dissolved oxygen; IP, inorganic phosphorus; OM, organic matter; OM CI, inert OM; OM CR, readily biodegradable OM; OM CS,
slowly biodegradable OM; PSCW, pilot-scale subsurface vertical flow
constructed wetland; SSP, Small-Scale Plot; TOC, total organic carbon;
XANb, autotrophic microorganisms, Nitrobacter; XANs, autotrophic
microorganisms, Nitrosomonas; XH, heterotrophic microorganisms.

Published in Vadose Zone Journal 4:924–938 (2005).
Original Researchdoi:10.2136/vzj2004.0166© Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

924

925

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

though much experience exists in constructing and operating such systems, their design is still based mostly
on empirical rules, such as requiring a certain area perperson or using simple first-order kinetic rates (e.g.,
Kadlec and Knight, 1996). However, Kadlec (2000) summarizes evidence that first-order models are inadequate
for the design of treatment wetlands.
Understanding how CWs function is difficult due tothe large number of physical, chemical, and biologicalprocesses that are active at the same time and mutually
influence each other. Numerical simulators can represent valuable tools for analyzing and improving ourunderstanding of complex systems such as CWs. Only
few simulation results describing the behavior of subsurface flow CWs have been published. Schwager and Boller(1997) simulated tracer experiments and O2 transport
in intermittent sand filters using HYDRUS (an older
version of HYDRUS-1D) and MOFAT programs, respectively. However, they did not consider reactive
transport. Other models focused mainly on a descriptionof seasonal trends (e.g., Wynn and Liehr, 2001) or used
simplifications such as first-order reaction rates to describe the degradation of a substance along the flow
path (Werner and Kadlec, 2000). A detailed reaction
model using a similar approach as described in this paper
was used by Tanner et al. (1999) to model reactions insurface flow CWs.
To our knowledge, CW2D is the first model that describes detailed reactions in subsurface flow CWs in
combination with a flow model capable of describing
variably saturated flow. Our main objectives in this article are (i) to present the development of a multicomponent reactive transport model, CW2D (i.e., ConstructedWetlands 2D) (Langergraber, 2001), which is an extension of the variably-saturated water flow and solutetransport program HYDRUS-2D and which preserves
all numerical and graphical capabilities of HYDRUS-2D(Šimůnek et al., 1999); (ii) to demonstrate the performance of the model for one- and two-stage subsurface
flow CWs; and (iii) to compare model simulations ofwater flow, tracer transport, and selected biochemicalcompounds against experimental observations.

THEORY
Variably Saturated Flow
The water flow regime in subsurface-flow CWs canbe highly dynamic because of intermittent loadings ofwastewater, and hence requires the use of a transientvariably saturated flow model. The computational module of the HYDRUS-2D software package numericallysimulates two-dimensional variably saturated water flowusing the Richards equation (Šimůnek et al., 1999). Although HYDRUS-2D offers various models for the soil
hydraulic properties, in this study we used only the formulation proposed by van Genuchten (1980). Thismodel has six parameters, namely the residual and saturated water contents, ␪r and ␪s, respectively; the saturated hydraulic conductivity, Ks, the van Genuchten
shape parameters ␣ and n, and Mualem’s pore connec-

tivity parameter l. The HYDRUS-2D model also offers
a wide range of both system-dependent and independent boundary conditions. For common applications to
CWs, two additional boundary conditions had to be
implemented. One concerns surface ponding during
wastewater loadings that exceed the infiltration capacityof a vertical flow CW, which typically occurs during
its operation. Surface ponding was implemented as anextension of the atmospheric boundary condition. For
applications to CWs to treat combined sewer overflow,
it was also necessary to implement a limited outflow
boundary condition at the bottom of the CW to enablemodeling of the reduction of peak flows. Reduction ofpeak flows is used to increase the treatment performance and to reduce the impact of combined sewer
overflow on the receiving water (Dittmer et al., 2004).

Solute Transport
The governing equation for the macroscopic transportof component i can be written in the form
⳵␪ci ⳵␳si⫹⫽ ⵜ(␪Diⵜci) ⫺ ⵜ(qⵜci) ⫹ Scs,i ⫹ ri
⳵t
⳵t

[1]

where i ⫽ 1,..., N (N is the number of components), ci
is the concentration in the aqueous phase (M L⫺3), si isthe concentration in the solid phase (M M⫺1), ␪ is thevolumetric water content (L3 L⫺3), ␳ is the soil bulk
density (M L⫺3), Di is the effective dispersion tensor (L2T⫺1), q is the volumetric flux density (L3 L⫺2 T⫺1), S is thesource–sink term (L3 L⫺3 T⫺1); cs,i is the concentration of
the source–sink (M L⫺3), and ri is the reaction term (ML⫺3 T⫺1). Components of the effective dispersion tensor
Di include molecular diffusion, and longitudinal andtransverse dispersion. The reaction term ri in the multicomponent reactive transport model CW2D is describedin detail below. Liquid- and solid-phase concentrations
at equilibrium can be related using either linear, Frendlich, or Langmuir adsorption isotherms (Šimůnek etal., 1999).
The concept of two-region, dual-porosity type solutetransport (van Genuchten and Wierenga, 1976) is implemented in HYDRUS-2D to permit consideration ofphysical nonequilibrium transport. The physical nonequilibrium transport model divides the liquid phase
into mobile (flowing) and immobile (stagnant) regions.
Solute exchange between the mobile and immobile regions is modeled as a first-order process.

The CW2D Multicomponent Reactive
Transport Module
The need for a numerical model describing biochemical transformation and degradation processes in subsurface flow CWs motivated the development of the CW2D
multicomponent reactive transport module (Langergraber, 2001) for incorporation into the HYDRUS-2D
variably saturated water flow and solute transport program. CW2D is able to model the biochemical transformation and degradation processes for OM, N, and P.
The reactive transport module CW2D considers 12 components and nine processes. The components include

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

926

VADOSE ZONE J., VOL. 4, NOVEMBER 2005

Table 1. Components defined in the reactive transport module.

Table 2. Kinetic parameters (at 20ⴗC, values in brackets at 10ⴗC).†i

Parameter

1
2
3
4
5
6
7
8
9
10
11
12

Symbol

Unit

O2
CR
CS
CI
XH
XANs
XANb
NH4N
NO2N
NO3N
N2N
IP

L⫺1

mgO2
mgCOD L⫺1
mgCOD L⫺1
mgCOD L⫺1
mgCOD L⫺1
mgCOD L⫺1
mgCOD L⫺1
mgN L⫺1
mgN L⫺1
mgN L⫺1
mgN L⫺1
mgP L⫺1

Description
dissolved oxygen
readily biodegradable chemical oxygen demandslowly biodegradable chemical oxygen demand
inert chemical oxygen demandheterotrophic microorganisms
Nitrosomonas spp. (autotrophic bacteria 1)
Nitrobacter spp. (autotrophic bacteria 2)
Ammonium, NH4ⴙNitrite, NO2⫺
Nitrate, NO3⫺
Dinitrogen gas, N2
Inorganic phosphorus

DO, OM (three fractions of different degradability),ammonium, nitrite, nitrate, and N2 gas, inorganic P, andheterotrophic and autotrophic microorganisms. Organic
N and organic P were modeled as nutrient contentsof the OM. The processes considered are hydrolysis,
mineralization of OM, nitrification (modeled as a twostep process), denitrification, and a lysis process forthe microorganisms.
The mathematical structure of CW2D was based on
the structure of the Activated Sludge Models (ASMs)
introduced by Henze et al. (2000). Major assumptions
in all ASMs are (a) a constant value of pH, (b) constantcoefficients in the rate equations, and (c) constant stoichiometric factors. The various ASMs are valid only
for domestic wastewater and are applicable within thetemperature range of 10 to 25⬚C (Henze et al., 2000).
ASMs assume that nitrification occurs in one step, i.e.,that ammonia is directly converted (nitrified) into nitrate. The two-step nitrification model used in CW2D
(i.e., sequential nitrification of ammonium into nitrite
and nitrate) is based on models presented by Nowak(1996) and Brouwer et al. (1998).

Hydrolysis
hydrolysis rate constant, 1 d⫺1saturation/inhibition coefficient for hydrolysis,
mgCOD,CS mgCOD,BM⫺1
Heterotrophic microorganisms (aerobic growth)␮Hmaximum aerobic growth rate on CR, 1 d⫺1rate constant for lysis, 1 d⫺1bH
Khet,O2
saturation/inhibition coefficient for O2, mgO2 L⫺1saturation/inhibition coefficient for substrate,
Khet,CR
mgCOD,CR L⫺1
Khet,NH4N
saturation/inhibition coefficient for NH4ⴙ(nutrient), mgNH4N L⫺1
Khet,IPsaturation/inhibition coefficient for P, mgIP L⫺1
Heterotrophic microorganisms (denitrification)
␮DN
maximum denitrification rate, 1 d⫺1
KDN,O2saturation/inhibition coefficient for O2, mgO2 L⫺1
KDN,NO3N saturation/inhibition coefficient for NO3⫺,
mgNO3N L⫺1
KDN,NO2N saturation/inhibition coefficient for NO2⫺,
mgNO2N L⫺1
KDN,CRsaturation/inhibition coefficient for substrate,
mgCOD,CR L⫺1
KDN,NH4N saturation/inhibition coefficient for NH4ⴙ(nutrient), mgNH4N L⫺1
KDN,IPsaturation/inhibition coefficient for P, mgIP L⫺1
Nitrification step 1 (Nitrosomonas)
␮ANs‡maximum aerobic growth rate on NH4ⴙ, 1 d⫺1rate constant for lysis, 1 d⫺1
bANs‡
KANs,O2‡
saturation/inhibition coefficient for O2, mgO2 L⫺1
KANs,NH4N‡ saturation/inhibition coefficient for NH4ⴙ,
mgNH4N L⫺1
KANs,IPsaturation/inhibition coefficient for P, mgIP L⫺1
Nitrification step 2 (Nitrobacter)␮ANb‡maximum aerobic growth rate on NO2⫺, 1 d⫺1bANb‡
rate constant for lysis, 1 d⫺1
KANb,O2‡ saturation/inhibition coefficient for O2, mgO2 L⫺1
KANb,NO2N‡ saturation/inhibition coefficient for NO2⫺,
mgNO2N L⫺1
KANb,NH4N saturation/inhibition coefficient for NH4ⴙ(nutrient), mgNH4N L⫺1
KANb,IPsaturation/inhibition coefficient for P, mgIP L⫺1
Kh
KX

Value
3 (2)
0.1

6 (3)
0.4 (0.2)
0.2
2
0.05
0.01
4.8 (2.4)
0.2
0.5
0.5
4
0.05
0.01
0.9 (0.3)
0.15 (0.05)
1
0.5
0.01
1 (0.35)
0.15 (0.05)
0.1
0.1
0.05
0.01

† Values: Henze et al. (2000), except where ‡ indicates Nowak (1996).

Components
The components defined in CW2D (Table 1) are:• dissolved oxygen, O2 (mgO2 L⫺1)
• organic matter, OM (mgCOD L⫺1)—three pools, expressed in units of chemical oxygen demand(COD), are considered: readily biodegradable OM(CR), slowly biodegradable OM (CS), and inert
OM (CI)• nitrogen (mgN L⫺1) (including ammonium, NH4⫹;
nitrite, NO2⫺; nitrate, NO3⫺; and dinitrogen gas N2)• inorganic phosphorus, IP (mgP L⫺1)
• heterotrophic microorganisms, XH (mgCOD L⫺1)
• autotrophic microorganisms (Nitrosomonas, XANs,
and Nitrobacter, XANb) (mgCOD L⫺1)
Organic N and organic P were modeled as N and P
contents of the COD using the composition parameters
as given below. We assumed that the OM was present
only in the aqueous phase and that all reactions occurred
only in the aqueous phase as well. Adsorption was considered for ammonium, NH4⫹, and IP. Adsorption was
assumed to be a kinetic process.
Biofilm in CW2D was assumed to be ideal, composedof a homogeneous matrix of microorganisms, with all
species uniformly distributed over the biofilm. All micro-organisms were assumed to be immobile and hence wereassociated exclusively with the solid phase. Lysis wasconsidered to be the sum of all decay and sink processesof all microorganisms involved, with the rate of lysis
being independent of environmental conditions.
Reactions
The biochemical degradation and transformation processes in CW2D were modeled using Monod-type rate
expressions as used in the ASMs (Henze et al., 2000).
Equations [4] through [12] give the reaction rates rcj
for the nine processes considered. Parameters in theseequations are defined in Table 2. The various rate parameters in the kinetic growth and decay expressions wereassumed to be temperature dependent. Values in Table 2
are given for a reference temperature of 20⬚C and fortemperature dependent parameters also for 10⬚C. The
temperature dependence of all process rates and diffusion coefficients was assumed to be described using theArrhenius equation.
Two factors are used in Eq. [4] through [12] as anabbreviation. The first factor describes the conversionof solid-phase concentrations into liquid-phase concentrations:927

www.vadosezonejournal.org

cXY ⫽

␳
sXY
␪

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

where the subscripts Y ⫽ H, ANs, and ANb are used
for the heterotrophic and two autotrophic organisms,
respectively. The second factor is the Monod term fornutrients:
fN,x ⫽

cNH4
CIP
Kx,NH4 ⫹ cNH4 Kx,IP ⫹ CIP[3]

where the subscript x for Het, DN, and ANb represents
the processes of heterotrophic growth, denitrificationand growth of Nitrobacter, respectively.
In this study we assumed heterotrophic bacteria tobe responsible for hydrolysis, mineralization of OM(aerobic growth), and denitrification (anoxic growth).
Hydrolysis. Hydrolysis describes the conversion ofslowly biodegradable OM CS into readily biodegradable
OM CR, with a small fraction being converted into inert
OM CI. Ammonium (NH4⫹) and IP are released during
this transformation process. We further assume that hydrolysis takes place independently of the O2 conditions,leading to the following rate equation for hydrolysis:
rc1 ⫽ KhcCS/cXH
cXH
KX ⫹ cCS/cXH[4]

Aerobic Growth of Heterotrophic Bacteria. This process consumes O2 and OM CR, while NH4⫹ and IP areincorporated in the biomass as follows:
cO2
cCR
fN,Het cXH
rc2 ⫽ ␮H
KHet,O2 ⫹ cO2 KHet,CR ⫹ cCR[5]

Nitrate-Based Growth of Heterotrophs on Readily Biodegradable COD (denitrification). This process consumes NO3⫺ and OM CR. Nitrate is reduced to dinitrogen (N2). Again, NH4⫹ and IP are incorporated in thebiomass as follows:
rc3 ⫽ ␮DN

KDN,O2cNO3KDN,O2 ⫹ cO2 KDN,NO3 ⫹ cNO3
cCR
KDN,NO2
fN,DN cXH
KDN,NO2 ⫹ cNO2 KDN,CR ⫹ cCR[6]

Inhibition due to NO3⫺–N is modeled as a noncompetitive process using a soluble nonbiodegradable inhibitor(Nowak et al., 1995).
Nitrite-Based Growth of Heterotrophs on Readily
Biodegradable COD (denitrification). Similarly as forNO3⫺, this process consumes NO2⫺, OM CR, NH4⫹, and
IP. Nitrite is reduced to N2. The inclusion of both denitrification processes is necessary because of the use of a
two-step nitrification model.
rc4 ⫽ ␮DN

[8]

Autotrophic bacteria are responsible for nitrification,
which is modeled as a two-step process. Nitrosomonas
spp. consume NH4⫹ and produce NO2⫺, while Nitrobacter
spp. transform NO2⫺ to NO3⫺. Both steps of the nitrification process are strictly aerobic and can therefore only
occur when DO is available.
Aerobic Growth of Nitrosomonas on Ammonium (First
Step of Nitrification). This process consumes NH4⫹ andO2 and produces NO2⫺. Inorganic P and a small portionof NH4⫹ are incorporated in the biomass.
rc6 ⫽ ␮ANs

cO2
cNH4
KANs,O2 ⫹ cO2 KANs,NH4 ⫹ cNH4
cIP
cXANs
KANs,IP ⫹ cIP[9]

Lysis of Nitrosomonas. This lysis produces OM (CS,CR, and CI), NH4⫹, and IP as follows:
rc7 ⫽ bXANscXANs[10]

Aerobic growth of Nitrobacter on Nitrite (Second Stepof Nitrification). This process consumes NO2⫺ and O2
and produces NO3⫺. Ammonium and IP are incorporated in the biomass as follows:
rc8 ⫽ ␮ANb

cO2
cNO2
fN,ANb cXANb
KANb,O2 ⫹ cO2 KANb,NO2 ⫹ cNO2

[11]

Lysis of Nitrobacter. This process produces OM (CS,CR, and CI), NH4⫹, and IP as follows:
rc9 ⫽ bXANbcXANb[12]

Stoichiometric Matrix
Table 3 represents the stoichiometric matrix for
CW2D. The concept of yield coefficients is used to describe biomass growth. Yield coefficients represent theamount of substrate that is used for biomass growth.
Stoichiometric coefficients (including the yield coefficients) that appear in Table 3 are defined, together with
their default values, in Table 4 (for OM and microorganisms), Table 5 (for NH4⫹), and Table 6 (for IP).
Table 4 defines the stoichiometric and composition
parameters used in Table 3. The composition parameters provide information about the N and P contents of
the OM and biomass, and therefore are measures of themass of organic N and organic P, respectively.
Stoichiometric parameters for NH4⫹–N, ␷1,N through␷9,N, and IP, ␷1,P through ␷9,P, in Table 3 can be calculated
from the mass balance equations. These mass balance
equations can be written for each process j as follows:N兺 ␷j,i ic,i ⫽ 0

KDN,O2cNO2KDN,O2 ⫹ cO2 KDN,NO2 ⫹ cNO2
cCR
fN,DN cXH
KDN,CR ⫹ cCR

rc5 ⫽ bHcXH

[2]

[13]

i⫽1

[7]

Lysis of Heterotrophic Organisms. Lysis produces OM(CS, CR, and CI), NH4⫹, and IP. Lysis is assumed torepresent the sum of all decay and sink processes.

where i ⫽ 1,..., N (N is the number of components),j ⫽ 1,..., R (R is the number of processes), ␷j,i is the
stoichiometric factor for component i and process j (mgimgj⫺1), c represents the material to which the mass balance is applied (either N or P), and ic,i is a conversion
factor that converts the unit of component i to the unit

928

␷8,P
␷9,P

␷7,P

␷6,P

␷5,P

␷4,P

␷3,P

(1 ⫺ YH)/
2.86YH(1 ⫺ YH)/
1.71YH

␷1,P
␷2,P

gP

Table 4. Stoichiometric parameters for organic matter and microorganisms and composition parameters.
Parameter⫺(1 ⫺ YH)/
2.86YH

1/YANb
⫺1/YANb

YH
YANs*
YANb*
iN,CR
iN,CS
iN,CI
iN,BM
iP,CR
iP,CS
iP,CI
iP,BM

␷8,N
␷9,N

␷7,N

1/YANs⫺1/YANs
iN,BM

␷5,N

⫺(1 ⫺ YH)/
1.71YH
␷3,N

0.1
0.63
0.24
0.24
0.03
0.04
0.01
0.07
0.01
0.01
0.01
0.02

fBM,CI

It is assumed that reactions take place only in the
aqueous phase. Therefore the overall reaction term ri
for component i can be calculated as follows:Rri ⫽ ␪ 兺 ␷j,i rcj

1 ⫺ fBM,CR ⫺fBM,CI[14]

j⫽1† See Table 4, Table 5, and Table 6 for definitions of the stoichiometric coefficients.

fBM,CR
Aerobic growth of Nitrobacter on nitrite
Lysis of Nitrobacter89⫺(1.14 ⫺ YANb)/
YANb

Aerobic growth of Nitrosomonas onammonium
Lysis of Nitrosomonas
7

0
0.05

Calculation of the Reaction Term

6

5

4

3

Stoichiometric parameters
Production of CI in hydrolysis, mgCOD,CI mgCOD,CS⫺1
Fraction of CR generated in biomass lysis, mgCOD,CR
mgCOD,BM⫺1
Fraction of CI generated in biomass lysis, mgCOD,CI
mgCOD,BM⫺1
Yield coefficient for heterotrophs, mgCOD,BMmgCOD,CR⫺1
Yield coefficient for N. somonas, mgCOD,BM mgNH4N⫺1
Yield coefficient for N. bacter, mgCOD,BM mgNO2N⫺1
Composition parameters
N content of CR, mgN mgCOD,CR⫺1
N content of CS, mgN mgCOD,CS⫺1
N content of CI, mgN mgCOD,CI⫺1
N content of biomass, mgN mgCOD,BM⫺1
P content of CR, mgP mgCOD,CR⫺1
P content of CS, mgP mgCOD,CS⫺1
P content of CI, mgP mgCOD,CI⫺1
P content of biomass, mgP mgCOD,BM⫺1

of material c. The resulting equations for the stoichiometric coefficients ␷j,N, and ␷j,P are given in Tables 5 and
6, respectively.

⫺1
1 ⫺ fBM,CR ⫺fBM,CI
fBM,CI
Autotrophic organisms 2—Nitrobacter
fBM,CR1⫺(3.43 ⫺ YANs)/
YANs

1 ⫺ fBM,CR ⫺fBM,CI⫺1
fBM,CI
Autotrophic organisms 1—Nitrosomonas

⫺1/YH

⫺1/YH

fBM,CR

1

1

1
1
2

Hydrolysis
Aerobic growth of heterotrophs onreadily biodegradable COD
Nitrate-based growth of heterotrophs
on readily biodegradable COD
Nitrite-based growth of heterotrophs
on readily biodegradable COD
Lysis of heterotrophs

1 ⫺ 1/YH

1 ⫺ fHyd,CR
⫺1/YH

gCOD
Heterotrophic organisms⫺1fHyd,CI
gO2

Value

Values: Henze et al. (2000) except * Nowak (1996).

1
⫺1

7
XANb

␷1,N
␷2,N

fBM,CI

␷4,N

gN

10
NO3N
9
NO2N
8
NH4N

fHyd,CI
fBM,CR

6
XANs
5
XH
4
CI
3
CS
2
CR
1
O2
N:
Components:
R

Table 3. Stoichiometric matrix of reactions in CW2D.†

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

11
N2N

12
IP

VADOSE ZONE J., VOL. 4, NOVEMBER 2005

where as before i ⫽ 1,..., N, j ⫽ 1,..., R, and where ri is
the reaction term for component i (M L⫺3 T⫺1), ␪ is thevolumetric water content (L3 L⫺3), ␷j,i is the stoichiometric factor for component i and process j (M M⫺1) (Table3), and rcj is the zero-order reaction rate for process jin the aqueous phase (M L⫺3 T⫺1) as defined in Eq. [4]through [12].
Exchange of Oxygen
The exchange of O2 from the gas phase into the aqueous phase is described with a model commonly used inwastewater treatment (e.g., Gujer and Boller, 1990):
dcO2⫽ akaer,O2 (cO2,sat ⫺ cO2)
dt

[15]

where a is the air content (L3 L⫺3), kaer,O2 is the O2
Table 5. Stoichiometric coefficients for NHⴙ4 ⫺N.†
␷1,N ⫽ iN,CS ⫺ (1 ⫺ fHyd,CI)iN,CR ⫺ fHyd,CIiN,CI
␷2,N ⫽ 1/YHiN,CR ⫺ iN,BM
␷3,N ⫽ 1/YHiN,CR ⫺ iN,BM
␷4,N ⫽ 1/YHiN,CR ⫺ iN,BM
␷5,N ⫽ iN,BM ⫺ (1 ⫺ fBM,CR ⫺ fBM,CI)iN,CS ⫺ fBM,CRiN,CR ⫺ fBM,CIiN,CI
␷6,N ⫽ ⫺1/YANs ⫺ iN,BM
␷7,N ⫽ iN,BM ⫺ (1 ⫺ fBM,CR ⫺ fBM,CI)iN,CS ⫺ fBM,CRiN,CR ⫺ fBM,CIiN,CI
␷8,N ⫽ ⫺iN,BM
␷9,N ⫽ iN,BM ⫺ (1 ⫺ fBM,CR ⫺ fBM,CI)iN,CS ⫺ fBM,CRiN,CR ⫺ fBM,CIiN,CI[1]† See Table 4 for definitions of composition and stoichiometric parameters.

929

www.vadosezonejournal.org

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Table 6. Stoichiometric coefficients for inorganic phosphorus.†
␷1,P ⫽ iP,CS ⫺ (1 ⫺ fHyd,CI)iP,CR ⫺ fHyd,CIiP,CI
␷2,P ⫽ 1/YHiP,CR ⫺ iP,BM
␷3,P ⫽ 1/YHiP,CR ⫺ iP,BM
␷4,P ⫽ 1/YHiP,CR ⫺ iP,BM
␷5,P ⫽ iP,BM ⫺ (1 ⫺ fBM,CR ⫺ fBM,CI)iP,CS ⫺ fBM,CRiP,CR ⫺ fBM,CIiP,CI
␷6,P ⫽ ⫺iP,BM
␷7,P ⫽ iP,BM ⫺ (1 ⫺ fBM,CR ⫺ fBM,CI)iP,CS ⫺ fBM,CRiP,CR ⫺ fBM,CIiP,CI
␷8,P ⫽ ⫺iP,BM
␷9,P ⫽ iP,BM ⫺ (1 ⫺ fBM,CR ⫺ fBM,CI)iP,CS ⫺ fBM,CRiP,CR ⫺ fBM,CIiP,CI[2]† See Table 4 for definitions of composition and stoichiometric parameters.

reaeration rate (T⫺1), and cO2,sat is the temperature dependent saturation concentration of DO (M L⫺3). This
term must be added to the reaction rate constant for O2
since this is the only way to increase the concentration ofDO (with the exception of water flow). Note that theexchange of O2 in Eq. [15] is proportional to the levelof disequilibrium in the liquid phase and that the effectof the gas phase is included only by multiplying the ratewith the air content.

Implementation of CW2D into HYDRUS-2D
The original HYDRUS-2D code considers consecutive first-order decay and degradation, involving thesequential degradation of one into another solute. Thisapproach allows one to solve the transport equations
for individual species of the decay chain sequentiallysince the second solute depends on the first one (andso on), but not vice versa. However, since various components in the new reactive CW2D module depend ina relatively complex manner on each other (e.g., degradation of one species produces multiple other species),
the solute transport equations for all components mustbe evaluated simultaneously at each iteration step of
the numerical solution. The iteration procedure wasmodified such that the iteration included all components
simultaneously, rather than being performed at the levelof each single component. If convergence was not reached
for any single component, all components were recalculated during the next iteration. A self-adjusting time
stepping scheme was required to provide a relatively
smooth and efficient numerical solution. We note that
Galerkin finite element method was still used to numerically solve both the variably saturated flow and solutetransport equations.

APPLICATIONS TO PILOT-SCALECONSTRUCTED WETLANDS
Pilot-Scale Subsurface Vertical Flow
Constructed Wetland
System Description
Figure 1 shows a schematic of a pilot-scale subsurface
vertical flow constructed wetland (PSCW). The PSCW
is based on Austrian standards for vertical flow CWs
for wastewater treatment (ÖNORM B 2505, 1997). The
surface area of the wetland was 1 m2, and the height ofthe main layer of the filter bed 60 cm. The main layer
consisted of sand (gravel size 0.06–4 mm or 1–4 mm in
two PSCWs, respectively). The sandy material used forthe main layer had a porosity of 0.30 and a saturated

Fig. 1. Schematic of the pilot-scale subsurface constructed wetlands(values are in centimeters).

hydraulic conductivity, Ks, of 117 cm h⫺1. An intermediate layer of 10-cm thickness with a gravel size of 4 to
8 mm prevents fine particles from being washed out
into the drainage layer (15 cm thick; gravel 16–32 cm)
where the effluent was collected by means of tile drains(Langergraber, 2003).
The PSCW was loaded intermittently with wastewater. Following Austrian standards, the daily loading
rate was 40 L at application intervals of 6 h (i.e., a single
loading comprised 10 L).
Tracer experiments using NaCl were run for daily
hydraulic loading rates of 40 and 60 L. The electricalconductivity of the effluent, which is directly related tothe tracer concentration, was measured automatically.
The measured influent conductivity of the tracer was
30 mS cm⫺1, while the background conductivity during
the tracer experiment was 1.5 mS cm⫺1. Influent andeffluent concentrations of various N species and totalorganic carbon (TOC) were measured periodically. Table 7 shows the median values and confidence intervalsof the measured influent and effluent concentrations
for NH4–N, NO3–N, and TOC.
The width of the transport domain in the numericalsimulations was 1 m and its depth 0.6 m (only the main
layer was simulated), while the transport domain itself
was discretized into 11 columns and 40 rows. This results
in a two-dimensional finite element mesh consisting of
440 nodes and 780 triangular finite elements. An atmospheric boundary condition was assigned to the top of
the system, a constant pressure head boundary condition (constant head of ⫺2 cm) to the bottom of the
0.06- to 4-mm main layer and a free drainage boundary
condition to the bottom of the 1- to 4-mm main layer.
Water flow and single-solute transport were evaluatedusing the standard HYDRUS-2D software, while multiTable 7. Median values and confidence intervals (in parenthesis)of the measured influent and effluent concentrations of NHⴙ4 ,NO⫺
3 , and total organic carbon (TOC) (Langergraber, 2003).
Component
Influent
Effluent 0.06–4 mm substrate
Effluent 1–4 mm substrate

NH4–N

NO3–N

TOC

60.0 (11.2)
0.15 (0.03)
1.20 (0.31)

mg L⫺13.0 (0.8)
38.5 (6.8)
50.0 (8.7)

82.9 (3.5)
3.9 (0.2)
6.4 (0.3)

930

VADOSE ZONE J., VOL. 4, NOVEMBER 2005

Table 8. Literature and fitted soil hydraulic parameters of the van Genuchten–Mualem model for the 0.06- to 4-mm substrate.
Set

Parameter

Residual watercontent ␪r

␣

0.43
0.30
0.30
0.30
0.289
0.41

m⫺1
0.145
0.145
0.152
0.15960.126
0.145

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

m 3 m ⫺3
1
2
3
4
5

Literature
Measured Ks ⫹ porosity
Fitted 1
Fitted 2
Langergraber (2003)
1–4 mm main layer

Parameter

Saturated watercontent ␪s

0.045
0.045
0.08660.07240.056
0.060

component reactive transport was calculated using the
CW2D/HYDRUS-2D software.
Flow Simulations
Table 8 shows various parameter sets of the van Genuchten–Mualem soil hydraulic property model used inflow simulations (Fig. 2 and 3). The literature ParameterSet 1 was based on default parameters provided by
HYDRUS-2D for sand using pedotransfer functions
published by Carsel and Parrish (1988). The shape factors for sand were also used for Parameter Set 2 inTable 8, but now with the measured values of the saturated water content ␪s and the saturated hydraulic conductivity Ks. Parameter Set 3 in Table 8 was obtainedby calibrating the model against cumulative effluent
data, while Parameter Set 4 represents parameters optimized using the effluent rates. Parameter Set 5 are thosereported by Langergraber (2003), who in addition to theeffluent rates, also used water content measurements atthe 25-cm depth, while simultaneously optimizing thesaturated hydraulic conductivity Ks. Finally the parameters used for the coarser gravel (1–4 mm) are also given.
Figures 2 and 3 show the measured and simulated
cumulated effluents and effluent flow rates, respectively. Simulation results are presented for the parameter sets given in Table 8. The simulation using ParameterSet 3 showed a good match with the measured cumulated effluent (Fig. 2). However, peak effluent rates
were significantly overpredicted (Fig. 3). Parameter Set
4 described the effluent rates much better. These results

n
2.68
2.68
2.472
2.214
1.92
5

l

Saturated hydraulicconductivity Ks

0.5
0.5
0.696
0.682
0.5
0.5

m h⫺ 1
0.297
1.17
1.17
1.17
0.84
30

indicate that a good calibration of the flow model requires measurements of at least the saturated watercontent ␪s and the saturated hydraulic conductivity Ksof the main layer, as well as the effluent flow rate between two successive loadings. However, the best fit was
obtained (Fig. 4) with the parameter set of Langergraber(2003); this indicates that additional water content data
and more degrees of freedom in the optimization (i.e.,
one additional parameter was fitted) can result in abetter calibration of the flow model.
Simulation of Single-Component Transport
Results for a single component solute transport werepreviously presented by Langergraber (2003) and are
given here for completeness. Table 9 shows parametersthat were used for the general equilibrium solute transport model, and Table 10 gives parameters for the physical nonequilibrium transport model. The equilibriumtransport model assumes that all pore water is mobile,whereas the physical nonequilibrium transport model
considers a fraction of the pore water to be immobile(stagnant). Figure 5 compares simulated and measured
breakthrough curves. It can be seen that for both hydraulic loading rates the simulation results matched themeasured data relatively well.
Simulation of Multicomponent Reactive Transport
Table 11 shows measured influent and effluent concentrations of selected components, as well as the influ-

Fig. 2. Measured and simulated cumulated effluent for a single application of 10 L (four loadings per day; daily hydraulic load of 40 L).

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

931

Fig. 3. Measured and simulated effluent rates for a single application of 10 L (four loadings per day; daily hydraulic load of 40 L).

ent concentrations used in the CW2D/HYDRUS-2D
simulation, and the resulting simulated effluent concentrations. The simulated and measured effluent concentrations of NH4⫹ are lower for the 0.06- to 4-mm gravel
than for the coarser material (gravel size of 1–4 mm).
The higher flow rates in the main layer with the coarse
material reduced the contact time between water andthe bacteria. The effluent concentrations of NH4⫹–N andthe OM CS were significantly higher for the coarser
material. Effluent NO3⫺ concentrations were also higher
because of less readily degradable OM CR producedby hydrolysis that could be utilized for denitrificationin deeper layers.
Results of the reactive transport simulations showed
an overall good match with the measured data. Measured higher NH4⫹–N and NO3⫺–N concentrations in thesystem with the 1- to 4-mm substrate could be welldescribed by the model. However, further investigations
regarding parameters of the biokinetic model are stillrequired. The close match could be partly explainedby the calibration of the flow model, which was madepossible by having data that provided a good definitionof the flow processes in the system.
Below we present several figures showing time seriesof the simulation results for two loadings at three locations (5-, 15-, and 50-cm depths) in the 60-cm main layerof the PSCW system filled with gravel size 0.06 to 4 mm.
Figure 6 shows time series for the pressure head and
DO. Oxygen depletion at depths of 5 and 15 cm canbe observed to occur immediately after loading due todegradation. However, the O2 concentration at 50 cm
stayed constant at a level close to the saturation concentration during the entire loading interval. After ⬍2h, the O2 concentration at depths of 5 and 15 cm reachedsaturation again.
Time series and vertical concentration profiles at theend of the simulation (12 h) are shown in Fig. 7 forheterotrophic organisms XH. During two loading intervals the concentrations of the microorganisms remained
constant, thus indicating quasi steady state conditions
when growth and lysis processes are in balance (Fig. 7,left). Almost all microorganisms were present in the top
20 cm of the filter (Fig. 7, right), where most transformations and degradation processes took place. Autotrophic bacteria were distributed in a very similar manner.
The maximum concentration of autotrophic bacteria
was simulated to be lower than the maximum concentration of heterotrophic bacteria (data not shown).

Fig. 4. Measured and simulated effluent rates (left) and cumulated effluent (right) for a single feeding of 10 L (four loadings per day, daily
hydraulic load of 40 L).

932

VADOSE ZONE J., VOL. 4, NOVEMBER 2005

Table 9. Parameters of the general equilibrium solute transportmodel (Langergraber, 2003).
Parameter
Unit
Value

DiffusioncoefficientLongitudinal
dispersivity

m 2 h⫺1
0.00050.0125

Transverse
dispersivitymTable 10. Parameters for the physical nonequilibrium transportmodel (Langergraber, 2003).
Fraction of sorption Immobile
Coefficient for solutesites in contact
water
exchange between mobile
Parameter with mobile watercontent
and immobile water
Value

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

m 3 m ⫺3
0.05

0.01

Figure 8 shows time series for NH4⫹–N and NO3⫺–N,
respectively. Effluent concentrations for NH4⫹–N andNO3⫺–N were about 0.01 and 41 mg L⫺1, respectively.
Vertical concentration profiles at simulation times of
1 min (end of loading), 6 min, 30 min, and 2 h (onlyfor NO3⫺) are shown in Fig. 9 for NH4⫹–N and NO3⫺–N.
Notice that the concentration peaks for NH4⫹ at 5 and15 cm coincide with minima in the NO3⫺ simulations.
By contrast, concentrations of NH4⫹ and NO3⫺ remained
very stable at 50 cm depth.

Two-Stage Subsurface Vertical Flow
Constructed Wetland
System Description
The Small-Scale Plot (SSP) system described by Grosse
et al. (1999) is a PSCW for the treatment of heavily
polluted surface water. The total surface area of the
SSP was 2 m2, divided into downflow and upflow chambers (each having a surface area of 1 m2). The inlet was
situated on top of the downflow chamber. The effluent
was collected in the upflow chamber by means of perforated pipes. A separating wall was built between thedownflow and upflow chambers. Figure 10 shows a schematic of the SSP (Perfler et al., 1999). Heights of thedownflow and upflow chambers were 85 and 60 cm,
respectively. The drainage layer at the bottom of bothchambers had a height of 15 cm. This layer was connected to the two chambers and contained two tile
drains. The main 55- and 45-cm-thick layers in the downflow and upflow chambers, respectively, were filled withsand (gravel size 0.06–4 mm). The top layer of the downflow chamber had a thickness of 15 cm (gravel size 4–8

0.95

h⫺1
2

Table 11. Measured and simulated influent and effluent concentrations.
Component

NH4–N

NO3–N

TOC

COD

mg L⫺1
Measured (median)Simulation (input)
0.06-4 mm substrate
Measured (median)Simulation (median)
Measured (median)Simulation (median)

60.0
60.0

Influent
3.0
3.0
Effluent

CR
mgCOD

82.9
–

0.15
38.5
3.9
0.01
41.1
–
1–4 mm substrate
1.20
50.0
6.4
0.16
63.0
–

–
300

CS

L ⫺1

–
150

–
130

–
23.8

–
0.6

–
0.4

–
30.1

–
0.3

–
5.7

mm). Algae growth in the main layer was prevented bykeeping the water table during application beneath thesoil surface.
The unstructured two-dimensional finite element meshused for the SSP simulations consisted of 1135 nodes
and 2057 elements. The inlet distribution pipe at the topof the downflow chamber was modeled as a time dependent atmospheric boundary condition. Outlet (effluent)drainage pipes were modeled as circular openings in
the transport domain, and subjected to seepage faceboundary conditions. The soil hydraulic parameters ofthe van Genuchten–Mualem model used for the simulations are presented in Table 12. Values for the sandwere the same as those in Table 8 (Parameter Sets 1
and 2), while the parameters for the top and drainage
layer were again the same as those obtained by Langergraber (2001).
Since simulations of single-component reactive transport using the physical equilibrium and nonequilibrium

Fig. 5. Measured and simulated effluent electrical conductivities for different daily hydraulic loads.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

Fig. 6. Pressure head h (left) and dissolved oxygen (DO) (right).

Fig. 7. Simulated time series (left) and vertical profile (right) of heterotrophic organisms XH.

Fig. 8. Simulated time series for NH4ⴙ–N (left) and NO⫺
3 –N (right).

933

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

934

VADOSE ZONE J., VOL. 4, NOVEMBER 2005

Fig. 9. Vertical profiles of NHⴙ4 –N (left) and NO⫺
3 –N (right) at t ⫽ 1 min (end of loading), 6 min, 30 min, and 2 h.

(mobile–immobile water) transport models were presented by Langergraber (2003), we will focus below on
multicomponent reactive transport. Langergraber (2003)
showed that the presence of stagnant (immobile) regions of the pore water produced a much better matchwith the measured concentration data, compared withconsidering only mobile pore water.
Simulation of Multicomponent Reactive Transport
Table 13 shows influent concentrations for the multicomponent reactive transport simulations. Since thissystem was used to treat polluted surface water, theconcentrations were much lower compared with wastewater. The hydraulic loading rate was 50 L every 2 h(i.e., 600 L d⫺1).
Figure 11 shows the distribution of heterotrophic organisms XH. The maximum XH concentration was foundto be in the upper part of the main layer of the downflow
chamber. Bacteria concentrations decreased with depth.
Concentrations were low in the upflow chamber becauseof a lack of substrate available.
Figure 12 shows the distribution of Nitrosomonas
XANs. The distribution of the second autotrophic organism, Nitrobacter XANb, was found to be very similar(not further shown here). Again, most autotrophic bacteria were located in the upper part of the downflow
chamber. The maximum concentration of autotrophic
bacteria was lower than the maximum concentration of
the heterotrophic bacteria. The concentration decreased
faster with depth for the autotrophic bacteria than forthe heterotrophic bacteria. This is a direct consequenceof nitrification being a strictly aerobic process, with O2
being available only in the upper parts of the downflow
chamber (see also Fig. 13).
Figure 13 shows time series of DO for two consecutive
loadings to the main layer in the downflow chamber (0cm) at two different depths (5 and 15 cm, located in theunsaturated and saturated zones, respectively). Afterloading, the DO concentration decreased at the surface
and at the 5-cm depth. At the 15-cm depth, however,the DO concentration increased due to advective DOtransport with the infiltrating water. Still, O2 decreased
quickly (within 15 min after loading) due to the consumption of O2. No O2 was found at the 15-cm depth(i.e., in the saturated zone) during the remainder of thesimulation. In the unsaturated zone (5-cm depth) the
DO concentration increased again due to re-aeration,and reaching O2 saturation after about 1 h after loading.
Figure 14 shows time series of NH4⫹–N and NO3⫺–Nat the surface of the main layer (0 cm) and at threedifferent depths (5 cm, located in the unsaturated zone,
and 15 and 30 cm, located in the saturated zone). After⫺loading, NH⫹
4 initially increased whereas NO3 decreased.
⫹ is nitrified, resulting in higher NO⫺ concenLater, NH4
3
trations. At the 15-cm depth, changes in the concentration occurred mainly due to advection since O2 was
already a limiting factor (Fig. 13). Notice, again, that
only small changes in the concentrations occurred at the
30-cm depth. Figure 15 shows simulated and measured
concentrations of NH4⫹–N and NO3⫺–N in a cross sectionof the downflow chamber. The simulation results matched
the measured data well.

DISCUSSION

Fig. 10. Schematic of the two-stage vertical flow constructed wetlands(distances are in centimeters).

Our simulation results showed a good match withthe measured data for both single-component (tracer
experiments, e.g., Fig. 5) and multicomponent reactive
transport (e.g., Table 11 and Fig. 14) when the hydraulicbehavior of the system was accurately described (e.g.,
effluent flow rates in Fig. 2 and 4). The importance of
the water flow model was also evident from a sensitivity
analysis that showed that results of the reactive transport simulations were very sensitive to model parameters describing the unsaturated properties of the sub-935

www.vadosezonejournal.org

Table 12. van Genuchten–Mualem soil hydraulic parameters for the two-stage constructed wetlands system.
Residual watercontent ␪rParameter

␣

0.41
0.30
0.15

m⫺1
0.145
0.145
0.145

m 3 m⫺3

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

Top layer
Main layer
Drainage layer†

0.045
0.045
0.056

Parameter

Saturated watercontent ␪slSaturated hydraulicconductivity Ks

5.0
2.68
1.92

0.5
0.5
0.5

m h⫺1
600
1.17
600

n

† Drainage layer was always saturated.

Table 13. Influent concentrations for the multicomponent reactive transport simulations.
Parameter

O2

CR

CS

CI

NH4N

NO2N

NO3N

N2N

IP

Value

1

15

1

1

6.5

0.1

2

10

0.4

strate. The good match of the experimental data withthe reactive transport simulations could be obtainedusing literature values for the CW2D model parameters(Langergraber, 2001).
The CW2D module considers the main constituentsof domestic wastewater. CW2D hence is only applicable
to CWs that treat domestic wastewater or water withproperties that can be described with processes andparameters currently incorporated into CW2D. If othercomponents and/or waters are to be modeled, an extension to the present model is needed.
Further limitations of CW2D result from assumptionsof the basic model, the Activated Sludge Model, used
here. These include consideration of (i) a constant valueof pH, (ii) constant coefficients in the rate equations,(iii) constant stoichiometric factors, and (iv) a limited
temperature range from 10 to 25⬚C (Henze et al., 2000).
Although simulations showed a good match with themeasured data, much research remains so as to produce
a reliable tool for CW design. The main research needs
regarding modeling of CW processes as listed by Langergraber (2003) are:Detailed Hydraulic Investigations of Full-Scale CWs.
A good match between simulation results and measured
data can be obtained when sufficient data are available
to describe the hydraulic behavior of the system, andthus when the flow model can be successfully calibrated.
For full-scale systems, the match between the modeland data was thus far relatively poor, most probably
due to a lack of data to accurately describe the hydraulicbehavior of the system (Langergraber, 2001) and dueto the inherent heterogeneity of large-scale systems. The
heterogeneity of the substrate and uneven distributionof influent water on the surface are the main reasonsfor needing detailed hydraulic investigations of fullscale systems.
Numerical Simulations of Outdoor CW Systems. Simulations of outdoor (full-scale) CW systems are alsoneeded to validate the temperature dependencies of thevarious processes considered in CW2D.
Investigation of Plant Uptake Models. Models describing nutrient (N and P) uptake by plants have to beconsidered to describe the contribution of plants to theoverall treatment process. Langergraber (2004) showed
the important influence of plant uptake on effluent concentrations and removal efficiency using plant uptake
models that are currently available in HYDRUS-2D(Šimůnek et al., 1999). These models describe nutrient
uptake associated with water uptake (i.e., active uptake). For high loaded systems (e.g., CWs that treat
wastewater) the model gives good results, whereas forwaters with lower nutrient contents the simulation results indicate that potential nutrient uptake will be substantially overestimated when using literature values tocalculate the potential uptake rate, compared with simulations using models that couple nutrient uptake to water
uptake.
Need for a Module Describing Substrate Clogging.
One limitation of CW2D is that up to now only dissolved
solutes are being considered. Wastewater also contains
suspended solids that may cause pore clogging. Clogging
is considered to be the biggest operational problem affecting vertical flow in CWs. It causes a reduction in
the infiltration capacity of the substrate surface, andtherefore often a rapid decline in the treatment performance of the system (e.g., Platzer and Mauch, 1997;
Winter and Goetz, 2003; Langergraber et al., 2003). The
clogging module should be able to describe pore size

Fig. 11. Simulated steady state distribution of heterotrophic organisms XH.

Fig. 12. Calculated steady state distribution of Nitrosomonas XANs.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

936

VADOSE ZONE J., VOL. 4, NOVEMBER 2005

Fig. 15. Simulated and measured concentrations of NHⴙ4 –N andNO⫺
3 –N in a cross section of the downflow chamber.
Fig. 13. Simulated time series of dissolved oxygen (DO).

reductions due to the settling of suspended solids anddue to bacteria growth. Both particle transport and bacteria growth models are currently available (e.g., Mackie
and Bai, 1993; Vandevivere et al., 1995; Kildsgaard andEngesgaard, 2001) and could be incorporated into
CW2D/HYDRUS-2D. The effects of pore size reductions on the hydraulic properties of the substrate may
also need to be taken into account (e.g., Taylor and Jaffé,1990; Magnico, 2000).
Improved Experimental Methods for Estimating
CW2D Model Parameters. There is also a need for
improving technologies to measure CW2D model parameters rather than always requiring calibration of thereactive transport part of the model. This is important
before the model can be used as a design tool. A numberof measurement techniques have been developed tocharacterize the parameters of Activated Sludge Models
as applied to activated sludge wastewater treatment
plants (Henze et al., 2000; Vanrolleghem et al., 1999).
However, no experimental techniques are currently available to measure CW2D model parameters. Detailed investigations on microbial biocoenosis in particular mustbe performed, with a need to adopt existing measurement
methods to soil and/or wastewater environments.

Fig. 14. Simulated time series of NHⴙ4 –N (left) and NO⫺
3 –N (right).

SUMMARY AND CONCLUSIONSWe presented a multicomponent reactive transportmodel, CW2D, which is an extension of the HYDRUS-2D
variably saturated water flow and solute transport program, to include biochemical transformation and degradation processes in subsurface flow CWs. Constructedwetlands represent a very complex mixture of water,
substrate, plants, litter, and a variety of microorganisms,
in which a large number of simultaneous processes takes
place. A complete description of such a complex system
involving mutually dependent physical, chemical, and biological reactions and chemical and biological compounds
is nearly impossible without the help of a numerical modelthat considers the majority of these interactions. The waterflow regime in subsurface vertical flow CWs is also highly
dynamic, adding to the complexity of the overall system
by requiring a transient variably saturated flow model. By
developing CW2D/HYDRUS-2D we have attempted todevelop a numerical tool that can capture at least the mainprocesses of such a complex system, and that can help
in analyzing and studying the various transient processes
operative in CWs.
In this paper we described the main processes andcomponents considered in CW2D and demonstrated
the performance of the model for one- and two-stageReproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

www.vadosezonejournal.org

subsurface vertical flow CWs. Model simulations forwater flow, tracer transport, and selected biochemicalcompounds were compared against experimental observations. We also discussed limitations of the model andneeds for further improvements of the model. The research needs are largely those listed by Langergraber(2003), but were updated to reflect results of our recent investigations.
CW2D/HYDRUS-2D was developed to model processes in subsurface flow CWs intended for wastewater
treatment. However, the model has already been successfully applied to CWs for the treatment of combined
sewer overflows (Dittmer et al., 2004). Possible futureapplications include (i) CWs for other types of wastewater, such as industrial and agricultural wastewaters
having compositions similar to municipal wastewater;(ii) infiltration of treated (waste)water into soil forgroundwater recharge or other applications; (iii) C andN cycling in soils; (iv) processes in natural wetlands;
and (v) processes in riparian areas.
ACKNOWLEDGMENTS
The senior author was partly funded by the Austrian research project “Characterisation of Microbial Biocoenosis to
Optimise Removal Efficiency and Design of Subsurface Flow
Constructed Wetlands for Wastewater Treatment” (fundedby the Austrian Science Fund, FWF, Project P16212-B06) andpartly by the Agricultural Experiment Station of the University of California Riverside. The work of Dr. Šimůnek was
supported in part by SAHRA (Sustainability of Semi-Arid
Hydrology and Riparian Areas) under the STC Program of
the National Science Foundation, Agreement EAR-9876800
and the Terrestrial Sciences Program of the Army Research
Office (Terrestrial Processes and Landscape Dynamics andTerrestrial System Modeling and Model Integration).

REFERENCES
Abrams, R.H., K. Loague, and D.B. Kent. 1998. Development andtesting of a compartmentalized reaction network model for redox
zones in contaminated aquifers. Water Resour. Res. 34:1531–1541.
Brouwer, H., A. Klapwijk, and K.J. Keesman. 1998. Identification ofactivated sludge and wastewater characteristics using respirometricbatch-experiments. Water Res. 32:1240–1254.
Carsel, R.F., and R.S. Parrish. 1988. Developing joint probability
distributions of soil water retention characteristics. Water Resour.
Res. 24:755–769.
Dittmer, U., D. Meyer, and G. Langergraber. 2004. Simulation of a
subsurface vertical flow constructed wetland for CSO treatment.
p. 511–518. In A. Lienard and H. Burnett (ed.) Proceedings of the
9th IWA Specialized Group Conf. Wetland Systems for Water
Pollution Control, Avignon, France. Vol. 2. 26–30 Sept. 2004.
Essaid, H.I., B.A. Bekins, E.M. Godsy, E. Warren, M.J. Baedecker,
and I.M. Cozzarelli. 1995. Simulation of aerobic and anaerobic
biodegradation processes at a crude oil spill site. Water Resour.
Res. 31:3309–3327.
Fryar, A.E., and F.W. Schwartz. 1994. Modeling the removal of metals
from groundwater by a reactive barrier: Experimental results. Water Resour. Res. 30:3455–3469.
Grosse, W., F. Wissing, R. Perfler, Z. Wu, J. Chang, and Z. Lei. 1999.
Biotechnological approach to water quality improvement in tropicaland subtropical areas for reuse and rehabilitation of aquatic ecosystems. Final report, INCO-DC Project Contract ERBIC18CT960059.
INCO-DC, Cologne, Germany.
Gujer, W., and M. Boller. 1990. A mathematical model for rotating
biological contactors. Water Sci. Technol. 22:53–73.
Haberl, R., S. Grego, G. Langergraber, R.H. Kadlec, A.R. Cicalini,
S. Martins Dias, J.M. Novais, S. Aubert, A. Gerth, H. Thomas,

937

and A. Hebner. 2003. Constructed wetlands for the treatment oforganic pollutants. J. Soils Sediments 3:109–124.
Henze, M., W. Gujer, T. Mino, and M.C.M. van Loosdrecht. 2000.
Activated sludge models ASM1, ASM2, ASM2D and ASM3. IWA
Scientific and Technical Rep. 9. IWA Publ., London.
Kadlec, R.H. 2000. The inadequacy of first-order treatment kinetic
models. Ecol. Eng. 15:105–119.
Kadlec, R.H., and R.L. Knight. 1996. Treatment wetlands. CRC Press,
Boca Raton, FL.
Kadlec, R.H., R.L. Knight, J. Vymazal, H. Brix, P. Cooper, and R.
Haberl (ed.) 2000. Constructed wetlands for pollution control—
Processes, performance, design and operation. IWA Scientific andTechnical Rep. 8. IWA Publ., London.
Kildsgaard, J., and P. Engesgaard. 2001. Numerical analysis of biological clogging in two-dimensional sand box experiments. J. Contam.
Hydrol. 50:261–285.
Langergraber, G. 2001. Development of a simulation tool for subsurface flow constructed wetlands. Wiener Mitteilungen 169, Vienna,
Austria.
Langergraber, G. 2003. Simulation of subsurface flow constructedwetlands—Results and further research needs. Water Sci. Technol. 48(5):157–166.
Langergraber, G. 2004. The role of plant uptake on the removalof organic matter and nutrients in subsurface flow constructedwetlands—A simulation study. p. 491–498. In A. Lienard and H.
Burnett (ed.) Proceedings of the 9th IWA Specialized Group Conf.
Wetland Systems for Water Pollution Control, Avignon, France.
Vol. 2. 26–30 Sept. 2004.
Langergraber, G., and R. Haberl. 2001. Constructed wetlands forwater treatment. Minerva Biotecnologica 13(2):123–134.
Langergraber, G., R. Haberl, J. Laber, and A. Pressl. 2003. Evaluationof substrate clogging processes in vertical flow constructed wetlands. Water Sci. Technol. 48(5):25–34.
Lichtner, P.C. 1996. Continuum formulation of multicomponent-multiphase reactive transport. p. 1–81. In P.C. Lichtner et al. (ed.)
Reactive transport in porous media. Reviews in Mineralogy. Vol.
34. Mineralogical Society of America, Washington, DC.
Mackie, R.I., and R. Bai. 1993. The role of particle size distribution
on the performance and modelling of filtration. Water Sci. Technol. 27(10):19–34.
Magnico, P. 2000. Impact of dynamic processes on the coupling between fluid transport and precipitate deposition. Chem. Eng. Sci.
55:4323–4338.
Nowak, O., K. Svardal, and P. Schweighofer. 1995. The dynamic
behaviour of nitrifying activated sludge systems influenced by inhibiting wastewater compounds. Water Sci. Technol. 31(2):115–124.
Nowak, O. 1996. Nitrification in activated sludge processes with significant influence of industrial wastewater. (In German.) Wiener Mitteilungen 135, Vienna, Austria.
ÖNORM B 2505. 1997. Subsurface-flow constructed wetlands—
Application, dimensioning, installation and operation. (In German.)
Vornorm, Österreichisches Normungsinstitut, Vienna, Austria.
Parton, W.J., J.W.B. Stewart, and C.V. Cole. 1988. Dynamics of C,
N, P, and S in grassland soils: A model. Biogeochemistry 5:109–131.
Perfler, R., J. Laber, G. Langergraber, and R. Haberl. 1999. Constructed wetlands for rehabilitation and reuse of surface waters intropical and subtropical areas—First results from small-scale plots
using vertical flow beds. Water Sci. Technol. 40(3):155–162.
Platzer, Ch., and K. Mauch. 1997. Soil clogging in vertical flow reed
beds– mechanisms, parameters, consequences and... solutions? Water Sci. Technol. 35(5):175–181.
Rittmann, B.E., and J.M. VanBriesen. 1996. Microbiological processes
in reactive modeling. p. 311–334. In P.C. Lichtner et al. (ed.) Reactive transport in porous media. Reviews in Mineralogy. Vol. 34.
Mineralogical Society of America, Washington, DC.
Rockhold, M.L., R.R. Yarwood, M.R. Niemet, P.J. Bottomley, andJ.S. Selker. 2002. Considerations for modeling bacterial-inducedchanges in hydraulic properties of variably saturated porous media.
Adv. Water Resour. 25:477–495.
Scanlon, B.R., M. Christman, R. Reedy, I. Porro, J. Šimůnek, and
G. Flerchinger. 2002. Intercode comparisons for simulating waterbalance of surficial sediments in semiarid regions. Water Resour.
Res. 38(12):1323. doi:10.1029/2001WR001233.

Reproduced from Vadose Zone Journal. Published by Soil Science Society of America. All copyrights reserved.

938

VADOSE ZONE J., VOL. 4, NOVEMBER 2005

Schwager, A., and M. Boller. 1997. Transport phenomena in intermittent filters. Water Sci. Technol. 35(6):13–20.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 1999. The
HYDRUS-2D software package for simulating two-dimensional
movement of water, heat, and multiple solutes in variably saturated
media. Version 2.0. IGWMC-TPS-53. International Ground Water
Modeling Center, Colorado School of Mines, Golden, CO.
Šimůnek, J., and A.J. Valocchi. 2002. Geochemical transport. p. 1511–
1536. In J.H. Dane and G.C. Topp (ed.) Methods of soil analysis.
Part 4. SSSA Book Ser. 5. SSSA, Madison, WI.
Tanner, C.C., J. D’Eugenio, G.B. McBride, J.P.S. Sukias, and K.
Thompson. 1999. Effect of water level fluctuation on nitrogen removal from constructed wetland mesocosms. Ecol. Eng. 12:67–92.
Taylor, S.W., and P.R. Jaffé. 1990. Biofilm growth and the relatedchanges in the physical properties of a porous medium 3. Dispersivity and model verification. Water Resour. Res. 26:2171–2180.
Tebes-Stevens, C., A.J. Valocchi, J.M. VanBriesen, and B.E. Rittmann. 1998. Multicomponent transport with coupled geochemical
and microbiological reactions: Model description and example simulations. J. Hydrol. (Amsterdam) 209:8–26.
VanBriesen, J. 1998. Modeling coupled biogeochemical processes inmixed waste systems. Ph.D. diss. Civil Engineering, Northwestern
University, Evanston, IL.
Vandevivere, P., P. Baveye, D. Sanchez de Lozada, and P. DeLeo.
1995. Microbiological clogging of saturated soils and aquifer materials: Evaluation of mathematical models. Water Resour. Res. 31:
2173–2180.
van Genuchten, M.Th. 1980. A closed-form equation for predicting

the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am.
J. 44:892–898.
van Genuchten, M.Th., and J. Šimůnek. 2004. Integrated modelingof vadose zone flow and transport processes. p. 37–69. In R.A.
Feddes et al. (ed.) Proc. Unsaturated Zone Modelling: Progress,
Challenges and Applications, Wageningen, The Netherlands. 3–5
Oct. 2004,
van Genuchten, M.Th., and P.J. Wierenga. 1976. Mass transfer studies
in sorbing porous media: I. Analytical solutions. Soil Sci. Soc. Am.
J. 40:473–481.
Vanrolleghem, P.A., H. Spanjers, B. Petersen, P. Ginestet, and I.
Takács. 1999. Estimating (combinations of) Activated Sludge
Model No. 1 parameters and components by respirometry. Water
Sci. Technol. 39(1):195–214.
Viswanathan, H.S., B.A. Robinsons, A.J. Valocchi, and I.R. Triay.
1998. A reactive transport model of neptunium migration frompotential repository at Yucca Mountain. J. Hydrol. (Amsterdam)
209:251–280.
Walter, A.L., E.O. Frind, D.W. Blowes, C.J. Ptacek, and J.W. Molson.
1994. Modeling of multicomponent reactive transport in groundwater 2. Metal mobility in aquifers impacted by acidic mine tailings
discharge. Water Resour. Res. 30:3149–3158.
Werner, T.M., and R.H. Kadlec. 2000. Wetland residence time modelling. Ecol. Eng. 15:77–90.
Winter, K.-J., and D. Goetz. 2003. The impact of sewage compositionon the soil clogging phenomena of vertical flow constructed wetlands. Water Sci. Technol. 48(5):9–14.
Wynn, M.T., and S.K. Liehr. 2001. Development of a constructedsubsurface-flow wetland simulation model. Ecol. Eng. 16:519–536.