M

Modeling of Processes in Subsurface Flow
Constructed Wetlands: A Review

S

S

:V

Z

Günter Langergraber*Constructed wetlands (CWs) provide a natural way for simple, inexpensive, and robust wastewater treatment. Detailed
understanding of CW func oning is diﬃcult, because a large number of physical, chemical, and biological processes
occur in parallel and influence each other. For this reason, CWs have long been seen as “black boxes” where wastewater
enters and treated water leaves the system. Numerical models describing the biochemical transforma on and degrada on processes in CWs are promising tools to be er understand CW func oning. The first part of this paper reviews
published mechanis c models for CWs. Horizontal flow systems can be simulated when only water flow saturated condi ons are considered; six models have been reviewed whereby a series or network of completely s rred tank reactors
is most frequently used to describe the hydraulics. For modeling ver cal flow CWs with intermi ent loading, transientvariably saturated flow models are required. Due to the intermi ent loading, these systems are highly dynamic, adding
to the complexity of the overall system. Five models of diﬀerent complexity have been reviewed; three use the Richards
equa on to describe variably saturated flow, whereas the two others use simplified approaches. In the second part ofthe paper, the mul component reac ve transport module CW2D is demonstrated. Simula on results for CWs trea ng
domes c wastewater, combined sewer overflow, and surface water are presented. In general, a good match between
simula on results and measured data could be achieved if the hydraulic behavior of the system could be described well.
Based on the experience from these examples, the need for further model development is determined.
A: AO7, acid orange 7; ASM, Activated Sludge Model; BOD, biochemical oxygen demand; COD, chemical oxygen demand;
CSO, combined sewer overﬂow; CSTR, completely stirred tank reactor; CW, constructed wetland; HF, horizontal ﬂow; HRT, hydraulicresidence time; PSCW, pilot-scale constructed wetland; RCB, RetrasoCodeBright; RTD, residence time distribution; SF, surface ﬂow; SIR,
substrate-induced respiration; SSF, subsurface ﬂow; VF, vertical ﬂow.

Dexperience constructing andoperating CWs, these systems are still often consideredﬁgurative black boxes in which water is treated. Until now, CW
design has been mainly based on rule of thumb approaches using
speciﬁc surface area requirements (e.g., Brix and Johansen, 2004;ÖNORM B 2505, 2005; DWA-A 262, 2006) or simple ﬁrstorder decay models (e.g., Kadlec and Knight, 1996; Rousseau etal., 2004). However, Kadlec (2000) documents that ﬁrst-ordermodels are inadequate for the design of treatment wetlands.
Constructed wetlands are engineered to optimize the treatment conditions found in natural wetlands by using a complex
mixture of water, substrate, plants, litter (fallen plant material),
and a variety of microorganisms (especially bacteria). The large
number of physical, chemical, and biological processes operating
simultaneously and inﬂuencing each other makes the CW systemdiﬃcult to understand. Numerical models have drawn increased
attention in recent years because they can provide insight into

Ins tute of Sanitary Engineering and Water Pollu on Control, Univ. of
Natural Resources and Applied Life Sciences, Vienna (BOKU), Muthgasse
18, A-1190 Vienna, Austria. Received 22 Mar. 2007. *Corresponding
author (guenter.langergraber@boku.ac.at).
Vadose Zone J. 7:830–842
doi:10.2136/vzj2007.0054© Soil Science Society of America
677 S. Segoe Rd. Madison, WI 53711 USA.
All rights reserved. No part of this periodical may be reproduced ortransmi ed in any form or by any means, electronic or mechanical,
including photocopying, recording, or any informa on storage andretrieval system, without permission in wri ng from the publisher.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008830

the black box and therefore improve our understanding of thesystem. Once reliable numerical models exist, they can be used
to evaluate and improve the existing design criteria.
Wetland treatment systems eﬀectively treat organic matterand pathogens. Water plants (macrophytes) have several relevant
properties, the most important being their physical eﬀects (Brix,
1997, Kadlec et al., 2000). The amount of nutrients removedby plant harvesting is generally insigniﬁcant compared to yearlyloadings with wastewater. If plants are not harvested, the nutrients will be returned to the water during plant decomposition(Brix, 1997; Tanner, 2001). Applications of CWs include treatment of domestic, agricultural, and industrial wastewater; storm
water; and landﬁll leachate. In general, the use of CWs provides
a relatively simple, inexpensive, and robust solution for treatment. As natural treatment systems, CWs require a larger speciﬁc
surface area compared to technical solutions such as activatedsludge. However, CWs usually have lower operation and maintenance expenses and other additional beneﬁts, including tolerance
against ﬂuctuations of ﬂow and pollution load, ease of water reuse
and recycling, provision of habitat for many wetland organisms,
and a more aesthetic appearance than technical treatment options(Kadlec et al., 2000; Haberl et al., 2003).
Constructed wetlands can be subdivided into two main types:
surface ﬂow (SF) and subsurface ﬂow (SSF) CWs (Kadlec andKnight, 1996; Haberl et al., 2003). Surface ﬂow, or free water surface CWs, are densely vegetated and have typical water depths ofless than 0.4 m. Subsurface ﬂow CWs, which have no visible freewater, are subdivided into horizontal ﬂow (HF) and vertical ﬂow(VF) systems depending on the direction of water ﬂow throughthe porous medium (sand or gravel).
To prevent clogging of the porous
ﬁlter material, the use of traditional
SSF CWs is limited to mechanically
pretreated wastewater, which contains
less particulate content. Compared toSF systems, the contact area of waterwith bacteria and substrate is muchlarger, decreasing the area requirement
for SSF CWs (Vymazal et al., 1998).
In an HF CW (Fig. 1), water F . 1. Longitudinal cross-sec on of a horizontal flow constructed wetland (Langergraber andenters the inlet and ﬂows slowly under Haberl, 2004).
the surface through the porous media
until it reaches the outlet zone, whereregression analysis and artiﬁcial neural networks when used to
it is collected and discharged. Oxygen for aerobic processes ispredict biochemical oxygen demand (BOD) eﬄuent concentraobtained mainly via diﬀusion from the atmosphere. The amounttions of HF CWs. Stein et al. (2006) estimated the temperatureof oxygen transported from the roots into zones under the water
dependence of the model parameters of the k–C* ﬁrst-ordertable is still under discussion; however, it is too little to facilitate
degradation model proposed by Kadlec and Knight (1996) for
aerobic processes. Therefore, anoxic and anaerobic processes playchemical oxygen demand (COD) removal in a batch-fed CW
the most important role in HF CWs. Organic matter is decomsystem using four plant species. Rousseau et al. (2004) composed both aerobically and anaerobically, resulting in eﬃcientpared diﬀerent design models for HF CWs and demonstrated
removal. Insuﬃcient oxygen supplies results in incomplete nitrithe output variability of these models. The k–C* ﬁrst-order degﬁcation (Langergraber and Haberl, 2001).
radation model has been shown to be the most suitable for design
Vertical ﬂow CWs with intermittent loading are widely used
purposes. However, the authors claimed that this suitability istoday because they eﬃciently remove ammonia nitrogen. Figure
only true when model parameters from systems operating under2 shows the typical setup of a VF CW. Water is loaded intermitsimilar conditions (i.e., climatic conditions, wastewater compositently, and the large amount of water from a single loading causes
tion, porous ﬁlter material, and plant species) are available.
ﬂooding of the surface. The water inﬁltrates into the substrate,
This review article considers only mechanistic models. The
then gradually drains down vertically and is collected by a drainﬁrst section presents the available models in three diﬀerent catage network at the base. Air re-enters the system until the next
egories: (i) models describing hydraulic behavior and single-soluteloading, and great oxygen transfer rates into the system are possible.
transport only; (ii) models describing reactive transport in satuVertical ﬂow CWs with intermittent loading are therefore suitable
rated conditions, such as in HF CWs; and (iii) reactive transport
when nitriﬁcation and other strictly aerobic processes are required
models for variably saturated conditions, the most complex(Langergraber and Haberl, 2001, Kayser and Kunst, 2005). In themodels which are used to model both VF and HF CWs. The
last decade, VF CWs for treating raw wastewater have also been
second part of the article demonstrates the multicomponent reacintroduced and successfully applied (e.g., Molle et al., 2005).
tive transport module CW2D (Langergraber, 2001; Langergraber
Few numerical models are available to describe treatmentand Šimůnek, 2005) and provides examples of its application.
processes in SSF CWs. Much of the literature on models for
CWs (mostly HF CWs) refers to simple ﬁrst-order decay models
Review of MechanisƟc Models for Subsurface(e.g., Pastor et al., 2003; Stein et al., 2006; Tomenko et al., 2007)
Flow Constructed Wetlands
where only eﬄuent concentrations are predicted based on inﬂuent
Models Describing the Hydraulic Behavior andconcentrations. Pastor et al. (2003) used hybrid neural network
Single-Solute Transport in Constructed Wetlands
models to model the rate constant of the ﬁrst-order decay model
to predict eﬄuent concentrations of HF CWs with the goal of
Several authors have conducted tracer experiments to evaluoptimizing their design. Using a similar approach, Tomenko
ate the hydraulic behavior of SSF CWs. Tracer experiments are
et al. (2007) compared the accuracy and eﬃciency of multiple
also used to calibrate the hydraulic behavior of the system, including the more complex models described in the next
sections.
Werner and Kadlec (2000) modeled the nonideal
ﬂow of CWs with a network of an inﬁnite number ofsmall stirred tanks distributed along a set of main plugﬂow channels. Although the method has been developed mainly for SF wetlands, it can also be applied
to subsurface HF CWs. The basic concept of themodel assumes that the CW has an inﬁnite number ofmicrozones of diminished mixing along a set of main
channels. These zones are not excluded “dead zones,”
F . 2. Typical setup of a ver cal flow constructed wetland (Langergraber andbut they only exchange water with the main ﬂows
Haberl, 2004).
on a limited basis. In total, 47 tracer studies fromwww.vadosezonejournal.org · Vol. 7, No. 2, May 2008831

four diﬀerent sites were investigated. The model was capable ofproducing realistic wetland residence time distributions (RTDs),and reproduced the experimental RTDs from impulse tracer studies performed on both steady and variable ﬂow systems for the
more reliable and carefully performed experiments. The authors
concluded that their method was more ﬂexible than tanks-inseries models and could therefore better reproduce the diﬀerent
anomalies observed at the four sites.
Chazarenc et al. (2003) applied mathematical tools fromchemical engineering to model the hydraulic residence time(HRT) distribution of a HF CW in France. The nonideal ﬂow
wetland was modeled by dispersion plug ﬂow and stirred-tanksin-series models. Six experimental HRTs were monitored andmodeled. In dry periods, signiﬁcant dispersion occurred; in wet
periods, plug ﬂow behavior could be observed. Results from dispersion plug ﬂow and stirred-tanks-in-series models showed that
evapotranspiration played a major role in summer HRT simulations. Evapotranspiration seemed to reduce dispersion. It was
shown that the design of the system was favorable for developing
preferential ﬂow paths, which also gave rise to dead and stagnant
zones. Additionally, short cut ﬂows occurred due to surface ﬂow
during periods of heavy rain. The conceptual models developed
could ﬁt all response curves.
García et al. (2004) modeled tracer tests performed at pilot
HF CWs with diﬀerent ﬁlter materials and diﬀerent length/width
ratios of 1:1, 1.5:1, and 2:1. Nonideal ﬂow models (plug ﬂow
with dispersion and completely stirred tank reactors (CSTRs)
in series with delay) have been applied. The CSTRs-in-serieswith-delay model represented the data more accurately whencompared with the plug ﬂow with dispersion model. García etal. (2004) concluded that the former ﬁts the asymmetry of the
tracer response more accurately. The number of CSTRs required
ranged from one to three. The main conclusion of their study wasthat the construction of a HF CW with a higher length/width
ratio and ﬁner medium improves the hydraulic behavior of thesystem and reduces internal dispersion.
Małoszewski et al. (2006) used tracer experiments to determine hydraulic parameters in three parallel inhomogeneous gravel
beds at an HF CW in Poland. Instantaneously injected bromide
and tritium tracers were used to obtain RTDs. The multiﬂow
dispersion model assumed the existence of several ﬂowpaths
with diﬀerent hydraulic properties; it was developed using therespective parallel combination of analytical solutions from the
one-dimensional advection–dispersion equation. The model
was successfully used to ﬁt the experimental tracer breakthroughcurves, to identify the diﬀerent ﬂow components, and to derive
wastewater volumes, water-saturated porosity, mean wastewater
travel times, longitudinal dispersivities, and hydraulic conductivity from model parameters. The application of the model, which
describes possible diﬀusion of tracer into the zones with stagnant
water, has demonstrated that the calibration of such a model ispossible. However, the diﬀusion properties of the tracer obtained
from the model parameters were unrealistic, which suggested thatthe zones with stagnant water played an insigniﬁcant role in theinvestigated system.
Schwager and Boller (1997) simulated tracer experiments
and oxygen transport in intermittent sand ﬁlters using an older
version of HYDRUS-1D (Šimůnek et al., 1998) and MOFAT(Katyal et al., 1991), respectively. However, Schwager and Boller
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

(1997) did not consider reactive transport. The eﬀect of biomass
accumulation on solute breakthrough was assessed experimentallyby tracer studies. HYDRUS was used to extrapolate to diﬀerent
hydraulic conditions. MOFAT was used simultaneously to include
the eﬀects of water and air ﬂow. The simulations have been used
to estimate oxygen ﬂuxes by advection and diﬀusion into the top
10 cm of the ﬁlter. At a hydraulic loading rate of 120 L m−2 d−1,
the oxygen ﬂux was found to be 9 and 30 mg O2 m−2 min−1
caused by advection and diﬀusion, respectively. The calculated
results could be conﬁrmed by gas tracer experiments. The results
strengthened the assumption that oxygen mainly enters the poresin the intervals between intermittent loadings via diﬀusion.
Reac ve Transport Models for Saturated Condi ons

To simulate HF CWs, several authors have been coupling
reactive transport models and ideal reactor models such as aseries of CSTRs and/or plug ﬂow reactors, respectively. In thesereactors, either ﬁrst-order decay rates or Monod kinetics areapplied to model the degradation processes for organic matter
and/or nitrogen.
Chen et al. (1999) used the mixing cell method to model
the BOD elimination process. The HF bed is subdivided into
a number of same-sized cells that are assumed to be completely
mixed. In each cell, the degradation is described by a ﬁrst-order
rate constant. By using this approach, it is possible to derive
an analytical solution for the eﬄuent BOD concentration as afunction of the number of cells. The model is a simpliﬁcation of
the advection–dispersion partial diﬀerential equation. The dispersion coeﬃcient, the length, and the ﬂow rate of the HF bed
determined the number of cells required to get good results. The
applications of the model presented by Chen et al. (1999) indicate that the model is able to produce a better match to measured
data than a plug ﬂow model and can be also applied to transientconditions such as variable ﬂow rates.
A mechanistic, compartmental simulation model was developed by Wynn and Liehr (2001) to model and predict seasonal
trends in the removal eﬃciencies of HF CWs. The model consists of six linked submodels representing the carbon cycle, thenitrogen cycle, an oxygen balance, autotrophic bacteria growth,
heterotrophic bacteria growth, and a water budget. Bacterial
growth rates are modeled using Monod-type reaction kinetics.
The carbon cycle considers several state variables, including plantbiomass, dissolved and particulate organic carbon, and refractorycarbon. The wetland is assumed to act as either a single CSTR
or a series of CSTRs. Data from an existing CW were used tocalibrate the model. In general, the model predicts eﬄuent BOD,
organic nitrogen, ammonium, and nitrate concentrations well.
The model also reproduces seasonal trends well. Interactions
between the carbon, nitrogen, and oxygen cycles were evident
in model outputs. Because little was known about the root zone
aeration by wetland plants, oxygen predictions were fair. The
model was generally insensitive to changes in individual parameters, assumably because of the complexity of both the ecosystemand the model, as well as the numerous feedback mechanisms.
The model was most sensitive to changes in parameters that aﬀect
microbial growth and substrate use directly. Wynn and Liehr(2001) concluded that, with further evaluation and reﬁnement,
the model could be a useful design tool for HF CWs. However,no further work has been published since 2001.
832

Mashauri and Kayombo (2002) developed a coupled model fora waste stabilization pond and HF CW based on the growth rate
kinetics described by Monod kinetic equations. The transformation of organic carbon considered heterotrophic growth and die-oﬀ
of bacteria, and mineralization and settling of organic matter. The
main goal of the modeling exercise was to verify how much of theinﬂuent organic carbon is degraded, used for biomass growth, andsettled. No description of the hydraulic model was given, but it
can be assumed that the model is only capable of dealing with constant ﬂow rates. Mashauri and Kayombo (2002) concluded thatthe processes in the pond and CW system may well be deﬁned
by the Monod approach on substrate utilization and growth. The
transformation of organic carbon was found to be dominated bymineralization, which led to high growth of heterotrophic bacteria in the systems. This was also shown by the amount of organiccarbon that accumulated in the system. The model also revealed
and quantiﬁed the amount of algae growing on utilization of carbon
dioxide produced by oxidation processes in the pond.
The model presented by Mayo and Bigambo (2005) was
developed to predict nitrogen transformation in HF CWs. The
mathematical model considers the activities of biomass suspended
in the water body and bioﬁlm on aggregates and plant roots. The
state variables modeled include organic, ammonia, and nitrate
nitrogen, which were sectored in water, plants, and aggregates.
The nitrogen transformation processes considered are mineralization, nitriﬁcation, denitriﬁcation, plant uptake, release after plant
decay, and sedimentation. Environmental conditions considered
are temperature, pH, and dissolved oxygen. The experimental
CWs used for calibration and validation of the model operated
with constant water ﬂow, as the model was designed only for constant ﬂow rates. A total nitrogen removal of 48.9% was achieved,with denitriﬁcation (29.9%), plant uptake (10.2%), and net
sedimentation (8.2%) as the major pathways. In a further study,
Bigambo and Mayo (2005) investigated the signiﬁcance of the
bioﬁlm biomass in nitrogen removal. The amount of biomass
present had no inﬂuence on the removal of organic nitrogen;however, it signiﬁcantly inﬂuenced ammonia nitrogen and nitrate
nitrogen transformation. The developed CW model was also coupled with a high rate pond (Mayo and Mutamba, 2005).
Marsili-Libelli and Checchi (2005) proposed the combinationof a set of ideal reactors with a robust identiﬁcation method toapproximate the dispersed ﬂow and pollution reduction dynamics in HF CWs. The models are based on combinations of series
and parallel CSTRs of unequal volumes in series with a plug ﬂow
reactor. In each reactor, carbon removal is modeled using either
ﬁrst-order or Monod kinetics. The motivation for such a simple
model was to avoid the diﬃculties in estimating the large numberof input parameters required by other models. The estimation technique computes the parameter’s conﬁdence regions on the basis of
two diﬀering approximations producing coincident regions onlyin the case of a consistent identiﬁcation and thus allows structural
discrimination on the basis of the agreement of these regions. This
test has been shown to be very sensitive to structural and parametric perturbations and can detect model criticality not revealed byother performance indexes. Diﬀerent model structures have beencalibrated with data sets from several CWs with widely diﬀering
hydraulics and pollution removal characteristics, drawn either fromthe literature or from ﬁeld experiments performed by the authors.
The identiﬁcation method can assist in the selection of the best comwww.vadosezonejournal.org · Vol. 7, No. 2, May 2008

bination of hydraulics and kinetics to obtain robust and yet simple
models for HF CWs. It was proposed that the volumes estimated by
the identiﬁcation method can be used for wetland design.
Rousseau (2005) developed a reaction model that is coupledwith a network of CSTRs for describing water ﬂow. The CSTR
approach assumes a vertical uniform distribution of substrates,intermediates, products, and bacteria, which may not be the casefor HF CWs. Vertical mixing between the CSTRs was therefore introduced to model vertical gradients in the ﬁlter bed. The
Activated Sludge Model (ASM; Henze et al., 2000) is used to model
microbial conversions. Particulate substances are incorporated
into the model, to investigate clogging and long-term assessmentof hydraulic characteristics. Additionally, long-term simulations
incorporate meteorological data. The developed model considers
microbiological and plant-related processes aﬀecting COD andnitrogen in HF CWs. Phosphorus removal is not considered, and
it is therefore assumed that P concentrations are nonlimiting for
microbial and plant growth. Aerobic and anoxic microbial carbon
and nitrogen conversion processes are based on ASM processes. As
oxygen transfer in HF CWs is limited, anaerobic microbial processes are considered. The competition between sulfate reducing
and methanogenic bacteria is modeled as described by Kalyuzhnyi
and Fedorovich (1998). To avoid microbial inhibition due to sulﬁde accumulation, an inverse pathway was foreseen by adding
sulﬁde-oxidizing bacteria to the model. Under normal operating
conditions, it is assumed that suspended solids are completely
removed near the inlet. Wash-out of solids proportional to the ﬂow
rate is foreseen only at higher ﬂow rates. It is assumed that detachedparts of the bioﬁlm are retained within the pores and metabolized
until washed out by a peak ﬂow. The plant growth and decay model
was deliberately kept simple, describing plant growth using growth
rates that depend on ammonium and nitrate availability. Plant
material is expressed as COD, which allows smooth integration
with the COD-based model description of the microbial processes.
Other plant-related processes include decay, senescence, and physical degradation, and root oxygen loss. A pilot-scale HF CW could
be successfully modeled using this approach, although uncertainties
have been noted regarding wastewater fractionation and low measurement frequencies at the experimental plant. Rousseau (2005)
mentioned that, at the current state, the model cannot be used
for design purposes, but it might provide wetland scientists witha framework to discuss experimental results. The modeling exercise helped show the importance of sulfate processes in HF CWs.
Research needs mentioned by Rousseau (2005) are (i) the physicalreaeration process, in relation to parameters such as water velocity,
porosity, water depth, and water temperature; and (ii) the behaviorof particulate substances in the gravel matrix (i.e., ﬁltration andsettling processes, and resuspension).
Reac ve Transport Models for Variably Saturated Condi ons

The multicomponent reactive transport module CW2D
(Langergraber, 2001) was developed to describe the biochemicaltransformation and degradation processes in SSF CWs. CW2D
was incorporated into the HYDRUS variably saturated waterﬂow and solute transport program (Langergraber and Šimůnek,2006; Šimůnek et al., 2006a, 2006b). The HYDRUS program
numerically solves the Richards equation for saturated–unsaturated water ﬂow and the convection–dispersion equation for heatand solute transport. The ﬂow equation incorporates a sink term833

to account for water uptake by plant roots. The solute transport
equations consider convective–dispersive transport in the liquid
phase, diﬀusion in the gaseous phase, and nonlinear nonequilibrium reactions between the solid and liquid phases (Šimůnek etal., 2006b). The CW2D module considers 12 components and
9 processes. The components include dissolved oxygen; organic
matter (three fractions of diﬀerent degradability, i.e., readily andslowly biodegradable, and inert); ammonium, nitrite, nitrate,
and nitrogen gas; inorganic phosphorus; and heterotrophic andtwo species of autotrophic microorganisms. Organic nitrogenand organic phosphorus are modeled as nutrient contents of theorganic matter and are calculated as a percentage of COD. The
biochemical elimination and transformation processes are basedon Monod-type expressions used to describe the process rates. Allprocess rates and diﬀusion coeﬃcients are temperature dependent.
The processes considered are hydrolysis, mineralization of organic
matter, nitriﬁcation (modeled as a two-step process), denitriﬁcation, and a lysis process (as the sum of all decay and loss processes)
for the microorganisms. CW2D assumes a constant concentrationof microorganisms (and other compounds) in each ﬁnite element.
The thickness of the bioﬁlm is not considered. The mathematicalstructure of CW2D is based on the mathematical structure of
the ASMs (Henze et al., 2000). Langergraber (2005) investigatedthe plant uptake models provided by HYDRUS, which describe
nutrient uptake coupled with water uptake, and concluded that it
was possible to simulate plant uptake in high-loaded systems (e.g.,systems treating mechanically pretreated municipal wastewater).
For low-strength wastewater, the simulation results indicate thatthese models overestimate potential nutrient uptake. Oxygen
release via roots can be modeled in a way similar to nutrient
uptake (Toscano et al., 2006). For a detailed discussion of the
CW2D module, see Langergraber and Šimůnek (2005). In thelast years, several groups have used CW2D for applications besides
treating domestic wastewater, such as CWs polishing the eﬄuentof a wastewater treatment plant for reuse purposes (Toscano et al.,2006) and for treatment of combined sewer overﬂow (Dittmer etal., 2005; Meyer et al., 2006; Henrichs et al., 2007).
The recently developed model by Ojeda et al. (2006) is basedon the two-dimensional ﬁnite element code RetrasoCodeBright(RCB), which has been applied in hydrogeological studies tosimulate reactive transport of dissolved and gaseous species fornonisothermal saturated or unsaturated ﬂow domains (Rezaei etal., 2005; Saaltink et al., 2003). In RCB, a ﬁrst module calculates the ﬂow properties which are used by a second module tocalculate reactive transport. Ojeda et al. (2006) modiﬁed RCB
to include the most signiﬁcant biochemical pathways for organic
matter removal in HF CWs. The reactive transport model, basedon work by Van Cappellen and Gaillard (1996), basically consistsof six microbial kinetic reactions: hydrolysis, aerobic respiration,
nitriﬁcation, denitriﬁcation, sulfate reduction, and methanogenesis. The degradation rate associated with each of the six kinetic
equations, except for hydrolysis, is described by multiplicative
Monod-type expressions. The hydrolysis process is modeled by
an exponential function describing the total suspended solids
removal that occurs near the inlet of the bed. Phosphorus transformations, bioﬁlm growth, and oxygen leaking from macrophytes
have not been considered. Until now, the model has been used
to evaluate the importance of the organic loading rate on theremoval eﬃciencies of experimental HF CWs.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

McGechan et al. (2005b) adapted a nitrogen-cycling submodel from a soil model for CWs. Pools in the model represent
organic material with high BOD, ammonium, and nitrate.
Microbiologically controlled transformations between pools are
represented by ﬁrst-order exponential kinetics, with nitrogenﬁnally lost to the atmosphere either by ammonia volatilization or
by denitriﬁcation to gaseous nitrogen or nitrous oxide. The model
was set up to represent an experimental system with one HF andthree VF CWs. A simple model was implemented for describing
the hydraulic behavior of the CWs. Based on the assumptionof plug ﬂow, the HF bed was considered to consist of a numberof equal imaginary vertical slices, each with the same residence
time. The individual vertical slices were subdivided into equalhorizontal layers, each containing a proportion of the plant roots.
For the HF bed slices, this represents a ﬂow rate over the time
step, so the rate constant is a true time constant. However, for
the VF beds, it is simply a ﬂow over the time that the eﬄuenttakes to pass through the individual bed, so the rate constant does
not include time units. Rate constants have been selected so thatsimulated results are a reasonable approximation of measurementsfrom the experimental system. The model has been developed toassist in optimizing design parameters for new systems (including the number and sequence of reedbed types, dimensions, andﬂow rates) for various incoming contaminant concentrations
and target water-purity standards. An extension of the model(McGechan et al., 2005a) simulates oxygen transport through thewater by diﬀusion and convection and through the macrophyte
plants to the microorganisms that reside on their roots. Parametervalues have been selected to ﬁt data from a functioning reedbed
processing sewage in a rural community. This reedbed achievesa large reduction in ammonium concentration in the eﬄuent.
Since it is unlikely that so much ammonium could have been lostby volatilization alone, this implies the rapid transport of oxygento microorganisms that convert the ammonium to nitrate. The
rate found was substantially higher than could be explained by
diﬀusion and convection through the water alone, suggesting animportant role for oxygen transport via the plants to the roots.
Wanko et al. (2006b) developed a model combining chemical
reaction engineering and the mechanistic approach to transport ina one-dimensional domain. The main purpose of the model was toinvestigate oxygen transport in vertical ﬂow ﬁlters with intermittent
loading. To achieve this, oxygen was modeled in both gaseous anddissolved form. The model considers the hydrodynamics of theporous medium, development of the active biomass, transport ofsubstrate, and oxygen transfer and consumption. The vertical ﬂow
in a nonsaturated porous medium is described by the Richards
equation. Macroscopic transport is considered for one substrate andactive biomass in suspension in the porous medium liquid phase(i.e., considering biomass that is detached from the ﬁlter sand andwashed out). Biomass growth is described by Monod-type kinetic
functions. The numerical solution is based on ﬁnite elements anduses a splitting operator technique that has the advantage of aseparate solution of the convection, dispersion, and kinetics equations; each using appropriate numerical techniques. The numerical
implementation has been veriﬁed by comparing simulation results
to analytical solutions. To date, only results for calibrating the ﬂow
model have been published (Wanko et al., 2006b). Experimentswith sand columns have been performed to gain data for calibrationand validation of oxygen transport (Wanko et al., 2006a).
834

Freire et al. (2006) used a chemical engineering approach
to model a VF CW. The main aim was to develop a model for
simulating the treatment of industrial wastewater from textile
industries, including dyes. The pilot-scale VF CW was operated
in a continuous feeding mode. The study had as its main objective
to model the pilot plant and to investigate the eﬀect of higher
oxygen transfer rates in the intermittent feeding mode on treatment eﬃciency. The mechanistic model, which is in an early stageof development, considers the hydraulic and kinetic propertiesof a pilot-scale VF CW used to treat the acid orange 7 (AO7)
azo dye. The model is able to deal with intermittent and continuous feeding modes. The hydraulic submodel considers threereactors: two CSTRs and a dead zone volume. Flooding eventsand changes in bed permeability have been taken into accountto describe clogging of the ﬁlter. Transport of AO7 considers themass balance, and degradation is modeled using ﬁrst-order kinetic
rates. It was shown that the model could match the results fromthe continuous feeding mode well. However, the interpretationof the results from simulations of the intermittent feeding mode
seemed to be more diﬃcult, as the simulated removal eﬃciency
in the continuous operation was higher compared to that inintermittent feeding. This simulated result was in contrast to the
expected results. The authors concluded that the model needed
to be developed further before it can be used in practice.
F . 3. Indoor pilot-scale constructed wetlands at the technicallaboratory hall at the University of Natural Resources and AppliedLife Sciences, Vienna (BOKU).

ApplicaƟons Using CW2D
Ver cal Flow Constructed Wetlands for Wastewater Treatment
Indoor Ver cal Flow Pilot-Scale Constructed Wetlands

F . 4. Schema cExperiments have been performed at indoor VF pilot-scale
representa on of
constructed wetlands (PSCWs, Fig. 3) in the technical laborathe indoor pilot-scale
tory hall of the Institute of Sanitary Engineering at BOKU sinceconstructed wetland
2000 (Langergraber, 2003; Langergraber and Šimůnek, 2005;(values are in cen meters).
Langergraber et al., 2007). Each of the up to 10 parallel operating
PSCWs had a surface area of 1 m2. The PSCWs were loaded fourtimes per day with mechanically pretreated municipal wastewaterand have either been planted with Arundo donax L., Miscanthus
gigantean, or left unplanted. Figure 4 provides a schematic representation of a PSCW. The 60-cm main layer consists of sandy
substrate with a grain size of either 0.06–4 mm (d10 = 0.2 mm;
d60 = 0.8 mm) or 1–4 mm (d10 = 2 mm; d60 = 3 mm). The interHYDRUS. The second parameter set, obtained by inverse simulamediate layer (gravel 4–8 mm) prevents ﬁne particles from beingtion, shows results reported by Langergraber (2003). Measured
washed out into the drainage layer (gravel 16–32 mm). Diﬀerent
eﬄuent ﬂow rates and water content measurements at a depth oforganic and hydraulic loads were applied.
25 cm were used to estimate the soil hydraulic parameters using
For the numerical simulations, a vertical cross-section of thethe inverse simulation routine provided by HYDRUS. Figure 5
main layer was considered. The transport domain was discretizedcompares the measured and simulated eﬄuent ﬂow rates using
into 11 columns and 40 rows. This resulted in a two-dimensional
these two parameter sets. The simulated eﬄuent ﬂow rates are
ﬁnite element mesh consisting of 440 nodes and 780 ﬁnite elements.
already close to the data when using measured values of porosity
An atmospheric boundary condition was assigned to the top of theand saturated hydraulic conductivity (the “measured Ks + porosity”
system (a constant pressure head boundary condition of −2 cm).
parameter set). This indicates that it is necessary to measure at least
Table 1 shows two parameter sets ofT
1. Fitted soil hydraulic parameters of the van Genuchten–Mualem model for thethe van Genuchten–Mualem model (van 0.06- to 4-mm main layer (shape parameters α, N, and L [van Genuchten, 1980]; subGenuchten, 1980) of the soil hydraulic script w and s: water and substrate, respectively).
properties used in simulations of water
Parameter
Residual water Saturated water Shape parameters Saturated hydraulicﬂow. The “measured Ks + porosity” paramcontent θ rcontent θ sconduc vity Ks
NLαeter set uses measured values of porosity———— mw3 ms−3 ————
ms−1
ms h−1and saturated hydraulic conductivity Ks
Measured Ks + porosity
0.045
0.30
14.5 2.68
0.5
1.17
and the default residual water content
Langergraber (2003)
0.056
0.289
12.6 1.92
0.5
0.84
and shape parameters for sand provided bywww.vadosezonejournal.org · Vol. 7, No. 2, May 2008835

the porosity and the saturated hydraulic conductivity of the ﬁlter
material to obtain reasonable simulation results for water ﬂow.
Simulation results for water ﬂow, single-solute transport(tracer experiments), and multicomponent reactive solutetransport have been presented by Langergraber (2003) and
Langergraber and Šimůnek (2005). Table 2 shows measured
and simulated concentrations of ammonia and nitrate nitrogenfor the two diﬀerent substrates used. The diﬀerent levels of theNO3–N eﬄuent concentrations for the two substrates could be
simulated. The eﬄuent concentrations of NO3–N for the 1- to
4-mm substrate are higher because hydrolysis produced less readily
degradable organic matter that could be used for denitriﬁcationin deeper layers. In general, the simulation results showed a good
match to the measured data when the hydraulic behavior of thesystem could be well described. The good match of experimental data to reactive transport simulations can be obtained using
literature values for the CW2D model parameters (Langergraber,
2003; Langergraber and Šimůnek, 2005).

T
2. Median values of measured and simulated influent andeﬄuent concentra ons for NH4 –N and NO3–N (adapted from
Langergraber and Šimůnek, 2005).
Parameter

NH4 –N
NO3–N

Eﬄuent 0.06–4 mm
Eﬄuent 1–4 mm
substrate
substrate
Measured Measured Simula on Measured Simula on———————————— mg L−1 ————————————
60.0
0.15
0.01
1.20
0.16
3.0
38.5
41.1
50.0
63.0

calculated (from SIR) and simulated microbial biomass COD
in diﬀerent depths of the main layer for a heterotrophic lysis
rate of 0.30 d−1. When comparing measured and simulated
biomass COD in diﬀerent depths of the main layer, simulations seem to overpredict biomass COD in 1- to 5-cm depth
and underpredict biomass COD in the 5- to 10-cm depth. This
could be an indication that the inﬂuence of biomass growthon the hydraulic properties has to be considered in the model(Langergraber et al., 2007).

Comparison of Measured and Simulated Distribu on of Microbial Biomass

Outdoor Ver cal Flow Pilot-Scale Constructed Wetlands

CW2D, as well as the underlying ASMs, are based on massbalances for COD. The C and N content of the biomass wasmeasured and converted to microbial COD using factors basedon stoichiometry (Langergraber et al., 2007).
Tietz et al. (2007) applied several methods to measure
microbial and bacterial biomass parameters, including conversion of bacterial abundance (determined by microscopic direct
counts) into biomass by measurement of the cell volume, fumigation–extraction for biomass C and N, adenosine triphosphatemeasurements for biomass C, and substrate-induced respiration(SIR) for biomass C. Converting the measured data into biomass
COD results in mean values of 3400 to 5100 mg COD gDW−1(where DW is dry weight of the substrate, i.e., sand) for the ﬁrst
centimeter of the main layer, 1100 to 2600 mg COD gDW−1
from 1 to 5 cm, and 640 to 1400 mg COD gDW−1 from 5 to 10
cm, respectively. Most of the biomass could be found in the upper
10 cm of the VF bed (Langergraber et al., 2007).
Simulated microbial biomass COD in the ﬁrst centimeterof the main layer is between 5600 and 3400 mg COD gDW−1(the range of the measured values) when using heterotrophic
lysis rates of between 0.25 and 0.35 d−1. Figure 6 compares

Outdoor experiments have been performed at experimental
CWs located at the wastewater treatment plant in Ernsthofen(Lower Austria, Fig. 7). The experimental plant, in operation
since 2003, consists of three parallel operating VF beds with asurface area of about 20 m2 each operated with intermittent loading (Langergraber, 2007). The organic loads for Beds 1 through3 have been 20, 27, and 40 g COD m−2 d−1, which correspondto hydraulic loading rates of 32.2, 43.0, and 64.7 mm d−1. The
main layer of the ﬁlter consists of 50 cm sandy substrate (gravel
size 0.06–4 mm, d10 = 0.2 mm; d60 = 0.8 mm). The beds areplanted with common reed (Phragmites australis).
The data gained from this experiment have been used
to verify the temperature model incorporated into CW2D
(Langergraber, 2007). The width of the transport domain in the
numerical simulations was 4 m, and its depth was 0.8 m. The
transport domain itself was discretized into 31 columns and
33 rows, resulting in a two-dimensional ﬁnite element meshconsisting of 1023 nodes and 1920 triangular ﬁnite elements.
An atmospheric boundary condition was assigned to the top of
the system representing the inﬂuent distribution system, and a

F . 5. Measured and simulated eﬄuent flow rates for a single
feeding of 10 L (“Measured Ks+porosity” and “Langergraber(2003)” represent the parameter sets from Table 1, according toLangergraber and Šimůnek, 2005).
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

Influent

F . 6. Calculated and simulated microbial biomass chemical oxygen
demand (COD) in diﬀerent depths of the main layer (adapted from
Langergraber et al., 2007) (DW: dry weight).
836

F . 7. Outdoor experimental constructed wetlands in Ernsthofen (Lower Austria).

F . 8. Measured and simulated chemical oxygen demand (COD) and NH4–N
eﬄuent concentra ons using the modified parameter set (Langergraber, 2007).

constant pressure head boundary condition (constant head of −4cm) was assigned to one side of the drainage layer.
Results for water ﬂow comparing measured and simulated
eﬄuent ﬂow rates and the cumulative eﬄuent of the experimental
CW are shown in Langergraber (2007). Simulated eﬄuent ﬂow
rate and cumulative eﬄuent ﬂow match the measured data well.
Using the calibrated ﬂow model, the eﬄuent concentrations during
summer could be simulated using the standard CW2D parameterset (Langergraber and Šimůnek, 2005). However, measured COD
and NH4–N eﬄuent concentrations at low temperatures could
not be simulated because hydrolysis and nitriﬁcation at low temperatures were overpredicted. The standard CW2D parameter set
considers temperature dependencies only for maximum growth,decay, and hydrolysis rates. Because experience from activatedsludge systems showed that considering temperature dependencies for the half-saturation constants of hydrolysis and nitriﬁcation
improves the model predictions, they have therefore been included
in ASM parameter sets (e.g., Bornemann et al., 1998). Figure 8
shows measured and simulated COD and NH4–N eﬄuent concentrations using a modiﬁed parameter set that includes temperaturedependencies for the half-saturation constants of hydrolysis andnitriﬁcation. Using this modiﬁed parameter set, it was possible tosimulate the COD and NH4–N eﬄuent concentrations at low
temperatures (Langergraber, 2007).

velocity in the ﬁlter itself. For simulation purposes, the maximum
allowed eﬄuent ﬂow rate had to be implemented (Langergraberand Šimůnek, 2006).
The ﬁrst simulations were performed for lab-scale columns
with a diameter of 19 cm (Dittmer et al., 2005). The 100-cm
main layer consisted of sand (grain size 0–2 mm, eﬀective size 0.1
mm, coeﬃcient of uniformity 2.8). The test site was comprised ofsix columns all fed identically. The columns were loaded once aweek with 15.7 L of synthetic sewer, which corresponds to a loading rate of 0.5 m per event. The ﬁltration rate of Columns 1 to
4 was controlled manually by a throttle valve, whereas Columns
5 and 6 had free drainage. For the simulation, a structured twodimensional ﬁnite element mesh with 105 rows and 3 columns
was used, resulting in 315 nodes and 416 elements. The storage
volume was modeled by a virtual layer with a pore volume of
100% and a residual water content of 0%. On the top of thevirtual layer, an atmospheric boundary condition was applied,
while a free drainage or seepage face with maximum allowed ﬂux
was used at the bottom (Dittmer et al., 2005).
Figure 9 compares simulated eﬄuent ﬂow rates for a column
with free drainage and a limited seepage face ﬂux boundary condition. The maximum ﬂow in the controlled eﬄuent rate case(right) corresponds to 1 L/h. While in the unrestricted ﬂow case,the bulk of the water passed through the ﬁlter bed in about 1 h,in the ﬂow restricted case this happened in about 15 h.
Figure 10 shows measured and simulated breakthroughcurves of a tracer experiment performed for Column 6 with freedrainage. Figure 11 shows breakthrough curves for Column 4
and Column 2, with controlled eﬄuent rates of 5 and 1 L/h,
respectively. Using the transport parameters obtained from thefree drainage simulations (Column 6), the tracer experimentswith controlled eﬄuent rate could be simulated well.
For multicomponent reactive transport simulations, goodresults could be achieved when considering the diﬀerent composition of CSO compared to domestic wastewater (Dittmer
et al., 2005) and considering adsorption for COD fractions(Henrichs et al., 2007). Figure 12 shows, as an example, thesimulated eﬄuent concentrations of COD and its fractions fora lab-scale experiment with a controlled eﬄuent rate. Using theparameters obtained from lab-scale experiments, it was also pos-

Ver cal Flow Constructed Wetlands for Combined
Sewer Overflow Treatment

CW2D was used to model VF CWs for the treatment of
combined sewer overﬂow (CSO) by Dittmer et al. (2005), Meyer
et al. (2006) and Henrichs et al. (2007). Some of the diﬀerences
between using wetlands for wastewater and for CSO treatments
are the loading regime and the quality parameters of the inﬂow.
For CSO treatment, the succession of loading events and dry
periods is characterized by the stochastic nature of rainfall and therunoﬀ behavior of the catchment area. On the one hand, extreme
cases involve a permanent loading for weeks, while several months
without any loading event are possible at the other extreme. To
reach the main treatment objective, the detention and reductionof peak ﬂows (Dittmer et al., 2005), throttle valves are applied
to limit the maximum eﬄuent ﬂow rate and therefore the ﬂow
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

837

F . 9. Comparison ofsimula ons withfree drainage
boundary condion (le ) andcontrolled eﬄuent rate (right).

F . 10. Measured and simulated breakthrough curves of Column 6
with free drainage (Di mer et al., 2005).

F . 11. Measured and simulated breakthrough curves of Column 4
and Column 2 for a controlled eﬄuent rate of 5 and 1 L h−1, respecvely (Langergraber and Šimůnek, 2006).

sible to successfully simulate single events of ﬁeld experiments.
The consideration of the adsorption process for COD fractions
is important for controlled eﬄuents as the ﬂow velocity in theﬁlter is low. For vertical ﬁlters with free drainage, adsorption
processes for COD fractions can be neglected due to the highﬂow velocities in the ﬁlter. Henrichs et al. (2007) further showedthat problems occur when simulating long-term experiments that
include long dry periods. These problems might be caused by
diﬀerent degradation processes during dry periods that may not
be represented well in CW2D.

simulate the HF beds, a longitudinal section (3.0-m length and
0.6-m depth, with 16 rows and 25 columns, resulting in 400
nodes and 720 ﬁnite elements) was used; the VF beds used a vertical cross-section (1.5-m width and 0.75-m depth; 25 rows and
16 columns; 400 nodes and 720 ﬁnite elements). Atmosphericboundary conditions have been used for the inﬂuent boundaries, and constant pressure head boundary conditions for theeﬄuent boundaries.

Two Stage Horizontal Flow and Ver cal Flow
Pilot-Scale Constructed Wetlands

The experimental plant used for these experiments is
located in San Michele di Ganzaria (Eastern Sicily, Italy)
and has been in operation since June 2004 (Toscano et al.,2006). Each of the four parallel lines consists of two SSF
beds in series (Fig. 13). The plant was designed to treat
eﬄuents from both the primary and secondary settlersof the conventional wastewater treatment plant (tricking
ﬁlter). For each line, the ﬁrst stage consists of a HF bed(HF1–4, depths 0.6 m, volcanic gravel, 10–15 mm), while
the second stage uses an intermittently loaded VF bed for
two lines (VF5–6, 0.5 m main layer, volcanic sand, 0.06–4mm) and a HF bed for the other two (HF7–8). Each bed
has a rectangular shape with a surface area of 4.5 m2 (1.5
× 3.0 m). Lines 1 and 3 are planted with reed (Phragmites
australis) whereas Lines 2 and 4 are unplanted.
Simulations were performed for Lines 1 and 2 ofthe pilot plant for secondary and tertiary treatment. To
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

F . 12. Eﬄuent concentra ons of chemical oxygen demand (COD) and its fracons with COD adsorp on and controlled eﬄuent rate (adapted from Henrichs
et al., 2007; O2: dissolved oxygen; CR, CS, CI: readily biodegradable, slowly
biodegradable, and inert organic ma er, respec vely).
838

the planted bed (HF1) were signiﬁcantly lower comparedwith the unplanted bed (HF2). By considering the eﬀects ofplants (i.e., oxygen release by roots and nutrient uptake), it
was possible to model the diﬀerences in the NH4–N eﬄuent
concentrations in planted and unplanted HF beds.
Two-Stage Pilot-Scale Downflow and Upflow
Constructed Wetland

Simulation results for a two-stage pilot-scale CW fortreating surface water (Langergraber, 2003; Langergraberand Šimůnek, 2005) are included here to show an application
for a more complicated ﬂow domain geometry. The
F . 13. Schema c sketch of pilot plant in San Michele di Ganzaria (Toscanosystem was designed for the treatment of heavily polluted
et al., 2006).
surface water from shallow lakes in China (Grosse et al.,The ﬂow model has been calibrated using eﬄuent ﬂow rate1999) with inﬂuent concentrations of 5 to 10 mg NH4–N
measurements of the VF beds as well as the measured porosity
L−1. The total surface area of the two-stage system was 2 m2,
and saturated hydraulic conductivities of the ﬁlter material. The
divided into downﬂow and upﬂow chambers (surface area is 1 m2times of the loadings of the VF beds have been recorded. Theeach). The inlet was situated on the top of the downﬂow chamber;inﬂuent ﬂow rate for the HF beds was calculated by providing
the eﬄuent was collected on the top of the upﬂow chamber by
the ﬁxed amount of a single loading for VF beds and the recordedmeans of perforated pipes. There was a separating wall between
interval between two successive loadings. Figure 14 comparesthe downﬂow and upﬂow chambers (Perﬂer et al., 1999; Fig. 16).
measured and simulated breakthrough curves for a tracer experiThe main layer (depths of 55 cm downﬂow and 45 cm upﬂow)
ment on Line 2 (HF2 + VF6). The simulation results matched
was ﬁlled with sand (gravel size 0.06–4 mm). The hydraulic loadthe measured data welling rate was 50 L every 2 h (i.e., 600 L d−1).
For secondary treatment, simulations were performed for
An unstructured two-dimensional ﬁnite element meshLine 1 (with vegetation) and Line 2 (without vegetation). The
was used for simulations consisting of 1135 nodes and 2057
eﬀects of plant roots, oxygen release, and nutrient uptake haveelements (Fig. 17). The inlet distribution pipe at the top ofbeen taken into account. The root zone in the HF bed was set tothe downflow chamber was modeled as a time dependent
be half the depth of the bed (30 cm) (Rousseau, 2005). Oxygen
atmospheric boundary condition. Outlet (eﬄuent) pipes wererelease was assumed to be 5 g O2 m−2 d−1) (Vymazal et al., 1998).
modeled as circular openings in the transport domain with a
The evapotranspiration rate was calculated to be 7.4 mm d−1,
seepage face boundary condition.
and it was assumed that ammonia and nitrate are taken up by
Langergraber (2003) showed that taking into account stagplants with the same preference. Figure 15 shows the simulated
nant (immobile) regions of the pore water resulted in a muchoxygen concentration proﬁle in HF1 with plants. The oxygen
better match of the simulated concentration data for the tracer
concentration increased up to about 0.2 mg L−1 in the upper
experiment, compared with instances when only mobile pore
half of HF1, with higher concentrations of up to 0.3 mg L−1
water was considered. Langergraber and Šimůnek (2005) showedvery close to the inlet.
several simulation results using CW2D. Figure 18 shows as an
Table 3 shows measured inﬂuent and eﬄuent, and simulated
example of simulated and measured concentrations of NH4–N
eﬄuent concentrations for HF1 and HF2. No diﬀerence could be
and NO3–N in a cross-section of the downﬂow chamber.
observed for the eﬄuent COD concentrations of the planted andunplanted bed. However, the eﬄuent NH4–N concentrations ofDiscussion
Until now, very few models to simulate removal processes
in SSF CWs have been developed. Of these models, several areapplicable to HF CWs, as they only consider saturated waterﬂow using, for example, a tanks-in-series approach (Chen etal., 1999; Wynn and Liehr, 2001; Marsili-Libelli and Checchi,
2005; Rousseau, 2005) or are only applicable for constant ﬂow
rates (Mashauri and Kayombo, 2002; Mayo and Bigambo,
2005). Three of these six models considered only carbon trans-

F . 14. Measured and simulated breakthrough curves for a tracer experiment on Line 2 (HF2 + VF6) (Toscano et al., 2006).
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

839

F . 15. Simulated oxygen concentra on profile (mg L−1) in
HF1 with plants (root depth is half of bed depths).

T
3. Measured influent and eﬄuent concentra ons (mean values and standard devia on in parentheses) and simulated eﬄuent concentra ons (mg L−1) for
HF1 and HF2.

oxygen transport, whereas the model of Freire et al.
(2006) only describes the removal of the dye AO7.
Like CW2D/HYDRUS, the models developed byParameter
COD
NH4 –N
NO3–N
NO2–N
Ojeda et al. (2006) and Wanko et al. (2006b) use the————————— mg L−1 —————————
Richards equation for describing variably saturated
Influent
Measured
63.9 (3.3) 22.1 (1.9)
0.18 (0.38) 0.34 (0.10)
ﬂow, whereas the other models use either diﬀerent
Eﬄuent HF1 (planted)
Measured
47.5 (2.8) 16.6 (1.2)
0.06 (0.05) 0.18 (0.17)
layers (McGechan et al., 2005a, 2005b) or a combinaSimulated
44.4
17.2
0.02
2.82
tion
of CSTRs and dead zones (Freire et al., 2006).
Eﬄuent HF2 (unplanted) Measured
48.2 (4.1) 20.4 (0.6)
0.02 (0.01) 0.12 (0.11)
The multicomponent reactive transport module
Simulated
50.6
21.5
6E-05
1E-04
CW2D (Langergraber, 2001) is able to describe thebiochemical transformation and degradation processes
formation processes (Chen et al., 1999; Mashauri and Kayombo,
for organic matter, nitrogen, and phosphorus in SSF CWs. The
2002; Marsili-Libelli and Checchi, 2005), one only nitrogen transexperience showed that simulation results match the measured data
formations (Mayo and Bigambo, 2005), and the remaining two
when the hydraulic behavior of the system can be described well.
both carbon and nitrogen transformations (Wynn and Liehr, 2001;
A good match of experimental data to reactive transport simulaRousseau, 2005). The most advanced of these models is the one
tions can then be obtained using literature values for the CW2D
developed by Rousseau (2005), who implemented a reaction model
model parameters (Langergraber, 2003; Langergraber and Šimůnek,presented in matrix notation based on the mathematical structure2005). For practical applications, it is advisable to measure at leastof the ASMs that includes the description of aerobic, anoxic, andthe porosity and saturated hydraulic conductivity of the ﬁlter mateanaerobic processes. However, until now Rousseau’s model (2005)
rial to obtain reasonable simulation results for water ﬂow.
has only been applied to one pilot-scale HF CW.
Being incorporated into the HYDRUS variably saturated
Besides CW2D (Langergraber, 2001), there have been fourwater ﬂow and solute transport program (Langergraber andmodels described in the literature that consider vadoze zone proŠimůnek, 2006), CW2D is easy to use and has been applied bycesses. The other models are in a very early stage of development.
other research groups, for example, for CWs polishing the eﬄuThe most advanced one is the model developed by Ojeda et al.
ent of a wastewater treatment plant for reuse purposes (Toscano(2006), which considers processes aﬀecting solids, organic matter,et al., 2006) and treatment of combined sewer overﬂow (Dittmer
nitrogen, and sulfur. The model was developed primarily for HF
et al., 2005; Meyer et al., 2006; Henrichs et al., 2007).
CWs, but, because of the underlying ﬂow model, it is also capable
Although CW2D has been has been applied several times,of simulating VF CWs. The model developed by McGechan etthere are still a number of research needs that must be addressedal. (2005a, 2005b) considers pools of organic matter, ammobefore it can be used as a reliable design tool for CWs. As alreadynium, and nitrate, as well as oxygen. Microbiologically controlledsuggested by Langergraber (2003), a module that will enable thetransformations are deﬁned between pools. The model developed
description of the inﬂuence of particulates and biomass growth on
by Wanko et al. (2006b) considers organic matter removal andthe hydraulic properties of the ﬁlter needs to be incorporated. The
module should be able to describe pore size reductions due to thesettling of suspended solids and bacteria growth. Langergraber etal. (2007) conclude that the inﬂuence of biomass growth on thehydraulic properties must be included for a better match of simulated and measured data on biomass distribution in a VF ﬁlter.
Another research need is the development of experimental
methods that estimate CW2D model parameters, as none are currently available (Langergraber and Šimůnek, 2005). Andreottola
et al. (2007) presented a respirometric technique to measure theF . 16. Schema c representa on of the two-stage ver cal flow
constructed wetland (values are in cm).

F . 18. Simulated and measured concentra ons of NH4–N andNO3–N in a cross-sec on of the downflow chamber (Langergraberand Šimůnek, 2005).

F . 17. Two-dimensional mesh of the two-stage ver cal flow constructed wetland.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008840

References

kinetics of organic matter oxidation and nitriﬁcation using columns
that simulate cores of VF CWs. The following kinetic parameters
can be extracted from the respirograms: maximum oxidation rateof readily biodegradable COD, maximum nitriﬁcation rate, andendogenous respiration. The method proposed by Andreottola etal. (2007) seems promising for use with the CW2D module.
For the simulation of VF beds for CSO treatment, thedegradation processes in dry periods should be diﬀerent thanthose modeled in CW2D for organic matter degradation andnitriﬁcation (Henrichs et al., 2007). These diﬀerences have to beinvestigated in more detail and possibly included within CW2D.
As anaerobic processes play a major role in HF CWs, theinclusion of anaerobic processes as described in the models ofRousseau (2005) and Ojeda et al. (2006) may have to be considered as well.

Andreottola, G., E. Oliveira, P. Foladori, R. Peterlini, and G. Ziglio. 2007. Respirometric techniques for assessment of biological kinetics in constructedwetlands. Water Sci. Technol. 56(3):255–261.
Bigambo, T., and A.W. Mayo. 2005. Nitrogen transformation in horizontal subsurface ﬂow constructed wetlands: I. Eﬀect of bioﬁlm. Phys. Chem. Earth
30:668–672.
Bornemann, C., J. Londong, M. Freund, O. Nowak, R. Otterpohl, and T. Rolfs.
1998. [A guide for the dynamic simulation of aeration plants with the
IAWQ’s Activated Sludge Model No. 1]. (In German.). Korrespondenz
Abwasser 45(3):455–462.
Brix, H. 1997. Do macrophytes play a role in constructed treatment wetlands?
Water Sci. Technol. 35(5):11–17.
Brix, H., and N.H. Johansen. 2004. Guidelines for vertical ﬂow constructedwetland systems up to 30 PE. (In Danish.) Økologisk Byfornyelse og
Spildevandsrensning No. 52. Miljøstyrelsen, Miljøministeriet, Copenhagen, Denmark.
Chazarenc, F., G. Merlin, and Y. Gonthier. 2003. Hydrodynamics of horizontal
subsurface ﬂow constructed wetlands. Ecol. Eng. 21:165–173.
Chen, S., G.T. Wang, and S.K. Xue. 1999. Modeling BOD removal in constructed wetlands with mixing cell method. J. Environ. Eng. 125:64–71.
Dittmer, U., D. Meyer, and G. Langergraber. 2005. Simulation of a subsurface
vertical ﬂow constructed wetland for CSO treatment. Water Sci. Technol.
51(9):225–232.
DWA-A 262. 2006. Grundsätze für Bemessung, Bau und Betrieb von Pﬂanzenkläranlagen mit bepﬂanzten Bodenﬁltern zur biologischen Reinigung
kommunalen Abwassers. (In German.) Deutsche Vereinigung für Wasserwirtschaft, Abwasser und Abfall e.V., Hennef, Germany.
Freire, F.G., L.C. Davies, A. Vacas, I. Pedro, J.M. Novais, and S. Martins-Dias.
2006. Continuous and intermittent loading of a vertical ﬂow constructedwetland for an azo dye treatment. p. 1501–1509. In Proc. IWA SpecializedGroup Conf. on Wetland Systems for Water Pollution Control, Vol. 3,
10th, Lisbon, Portugal. 23–29 Sept. 2006. Int. Water Assn., London.
García, J., J. Chiva, P. Aguirre, E. Alvarez, J.P. Sierra, and R. Mujeriego. 2004. Hydraulic behaviour of horizontal subsurface ﬂow constructed wetlands with
diﬀerent aspect ratio and granular medium size. Ecol. Eng. 23:177–187.
Grosse, W., F. Wissing, R. Perﬂer, Z. Wu, J. Chang, and Z. Lei. 1999. Biotechnological approach to water quality improvement in tropical and subtropical areas
for reuse and rehabilitation of aquatic ecosystems. Final report, INCO-DC
Project Contract ERBIC18CT960059. INCO-DC, Cologne, Germany.
Haberl, R., S. Grego, G. Langergraber, R.H. Kadlec, A.R. Cicalini, S. Martins
Dias, J.M. Novais, S. Aubert, A. Gerth, H. Thomas, and A. Hebner. 2003.
Constructed wetlands for the treatment of organic pollutants. J. Soils Sediments 3:109–124.
Henrichs, M., G. Langergraber, and M. Uhl. 2007. Modelling of organic matter degradation in constructed wetlands for treatment of combined sewer
overﬂow. Sci. Total Environ. 380:196–209.
Henze, M., W. Gujer, T. Mino, and M.C.M. van Loosdrecht. 2000. Activatedsludge models ASM1, ASM2, ASM2D and ASM3. IWA Scientiﬁc andTechnical Rep. 9. Int. Water Assn., London.
Kadlec, R.H. 2000. The inadequacy of ﬁrst-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 Scientiﬁc and Technical Rep. 8. Int. Water
Assn., London.
Kalyuzhnyi, S.V., and V.V. Fedorovich. 1998. Mathematical modelling of competition between sulphate reduction and methanogenesis in anaerobic reactors. Bioresour. Technol. 65:227–242.
Katyal, A.K., J.J. Kaluarachchi, and J.C. Parker. 1991. MOFAT: A two-dimensional ﬁnite element program for multiphase ﬂow and multi-component
transport: Program documentation and user’s guide. Center for Environ.
and Hazardous Material Studies, Virginia Polytechnic Inst. and State
Univ., Blacksburg.
Kayser, K., and S. Kunst. 2005. Processes in vertical-ﬂow reed beds: Nitriﬁcation,oxygen transfer, and soil clogging. Water Sci. Technol. 51(9):177–184.
Langergraber, G. 2001. Development of a simulation tool for subsurface ﬂow
constructed wetlands. Wiener Mitteilungen 169, Vienna, Austria.

Summary and Conclusions
The previous sections have reviewed numerical models
describing processes in subsurface ﬂow constructed wetlands.
Several authors aimed for modeling the hydraulic behavior of thesystem including tracer experiments only. There are few publishedmodels that are able to model the biochemical transformationand degradation processes that occur in the complex system subsurface ﬂow constructed wetland.
Models describing reactive transport under water ﬂow saturated conditions can only be used to model horizontal ﬂow
constructed wetlands. Six models have been reviewed in thiscategory that use either a series or a network of ideal reactors(continuously stirred tank reactors and/or plug ﬂow reactors), or
are only applicable for constant ﬂow rates to model water ﬂow.
Transformation and removal processes are described for organic
matter and/or nitrogen.
Modeling vertical ﬂow constructed wetlands with intermittent loading requires transient variably saturated ﬂow models asthese systems are highly dynamic, which adds to the complexityof the overall system. Five models of diﬀering complexities havebeen published that consider variably saturated ﬂow and reaction models for constructed wetlands. Four of these models arein the rather early stages of development. For the multicomponent reactive transport module CW2D, several applications havebeen published that include treatment of domestic wastewaterwith vertical ﬂow constructed wetlands, treatment of combined
sewer overﬂow, and polishing of treated wastewater. The usage
of CW2D, an extension of the variably saturated water ﬂow andsolute transport program HYDRUS, was also described.
Based on the experience gained by applying CW2D, the following conclusions can be drawn: (i) a good match to measured
data can be achieved when the hydraulic behavior of the system iscorrectly described; (ii) it is thus suggested, for practical applications, to measure at least the porosity and the saturated hydraulicconductivity of the porous ﬁlter material; (iii) the description of
the inﬂuence of particulates and biomass growth on the hydraulicproperties of the ﬁlter in the model would increase its prediction
quality; (iv) research is needed to develop experimental methods for
estimating the kinetic parameters of CW2D; (v) degradation processes in very dry periods occurring when treating combined sewer
overﬂow have to be investigated; and (vi) anaerobic processes playa major role in horizontal ﬂow constructed wetlands and should
therefore be considered when modeling horizontal ﬂow beds.
www.vadosezonejournal.org · Vol. 7, No. 2, May 2008841

Langergraber, G. 2003. Simulation of subsurface ﬂow constructed wetlands: Results and further research needs. Water Sci. Technol. 48(5):157–166.
Langergraber, G. 2005. The role of plant uptake on the removal of organic matter and nutrients in subsurface ﬂow constructed wetlands: A simulation
study. Water Sci. Technol. 51(9):213–223.
Langergraber, G. 2007. Simulation of the treatment performance of outdoor
subsurface ﬂow constructed wetlands in temperate climates. Sci. Total Environ. 380:210–219.
Langergraber, G., and R. Haberl. 2001. Constructed wetlands for water treatment. Minerva Biotecnol. 13:123–134.
Langergraber, G., and R. Haberl. 2004. Application of constructed wetland
technology in EcoSan systems. In Proc. IWA World Water Congress, [CDROM], 4th, Marrakech, Morocco. 19–24 Sept. 2004. Int. Water Assn.,
London.
Langergraber, G., and J. Šimůnek. 2005. Modeling variably saturated water ﬂow
and multicomponent reactive transport in constructed wetlands. Vadose
Zone J. 4:924–938.
Langergraber, G., and J. Šimůnek. 2006. The multi-component reactive transport module CW2D for constructed wetlands for the HYDRUS software
package. Hydrus Software Series 2. Dep. of Environ. Sciences, Univ. of
California, Riverside.
Langergraber, G., A. Tietz, and R. Haberl. 2007. Comparison of measured andsimulated distribution of microbial biomass in subsurface vertical ﬂow
constructed wetlands. Water Sci. Technol. 56(3):233–240.
Małoszewski, P., P. Wachniew, and P. Czupryński. 2006. Study of hydraulicparameters in heterogeneous gravel beds: Constructed wetland in Nowa
Słupia (Poland). J. Hydrol. 331:630–642.
Marsili-Libelli, S., and N. Checchi. 2005. Identiﬁcation of dynamic models forhorizontal subsurface constructed wetlands. Ecol. Modell. 187:201–218.
Mashauri, D.A., and S. Kayombo. 2002. Application of the two coupled models
for water quality management: Facultative pond cum constructed wetlandmodels. Phys. Chem. Earth 27:773–781.
Mayo, A.W., and T. Bigambo. 2005. Nitrogen transformation in horizontal subsurface ﬂow constructed wetlands: I. Model development. Phys. Chem.
Earth 30:658–667.
Mayo, A.W., and J. Mutamba. 2005. Modelling nitrogen removal in a coupled
HRP and unplanted horizontal ﬂow subsurface gravel bed constructedwetland. Phys. Chem. Earth 30:673–679.
McGechan, M.B., S.E. Moir, K. Castle, and I.P.J. Smit. 2005a. Modelling oxygen transport in a reedbed-constructed wetland puriﬁcation system fordilute eﬄuents. Biosystems Eng. 91:191–200.
McGechan, M.B., S.E. Moir, G. Sym, and K. Castle. 2005b. Estimating inorganic and organic nitrogen transformation rates in a model of a constructed wetland puriﬁcation system for dilute farm eﬄuents. Biosystems
Eng. 91:61–75.
Meyer, D., G. Langergraber, and U. Dittmer. 2006. Simulation of sorption processes in subsurface vertical ﬂow constructed wetlands for CSO treatment.
p. 599–609. In Proc. IWA Specialized Group Conf. on Wetland Systems
for Water Pollution Control, Vol. 1 , 10th, Lisbon, Portugal. 23–29 Sept.
2006. Int. Water Assn., London.
Molle, P., A. Liénard, C. Boutin, G. Merlin, and A. Iwema. 2005. How to treat
raw sewage with constructed wetlands: An overview of the French systems.
Water Sci. Technol. 51(9):11–21.
Ojeda, E., J. Caldentey, and J. García. 2006. 2D simulation model for evaluating
biogeochemical pathways involved in organic matter removal in horizontal
subsurface ﬂow constructed wetlands. p. 1405–1413. In Proc. IWA Specialized Group Conf. on Wetland Systems for Water Pollution Control, Vol. 3,
10th, Lisbon, Portugal. 23–29 Sept. 2006. Int. Water Assn., London.
ÖNORM B 2505. 2005. [Subsurface-ﬂow constructed wetlands: Application,
dimensioning, installation, and operation]. (In German.) Österreichisches
Normungsinstitut, Vienna, Austria.
Pastor, R., C. Benqlilou, D. Paz, G. Cardenas, A. Espuña, and L. Puigjaner.
2003. Design optimisation of constructed wetlands for wastewater treatment. Resour. Conserv. Recycling 37:193–204.
Perﬂer, R., J. Laber, G. Langergraber, and R. Haberl. 1999. Constructed wetlands for rehabilitation and reuse of surface waters in tropical and subtropical areas: First results from small-scale plots using vertical ﬂow beds. Water
Sci. Technol. 40(3):155–162.
Rezaei, M., E. Sanz, E. Raeisi, C. Ayora, E. Vázquez-Suñé, and J. Carrera. 2005.
Reactive transport modeling of calcite dissolution in the fresh-salt water

www.vadosezonejournal.org · Vol. 7, No. 2, May 2008

mixing zone. J. Hydrol. 311:282–298.
Rousseau, D.P.L. 2005. Performance of constructed treatment wetlands: Modelbased evaluation and impact of operation and maintenance. Ph.D. diss.,
Ghent Univ., Belgium.
Rousseau, D.P.L., P.A. Vanrolleghem, and N. De Pauw. 2004. Model-based
design of horizontal subsurface ﬂow constructed treatment wetlands: A
review. Water Res. 38:1484–1493.
Saaltink, M.W., C. Ayora, P.J. Stuyfzand, and H. Timmer. 2003. Analysis of adeep well recharge experiment by calibrating a reactive transport modelwith ﬁeld data. J. Contam. Hydrol. 65:1–18.
Schwager, A., and M. Boller. 1997. Transport phenomena in intermittent ﬁlters.
Water Sci. Technol. 35(6):13–20.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 1998. The HYDRUS-1D
software package for simulating the one-dimensional movement of water,
heat, and multiple solutes in variably-saturated media. Version 2.0. U.S.
Salinity Lab., USDA-ARS, Riverside, CA.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 2006a. The HYDRUS software package for simulating the two- and three-dimensional movementof water, heat, and multiple solutes in variably-saturated media: Technicalmanual. Version 1.0. PC-Progress, Prague, Czech Republic.
Šimůnek, J., M. Šejna, and M.Th. van Genuchten. 2006b. The HYDRUS software package for simulating the two- and three-dimensional movement ofwater, heat, and multiple solutes in variably-saturated media: User manual.
Version 1.0. PC-Progress, Prague, Czech Republic.
Stein, O.R., J.A. Biederman, P.B. Hook, and W.C. Allen. 2006. Plant species
and temperature eﬀects on the k–C* ﬁrst-order model for COD removal
in batch-loaded SSF wetlands. Ecol. Eng. 26:100–112.
Tanner, C.C. 2001. Plants as ecosystem engineers in subsurface-ﬂow treatment
wetlands. Water Sci. Technol. 44(11–12):9–17.
Tietz, A., G. Langergraber, K. Sleytr, A. Kirschner, and R. Haberl. 2007. Characterization of microbial biocoenosis in vertical subsurface ﬂow constructedwetlands. Sci. Total Environ. 380:163–172.
Tomenko, V., S. Ahmed, and V. Popov. 2007. Modelling constructed wetland
treatment system performance. Ecol. Modell. 205:355–364.
Toscano, A., G. Langergraber, and G.L. Cirelli. 2006. Simulation of hydraulics and pollutant removal of a pilot-scale two-stage constructed wetlands
functioning as secondary or tertiary treatment. p. 1303–1311. In Proc.
IWA Specialized Group Conf. on Wetland Systems for Water Pollution
Control, Vol. 2, 10th, Lisbon, Portugal. 23–29 Sept. 2006. Int. Water
Assn., London.
Van Cappellen, P., and J.-F. Gaillard. 1996. Biogeochemical dynamics in aquatic
sediments. p. 335–376. In P.C. Lichtner, C.I. Steefel, and E.H. Oelkers(ed.) Reviews in mineralogy and geochemistry. Vol. 34. Mineralogical Soc.
of America, Washington, DC.
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.
Vymazal, J., H. Brix, P.F. Cooper, M.B. Green, and R. Haberl (ed.). 1998. Constructed wetlands for wastewater treatment in Europe. Backhuys, Leiden,
The Netherlands.
Wanko, A., N. Forquet, R. Mose, and A.G. Sadowski. 2006a. Oxygen transfer
rates estimation for sand bed design: Model calibration and validation. p.
1595. In Proc. IWA Specialized Group Conf. on Wetland Systems for Water Pollution Control, Vol. 3, 10th, Lisbon, Portugal. 23–29 Sept. 2006.
Int. Water Assn., London.
Wanko, A., R. Mose, J. Carrayrou, and A.G. Sadowski. 2006b. Simulation of
biodegradation in inﬁltration seepage: Model development and hydrodynamic calibration. Water Air Soil Pollut. 177:19–43.
Werner, T.M., and R.H. Kadlec. 2000. Wetland residence time modeling. Ecol.
Eng. 15:77–90.
Wynn, M.T., and S.K. Liehr. 2001. Development of a constructed subsurfaceﬂow wetland simulation model. Ecol. Eng. 16:519–536.

842