S CI EN C E OF TH E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3a v a i l a b l e a t w w w. s c i e n c e d i r e c t . c o mw w w. e l s e v i e r. c o m / l o c a t e / s c i t o t e n v

Recent developments in numerical modelling of subsurface
flow constructed wetlands
Günter Langergraber a,⁎, David Giraldi b , Javier Mena c , Daniel Meyer d , Miguel Peña e ,Attilio Toscano f , Alessandro Brovelli g , E. Asuman Korkusuz aaInstitute of Sanitary Engineering and Water Pollution Control, University of Natural Resources and Applied Life Sciences, Vienna (BOKU),
Muthgasse 18, A-1190 Vienna, Austriab
Department of Civil Engineering, University of Pisa, Via Gabba 22, 56122 Pisa, Italyc
Department R+D+i, Water Research and Development Centre for the Recuperation of the Environment (CIDRA), C/Calidad 3,
Poligono Industrial Sur, 13250 Daimiel (Ciudad Real), SpaindInstitute of Urban Water Management, Technical University of Kaiserslautern, Paul-Ehrlich-Strasse, 14-316, D-67663 Kaiserslautern,
GermanyeUniversidad del Valle, Instituto Cinara, A.A 25157, Cali, ColombiafDepartment of Agricultural Engineering-Hydraulics division, University of Catania, Via S. Sofia 100, 95123 Catania, Italyg
Laboratoire de technologie écologique, Institut des sciences et technologies de l'environnement, Station No. 2,
Ecole Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland

AR TIC LE I N FO

ABS TR ACT

Article history:

Numerical modelling of subsurface flow constructed wetlands (CWs) gained increasing

Received 29 February 2008

interest during the last years. The main objective of the modelling work is, on the one hand, toReceived in revised form 21 July 2008

increase the insight in dynamics and functioning of the complex CW system by using

Accepted 30 July 2008

mechanistic or process based models that describe transformation and degradation processes

Available online 7 September 2008

in detail. As these mechanistic models are complex and therefore rather difficult to use there is,on the other hand, a need for simplified models for CW design. The design models should be

Keywords:

premium to the currently used design guidelines that are mainly based on rules of thumb or

Constructed wetlands

simple first-order decay models. This paper presents an overview of the current developments

Design models

on modelling of subsurface flow CWs based on the modelling work and model developments

Mechanistic models

presented at the WETPOL 2007 symposium. Three kinds of models have been presented:

Numerical modelling

simple transport and first-order decay models, complex mechanistic models, and a simplified

Subsurface flow

model that has been developed for design of CWs. The models are presented and selected
results are shown and discussed in relation to the available literature.
© 2008 Elsevier B.V. All rights reserved.

1.

Introduction

During the last years an increasing interest in numericalmodelling of subsurface flow constructed wetlands (CWs) has
been observed. CWs are engineered systems designed to
optimize the treatment conditions found in natural environments; consequently, CWs are complex systems that are
difficult to understand. A large number of physical, chemical,

and biological processes are active in parallel and mutually
influence each other. Therefore, for a long time CWs havebeen often considered as “black boxes” and only little effort
has been made to understand the main processes leading towastewater purification. For the same reason, almost all theavailable design guidelines are based on empirical rules of
thumb, such as those using the specific surface area requirements (e.g. Brix and Johansen, 2004; ÖNORM B 2505, 2005;⁎ Corresponding author.
E-mail address: guenter.langergraber@boku.ac.at (G. Langergraber).
0048-9697/$ – see front matter © 2008 Elsevier B.V. All rights reserved.
doi:10.1016/j.scitotenv.2008.07.057

3932

S CI EN CE OF T H E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

DWA-A 262, 2006) or simple first-order decay models (e.g.
Kadlec and Knight, 1996; Rousseau et al., 2004). The main
objective of numerical modelling is to obtain a better understanding of the processes governing the biological and chemical transformation and degradation processes occurring in
CWs, thus providing insights into the “black box” Once reliable
numerical models are developed and validated against experimental data, they can be used for evaluating and improving existing design criteria.
Two main types of CWs are distinguished: surface flow andsubsurface flow CWs (e.g. Kadlec and Knight, 1996; Kadlec
et al., 2000; Haberl et al., 2003). Surface flow (SF) CWs aredensely vegetated systems with open water surfaces andtypically have water depths of less than 0.4 m. In subsurface
flow (SSF) CWs, no free water level is visible. SSF CWs arefurther subdivided into horizontal flow (HF) and vertical flow(VF) systems depending on the direction of water flow throughthe porous medium (sand or gravel). To prevent clogging of theporous filter material, the use of traditional SSF CWs is limited
to mechanically pre-treated wastewater, having low content
of particulate matter. Compared to SF systems, the contact
area of water with bacteria and substrate is much larger,
which decreases the area requirement of SSF CWs (Vymazal
et al., 1998; Haberl et al., 2003) while increasing the removal
efficiency of the system. This work focuses on model developments for SSF CWs. Langergraber (2008) reported that todate only few numerical models are available and able todescribe treatment processes in SSF CWs. Most of the literature on “models” for CWs is relevant to HF CWs and refers
to simple first-order decay models (e.g. Rousseau et al., 2004;
Stein et al., 2006) or still describes the CWs as a “black box” (e.g.
Pastor et al., 2003; Tomemko et al., 2007) thus providing only alimited understanding of the studied facility. The number of
mechanistic or process based models is instead very limited.
Such models can be further distinguished between models
describing only the hydraulic behaviour and single-solutetransport, models describing reactive transport in saturated
conditions such as in HF CWs, and reactive transport models
for variably-saturated conditions that have the highest complexity and that can be used for modelling VF CWs as well as
HF CWs. In this paper we summarize different contributions
presented at the WETPOL 2007 symposium relevant to SSF CW
model development and validation. We try to show some of
the most recent findings in this area, and identify some of the
aspects where further research is needed. Three different
kinds of models are presented in the subsequent discussion:
simple transport and first-order decay models, complex mechanistic or process based models, and a simplified model thathas been developed for design of CWs. After presenting thebasis of the models, results of the modelling work arepresented, summarized and discussed.

laws governing water flow and solute transport (such as themass balance and the advection–dispersion equation). By
defining appropriate boundary conditions and making additional assumptions on system geometry and behaviour, thegoverning equations are simplified and reduced to simple
closed-form formulas. These models have the obvious advantage that no special mathematical or numerical treatment isrequired to solve the equations for the variable of interest.
However, the important drawback is that the validity of theresults is limited to situations where the underlying assumptions hold.
While different works pointed out that such models are
generally not suitable to model reactive and degradation
processes in wetlands (e.g. Kadlec, 2000, 2003), simplified
models are often useful tools for evaluating water flow andtransport of inert solutes. Several of such methods exist andare routinely applied by practitioners to estimate residencetime, degree of mixing and hydrodynamic dispersivity. Amongthese methods, some of the best known are (i) the method of
moments, (ii) the dispersed plug-flow (DPF), (iii) the tanks-inseries (TIS) model and (iiii) the detention time gamma distribution (DTGD). Details of each method, strength and limitations can be found in the literature (e.g. Kadlec and Knight,
1996; Levenspiel, 1999; Kadlec, 2003). However, we would liketo point out that all the approaches listed assume that the flow
field can be somehow simplified to one-dimensional behaviour, and the effects of dispersion, heterogeneity and deadzones are lumped together in the hydrodynamic dispersivity
term. This may introduce a significant error in the computation e.g. of the effective residence time, which also can havetremendous impact on the estimation of the degradation
efficiency (Kadlec and Knight, 1996). Recovery rate can besometimes a useful indicator of the effectiveness of the applied model; however, this parameter alone does not ensurethat modelling results are correct.
Sandoval-Cobo and Peña (submitted for publication), Giraldi et al. (in press) and Mena et al. (2007) evaluated thereliability of these methods on both horizontal and verticalflow conditions, and compared the results obtained using different approaches on the same data set. Moreover, SandovalCobo and Peña (submitted for publication) in parallel with thetracer experiment measured COD concentration profiles
across the length of their pilot-scale HF CW, and applied
some simplified first-order degradation equations to fit theirexperimental data. The closed-form equations used by
Sandoval-Cobo and Peña (submitted for publication) are thesame as those used for non-reactive transport, with anadditional sink term to account for degradation and transformation. Additionally, the modified k–C⁎ model was also tested(e.g. Kadlec and Knight, 1996; Persson et al., 1999).

2.2.

2.

Model description

2.1.

Simple transport and first-order decay models

The first class of models we present in this paper consists ofanalytical, closed-form equations. The derivation of all themodels belonging to this class is based on the fundamental

Mechanistic or process based models

The multi-component reactive transport module CW2D (Langergraber, 2001; Langergraber and Šimůnek, 2005) was used byseveral authors (Korkusuz et al., 2007; Langergraber, 2007b;
Mena et al., 2007; Toscano et al., in press). CW2D was developed to model transport and reactions of the main constituents of municipal wastewater in subsurface flow constructedwetlands and is able to describe the biochemical elimination

S CI EN C E OF TH E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

and the transformation processes for organic matter, nitrogenand phosphorus. CW2D was incorporated into the HYDRUS
variably-saturated water flow and solute transport program(Langergraber and Šimůnek, 2006; Šimůnek et al., 2006a,b). The
HYDRUS program numerically solves the Richards equation
for saturated/unsaturated water flow and the convection–
dispersion equation for heat and solute transport. The flow
equation incorporates a sink term to account for water uptake
by plant roots. The solute transport equations consider
convective–dispersive transport in the liquid phase, diffusionin the gaseous phase, as well as non-linear non-equilibrium
reactions between the solid and liquid phases (Šimůnek et al.,2006b). The CW2D module considers 12 components and 9
processes. The components include dissolved oxygen, organic
matter (three fractions of different degradability, i.e. readilyand slowly-biodegradable, and inert), ammonium, nitrite,
nitrate, and nitrogen gas, inorganic phosphorus, and heterotrophic and two species of autotrophic micro-organisms.
Organic nitrogen and organic phosphorus are modelled as
nutrient contents of the organic matter, i.e. they are calculated
as a percentage of COD. The biochemical elimination andtransformation processes are based on Monod-type expressions used to describe the process rates. All process rates anddiffusion coefficients are temperature dependent. The processes considered are hydrolysis, mineralization of organic
matter, nitrification (modelled as a two-step process), denitrification, and a lysis process (as the sum of all decay and loss
processes) for the micro-organisms. CW2D assumes a constant concentration of micro-organisms (and other compounds) in each finite element. The thickness of the biofilm
is not considered (Langergraber, 2001). The mathematicalformulation of CW2D is based on the mathematical formulation of the Activated Sludge Models (ASMs, Henze et al., 2000).
Brovelli et al. (2007) implemented a set of biological andgeochemical reactions into the three-dimensional numericalsimulator PHWAT (Mao et al., 2006). PHWAT has been developed for reactive transport in porous media and consists ofthree different modules (flow, transport and biogeochemistry)
coupled with an operator-splitting technique. The flow module is based on MODFLOW (MacDonald and Harbaugh, 1988),
which uses a finite difference scheme to solve the saturated
water flow equation within the domain. The transport moduleis based on MT3DMS (Zheng and Wang, 1999). Several numerical solvers are available, thus making the code robust andflexible. The biogeochemistry module is based on PHREEQC-2(Parkhurst and Appelo, 1999), a general purpose model for
aqueous geochemistry. PHREEQC is extremely powerful andflexible, thus making PHWAT potentially able to handle all thereactions of interest for CW modelling. The biogeochemical
reaction network is similar to that of CW2D, based on the ASM
structure (Henze et al., 2000). Kinetic oxidation of carbon
sources, organic matter hydrolysis, nutrient transformations
and assimilation (nitrogen and phosphorous primarily) aremodelled via Monod-type kinetic equations. Oxygen dissolution is accounted for via a first-order mass-transfer equation.
Temperature dependence is included via the Arrhenius
equation (e.g. Langmuir, 1997). While the ASM model consistsof kinetic equations only, full water chemistry and sediment–water interactions can be modelled with PHREEQC and thus
PHWAT. This allows modelling of the pH variations that in-

3933

fluence biodegradation kinetics (Kadlec and Knight, 1996), as
well as redox and surface complexation reactions. An important advantage of PHWAT is that, thanks to its modular
structure, it can be extended to account for the coupling between flow and transport, e.g., pore clogging due to suspended
solid deposition or biomass growth. The model has beendeveloped very recently, and the validation of PHWAT withexperimental data is currently on-going. Some preliminary
results have been presented (Brovelli et al., 2007), showing
reasonably good agreement with experimental data. However,model calibration has been found extremely difficult, due tothe strong non-linearities of the mathematical model.
FITOVERT (Giraldi et al., 2008) was developed specifically tosimulate VF CWs and considers therefore only a one-dimensional flow from the top to the bottom of the system. The
dynamic formulation of the model allows simulating thetypical non-stationary behaviour of VF CWs. Boundary conditions are provided to simulate both variable inflows andponding at the surface and free drainage, drainage with apressure head and no drainage at the bottom of the system.
The hydraulic flow for unsaturated condition is described bymeans of the Richards equation, while the constitutive
relationships among pressure head, hydraulic conductivity
and water content are described using the van Genuchten
parameterization. Evaporation at the CW surface and transpiration through plants are also considered. FITOVERT
describes biochemical processes of organic matter and nitrogen by means of the standard ASM (Henze et al., 2000).
Advective and diffusive transport of contaminants in theliquid phase is implemented according to a mass conservation
equation with a limited change to the dispersion term. Gas
transport and transfer to the liquid phase are described for theoxygen by a mass balance with advection and dispersion.
Settling and filtration of particulate matter is included thus
allowing porosity reduction due e.g. to bacteria growth and todeposition of particulate components, so that the clogging
process could be also simulated; the effect of pore size reduction on the saturated hydraulic conductivity is also considered. Giraldi et al. (in press) used the DPF model implemented
in FITOVERT.

2.3.

Simplification of complex models

Meyer et al. (2007) developed a simplified reaction model for
CWs treating combined sewer overflow (CSO), also known asretention soil filters (RSFs), based on experiences with CW2D
(Langergraber and Šimůnek, 2005). The background for thiswork was that the German design- and operation-guideline for
RSFs (DWA-M 178, 2005) introduced the definition of theannual hydraulic load (in m3/m2) by calculating CSO volumes
as result of long-term simulations of combined sewer-systems
as a main criteria for filter design. Therefore the objective was
the development of a tool that allows long-term simulations
and can be combined with pollution-load-models for sewersystems to optimize RSF design and operation.
Based on their experience Meyer et al. (2007) concluded
that modelling of CSO treatment in CWs to optimize designand operation requires a reasonable balance between detailed
description and practicable handling. The use of complex
models like CW2D has been shown to be not practicable for

3934

S CI EN CE OF T H E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

long-term simulations in combination with sewer-systemmodels (Dittmer et al., 2005; Meyer et al., 2006; Henrichs
et al., 2007). Based on this experience and further experiencefrom long-term simulations with KOSMO (Kaufmann andSchmitt, 2005) a RSF module for KOSMO has been developed.
The module RSF_Sim (Retention Soil Filter Simulation) is notdirectly embedded into the KOSMO simulation program.
RSF_Sim input data either can be generated manually usinga graphical user-interface or can be imported directly fromKOSMO.
RSF_Sim is able to describe filter layers whereby water flow
is described using a sequence of stirred reactors with variable
water contents (Fig. 1). At the beginning of loading events theretention layer (air) on top is filled as long as the loading rate ishigher than the infiltration rate into the filter body. The RSF
filter body consists of cover layer (gravel), filter layer (sand)and drainage layer (gravel). The infiltration rate into the RSF isdetermined by the hydraulic conductivity of the sand. The
outflow of the drainage layer is limited to a maximum flowrate (lower than the hydraulic conductivity of the filter layer).
In the case of very big rain events the rising water level in thesystem (dotted line) can lead to overflow. As a special feature
the filter layer itself is divided into two layers thus allowing todescribe different processes at different depths of the filter
separately.
Meyer et al. (2007) assumed that all processes occur only inthe sand layers, i.e. filter layers 1 and 2. Processes consideredare:“F” — Filtration: filtration of large particles on layer surface(filter layer 1) and of small particles in layer space (filter
layer 1 + 2),“S” — Sorption described by Freundlich-isotherms,“D” — Degradation described by a decay rate.
Combinations of the processes are possible (e.g. “FD” and“SD”). Additionally, the process “SP” describes degradation ofammonium nitrogen only in dry periods between loadings(dry periods between loading events can be extensive, e.g. Uhl
and Dittmer, 2005). Filter sand performance can be increased

Fig. 1 – RSFP Sim hydraulics as a sequence of stirred reactors
with variable water contents (Meyer et al., 2007).

Fig. 2 – Characteristic TIS model fitting to tracer data
(Sandoval-Cobo and Peña, submitted for publication).

or decreased for each pollutant separately, and all degradation
processes are temperature dependent.
According to Meyer et al. (2007) first simulation runs with
RSF_Sim gave plausible results. However, no calibration andvalidation of the model has been carried out up to now. As anext step Meyer et al. (2007) plan to calibrate the RSF_Sim
module parameters for single events and long-term scenarios.
Measured inflow- and outflow-curves shall be used forvalidation with manual data input before using RSF_Sim in
combination with KOSMO as a design tool.

3.

Results

3.1.

Simple flow and transport models

3.1.1.

Modelling of HF CWs

Sandoval-Cobo and Peña (submitted for publication) presented a study on HF CWs located in Columbia. The pilotscale HF CW (9 m length, 3 m width and 0.6 m depth) received
the effluent of an anaerobic pond treating domestic sewage
and was planted with Phragmites australis. The unit was
designed to receive a hydraulic loading rate (HLR) of 0.19 m
d− 1 of domestic wastewater with mean influent concentrations of 136 mg COD L− 1, 116 mg BOD5 L− 1, 45 mg NH4–N L− 1
and 5 mg PO4–P L− 1. Tracer experiments with an organic
fluorescent tracer Rhodamine WT (RWT) were performed todescribe the hydraulic behaviour of the system.
For wastewater characterisation samples were analyzed for
COD, BOD5, NH4–N, NO3–N, NO2–N, PO4–P, pH and temperature.
In order to evaluate longitudinal COD concentration profiles
and to obtain removal rate constants (k) across the wetland
length, samples were taken from four sampling locations along
the experimental unit centre line (at the influent, at 1/3 and 2/3of the length, and at the effluent).
Six tracer experiments were conducted, and results wereanalyzed using three simplified models. One of the characteristic residence time distributions (RTDs) and the corresponding fitted curve for the TIS model is shown in Fig. 2. Basically,
all six outflow tracer response curves produced skewed bellshaped RTDs with long tails, as observed in previous tracer
studies carried out on HF CWs (Chazarenc et al., 2003; MarsilliLibelli and Checchi, 2005).

S CI EN C E OF TH E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

3935

Table 1 – Hydraulic parameters for different models(SD = standard deviation).
Tracerrunt0(h)Method ofmoments
τ(h)

R (%)
a

E1
E2
E3
E4
E5
E6

35.6
35.6
33.8
33.9
33.7
33.5

41.3
42.9
44.0
49.2
48.3
53.9

N.A
61
87
91
81
75

Average
SD

34.4
1.0

46.6
4.7

–
–

DPF model

TIS model

δN

τTIS (h)

0.07
0.05
0.11
0.11
0.10
0.16

8
11
6
6
6
4

37.4
39.6
38.4
41.2
40.5
43.3

82
73
84
91
76
75

0.10
0.04

7
2

40.1
2.1

–
–

R (%)

a

R = The whole mass recovery of tracer from each test (the value N.
A accounts for unrealistic tracer recovery estimation due toinstabilities experienced with the LIF meter used to determine
Rd-WT concentrations in the first experiment).

Table 1 shows the hydraulic parameters calculated for eachof the tracer experiments using different models: Method of
moments, DPF and TIS models. The average real hydraulicresidence times (HRTs) obtained were above the theoretical
HRT value (t0) for all models showing the inadequacy of
hydraulic behaviour expected from the plug-flow basedmodels. This can be explained by the fact that homogeneous
velocity profiles (τ = t0) in HF CWs never occurs. Key factors
such as porosity change due to biofilm accumulation, plantcolonization and evapotranspiration (ET) rates must be taken
into account. For example, plant density increased from 6 to
150 plants m− 2 during the 5 months of experimental work
while plant roots occupied the whole wetland depth. Average
ET in the planted unit was calculated from water massbalance. Sandoval-Cobo and Peña (submitted for publication)found that about 29% of the water was lost due to ET, thus
dramatically changing the hydrodynamics of the system.
These two aspects must be considered when designing HF
CWs under particular environmental conditions.
First-order degradation models (PF, modified k–C⁎, DPF and
TIS models) were fitted to data of COD concentrations versus
the fractional distance (Table 2 and Fig. 3). The parametersdescribing the hydraulic behaviour, the dispersion number δ
for the DPF model, and the number of tanks (N) for the TIS
model have been taken from Table 1. Values for COD removal
rate constants (Table 2) were obtained by applying each model
to average COD concentration profiles measured. The optimization goal was firstly to fit the effluent concentration andsecondly the concentration profile.
A good match to the experimental data could be obtained
only for the k–C⁎ model. Plug-flow based models, i.e. PF, DPF

Table 2 – Parameters for first-order decay models obtained
for the data from Sandoval-Cobo and Peña (submitted for
publication).
Model
Parameter
Value

k–C⁎
−1

k (d )
4.21

PF−1C⁎ (mg L )55−1

k (d )
0.88

Influent concentration C0 = 136.5 mg COD L− 1.

DPF−1k (d )
0.90

TIS
k (d− 1)
0.96

Fig. 3 – Comparison of measured (averages and standard
deviation, n = 6) and modeled COD concentrations.

and TIS models, k resulted in similar values when matching
the effluent concentrations but did not fit the COD profiles
measured. Sandoval-Cobo and Peña (submitted for publication) fitted all six measured concentration profiles separately
and found the calibration of the k–C⁎ model parameters lead to
obtaining a good correlation (r2 ≥ 0.98) for all six data sets.
Moreover, the average values of the k and C⁎ parameters were
between those reported for HF CWs (Rousseau et al., 2004).
The values of k for CWs vary widely when compared toother wastewater treatment systems (Stein et al., 2006). The
raw data series from which Table 2 mean figures were
calculated show that k values for all models varied consistently throughout the six experiments. This feature has strong
implications in the design of full-scale HF CW systems andreinforces the idea expressed by some authors about the convenience of using design parameters obtained from previously
tested pilot units or HF CW operating under similar conditions(Rousseau et al., 2004). The combined use of dispersion studies
and simplified mathematical models to evaluate the performance of organic matter removal in HF CW allowed obtaining
more realistic key design parameters. The consideration ofboth flow dynamics and first-order kinetics leads to moresensible rate constant values that may have strong implications for the design and evaluation of HF CWs treatment units
located in tropical regions.
The results obtained herein show that the k–C⁎ model
reproduced well the performance of HF CWs in tropical regions,thus, the model may be used for design purposes provided that
overall rate constants and characteristic hydraulic parameters
are properly derived from experimental data.

3.1.2.

Modelling of VF CWs

Giraldi et al. (in press) analyzed a VF bed that was part of a CW
pilot plant operating as a tertiary treatment system after the
Soreq Mechanical Biological Wastewater Treatment Plant serving the metropolitan area of Tel Aviv (Israel). The plant has
been in operation since February 2005. The bed analyzed was a
33 m2 VF CW with a total depth of 68 cm and 6 gravel layers(listed from the top to the bottom): (1st) 3 cm depth, particle size:
20–30 mm, (2nd) 9 cm depth, particle size: 2–3 mm, (3rd) 8 cm
depth, particle size: 5–10 mm, (4th) 8 cm depth, particle size:
8 mm, (5th) 15 cm depth, particle size: 3–6 mm, (6th) 20 cm
depth, particle size: 50–60 mm. The system was planted with

3936

S CI EN CE OF T H E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

Fig. 4 – Experimental tracer breakthrough curves (dots) and numerical simulations (lines) for complete saturation (Test C, left)
and free drainage (Test E, right); Giraldi et al. (in press).
different plant species: Cyperus papyrus, Canna sp., Iris pseudoacorus, P. australis and Juncus ensifolius.
Five tracer tests were performed with RWT for different
hydraulic loads and saturation levels of the VF bed. The tracer
tests were carried out in steady state conditions, feeding thesystem with controlled tap water flow and dosing RWT as
square input signals. Each tracer test was temporally limited
within a single workday for technical and safety reasons. No
rain occurred during the tests.
Two experimental tracer breakthrough curves and therespective simulation results are reported in Fig. 4, while
complete results of the RTD analysis are reported in Table 3,along with the main experimental parameters of the tracer
tests. Background tracer concentrations were always negligible and the square input signals allowed obtaining welldefined RTD curves within single workday intervals.
The mean residence time was always lower than thenominal retention time (when calculated), so the presence ofdead zones in the bed could be highlighted. RTD analysis
indicates that the HRT and the degree of apparent mixing (i.e.
variance) increase as the water content in the system increases and as the water flow decreases. On the other handthe degree of local mixing (i.e. dimensionless variance)increases as the saturation degree of the bed decreases andslightly increases with the hydraulic load. The mono-dimensional plug-flow with dispersion model for unsaturated
conditions resulted very satisfying in simulating all the tracer
tests, both for saturated and unsaturated conditions. Giraldi
et al. (in press) found that the dispersivity (λ) decreases as thesaturation of the system increases, reflecting the behaviour ofthe dimensionless variance. The correlation between dispersivity and water content is discussed widely in literature
dealing with the mechanisms of water and solute transport in
the vadose zone. For a saturated system, λ is a function of thecharacteristics of the porous medium (Ogata, 1970). Straub andLynch (1982) reported that λ describes the variation of fluid
velocity from its average value and that λ increases with theeffects of channelling and non-uniformity of the porous medium. Therefore Giraldi et al. (in press) assumed a dependenceof λ on the degree of saturation of the system, since the watercontent affects the uniformity of the flow inside the porous
medium. In fact, for unsaturated conditions local mixing due
to different flow paths should increase, while saturation
conditions provide more uniform flow paths, thus reducing
local mixing. For their case study Giraldi et al. (in press) foundthe dependency between dispersivity and saturation degree tobe linearly decreasing, but more data on different unsaturated
conditions need to be analyzed to confirm this statement or tovalidate any other non-linear relationships. This is importantfor future research: numerical models describing the behaviour of VF CWs should implement this kind of dependency inorder to obtain meaningful simulations for the transport (andconsequently the fate) of solute compounds inside these kindsof systems. This is required to simulate both standard operation of VF CWs (periodic pulse loadings with draining conditions) and innovative operation schemes like the Tidal Flowmode (percolating with progressive saturation and sudden
outflow).

Table 3 – Main experimental parameters and results of RTD analysis (Giraldi et al., in press).
Test Saturation Hydraulic
Inlet tracer
load
concentration
Label
A
B
C
D
E

PartialComplete
Complete
Draining
Draining

DosingperiodMonitoring Tracer mass
periodrecoveryMeanresidencetime

Dimensionless Dispersivity
variance
(λ)

L/min

μg/L

min

h

%

min

–

cm

24
48
24
24
24

56
22
110
57
56

180
122
169
98
230

8
6
14
6
10

88
58
59
63
67

128
113
341
52
51

0.38
0.24
0.18
0.60
0.53

10
4.5
4.5
14
14

S CI EN C E OF TH E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

Fig. 5 – Comparison of experimental tracer concentrations tosimulation results for CW1 (planted with Phragmites australis)using different models (according to Mena et al., 2007).

3.2.

Simulations with CW2D

3.2.1.

Modelling of HF CWs

Mena et al. (2007) presented simulation results of tracer experiments and multi-component reactive transport obtained
for a pilot-scale experimental installation situated in southern
Spain. The installation consisted of three parallel operated HF
beds. Each bed consisted of a channel (size 2.5 m × 0.65 m witha depth of 0.6 m) with a longitudinal slope of 1%. Different
species of macrophytes were planted: in CW1 P. australis(common reed), in CW2 Lythrum salicaria (purple loosestrife)
whereas CW3 remained without plants. All beds were filledwith gravel with a diameter of 6–9 mm (porosity 0.38), apartfrom the first and last 10 cm, for which the diameter was
9–12 mm (porosity 0.40) to improve the distribution of the3937

wastewater in the inlet zone and to collect the treated wastewater in the outlet zone, respectively. The CWs were continuously fed with artificial urban wastewater. The mean
influent concentrations have been 127 mg COD L− 1 and 10 mg
NH4–N L− 1 (Mena et al., 2007; Mena, 2008).
Tracer experiments were carried out using bromide. One
litre solution of sodium bromide with a concentration of
5000 mg L− 1 was fed by means of an impulse in each CW. So,
the total tracer mass initially added was 5000 mg. After thepulse, wastewater feeding was kept running continuously andeffluent samples were collected. Fig. 5 shows the experimental
tracer concentration and the fit of the different hydraulicmodels for CW1. Mena et al. (2007) showed that for all experiments simulation with HYDRUS fitted better than DTGD
and DPF models. Especially the DPF model did not fit very wellto the experimental data mainly because the long skewed“tail” increases the mean HRT and this fact makes the DPF to
be delayed in relation to the experimental data.
The experiments by Mena (2008) for pollutant removal
have not been designed for modelling purposes. Thus, only
the mean effluent concentrations of the different parallel beds
could be modelled by introducing different oxygen release
rates for the different plant species. Oxygen release rates
caused by P. australis (CW1) and L. salicaria (CW2) have beenfound to be 2.06 and 5.92 g O2 m− 2 d− 1, respectively; oxygen
transport by diffusion into CW3 amounted to 5.38 g O2 m− 2 d− 1.
Differences between the efficiency of nutrient and organic
matter removal have been found not only because of thepresence of plants but also because of the plant species.
Mena et al. (2007) concluded that plant species influenceboth the hydraulics in and the removal efficiency of HF CWs.
The hydraulic behaviour is influenced due to plant growthFig. 6 – Simulated water flow (top) and tracer experiments (bottom) for Turkish (left) and Austrian (right) zeolite (Korkusuz et al.,2007).

3938

S CI EN CE OF T H E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

increases ET and decreases dispersion, respectively. Tracer
studies determined that the effective volume ratio can behigher than 1 caused by dead zones or hydraulic phenomena
in which the tracer is accumulated and is being released
gradually. The removal performance is influenced by plants asthey can increase dissolved oxygen concentration and take-up
nutrients but also because rhizospheric growth increases thesurface area for biofilm development.

3.2.2.

Modelling of VF CWs

Korkusuz et al. (2007) investigated 10 different natural andartificial substrates to be used for the main layer of VF CWs.
Simulation results have been presented for lab-scale studies
carried out at the technical lab of the Institute of SanitaryEngineering at BOKU (University of Natural Resources andApplied Life Sciences, Vienna). Each lab-scale column had adiameter of 20 cm and was filled with a 50 cm main layer using
different substrates. The substrates (gravel size: 0–4 mm)
tested were originally from Turkey (zeolite, sand, pumice,
perlite, blast furnace granulated slag) and Austria (zeolite,
sand, ferrosorp, crushed-concrete) and have been selected
because they have a potential higher nitrogen and phosphorous adsorption capacity. The following parameters have been
determined for all materials: particle and bulk density,
porosity, saturated hydraulic conductivity and grain size
distribution.
The lab-scale columns have been loaded four times a daywith mechanically pre-treated municipal wastewater. Influent and effluent flow rates of 10 lab-scale columns have beenmeasured. Tracer studies (with KCl dissolved in tap water)
have been conducted to describe single-solute transportcharacteristics. Electrical conductivity of influents and effluents have been measured online and recorded with a datalogger.
Fig. 6 shows simulation results for effluent flow rates froma single loading and from tracer experiments for twomaterials: Turkish and Austrian zeolite. The simulations
have been carried out as follows: Measured values of bulk
density, porosity and saturated hydraulic conductivity havebeen used as input parameters. Water flow has beencalibrated using effluent flow rate measurements with theinverse simulation procedures implemented in HYDRUS(Šimůnek et al., 2006b). The parameters obtained for the flow
model have been used for single-solute transport modelling of
the tracer experiments. The simulations matched the measured data well for all columns. Korkusuz et al. (2007) pointed
out that one main result of the study was the generation of
input data sets for CW2D for different substrates. Further
studies on reactive transport simulations regarding theremoval efficiency are planned.
Langergraber (2007b) presented results from a simulation
study using data from an experimental site located at thewastewater treatment plant Ernsthofen in Lower Austria (e.g.
Langergraber, 2007a). A two-stage CW system that consists oftwo vertical flow beds with intermittent loading operated inseries was investigated (Langergraber et al., in press). The firststage uses a grain size of 2–3.2 mm for the main layer and has adrainage layer that is impounded; the second stage a grain sizeof 0.06–4 mm and a drainage layer with free drainage. The twostage system was operated with an organic load of 80 g COD

m− 2 d− 1 for the first stage (1 m2/person equivalent), i.e. 40 g
COD m− 2 d− 1 for the whole system (2 m2/person equivalent).
The filter beds are planted with common reed (P. australis).
Langergraber et al. (in press) showed that the two-stage CW
system performed better compared to a single-stage VF bed
designed and operated according to the Austrian design
standards (ÖNORM B 2505, 2005) with 4 m2/person equivalent.
The better performance resulted in lower effluent concentrations for low temperatures (below 6 °C) and an increased
nitrogen removal of 53% at high average nitrogen elimination
rates of about 1000 g N m− 2 yr− 1, respectively.
Simulation results by Langergraber (2007b) showed that amore stable performance of two-stage CW system (compared
to a single-stage VF bed) can be also predicted when the experimental plant is operated according to the testing procedures required for small wastewater treatment plants (ÖNORM
EN 12556-3, 2001). Fig. 7 shows the simulated NH4–N effluent
concentrations of stage 1 and 2 of the two-stage CW systemand of the single-stage CW system for an experiment with
2 days overload (2 days with 150% load within a 14 day period).
The influent concentrations for the simulations have been
470 mg COD L− 1 and 55 mg NH4–N L− 1. The simulated NH4–N
effluent concentration of the two-stage CW system did notexceed 6 mg NH4–N L− 1 whereas effluent concentration of thesingle-stage VF bed exceeded 15 mg NH4–N L− 1 during theoverload period. The results of these simulations still have tobe verified with experiments at the experimental plant which
are planned for 2008.

3.2.3.

Modelling of hybrid systems

Toscano et al. (in press) presented simulation results for thehydraulic behaviour and effluent pollutant concentrations in apilot-scale two-stage SSF CW for treatment of municipal
wastewater. The experimental plant is located in San Michele
di Ganzaria (Eastern Sicily). It was used for secondary or tertiary
treatment and consists of four parallel lines each with two SSF
beds in series. Each bed has a rectangular shape with a surface
area of 4.5 m2 (1.5 m × 3.0 m). For their study Toscano et al. (inpress) used data from line 1 (planted with Phragmite sp.) and line
2 (unplanted), each comprising a HF and a VF bed in series. The
HF beds were filled, to an average depth of 0.6 m, with volcanic
gravel having a size of about 10–15 mm. The water depth in HF
systems is about 0.55 m. The VF beds were filled with two
different medium layers: volcanic sand as main layer for a depth

Fig. 7 – Simulated NH4–N effluent concentrations for a periodof 2 days overload.

S CI EN C E OF TH E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

3939

in the effluent of the HF bed for both lines (Toscano et al., inpress).
Toscano et al. (in press) concluded that the capability of
HYDRUS/CW2D to model water flow and tracer experiments inboth horizontal and vertical subsurface flow beds as well as tomodel the removal of some solute wastewater compounds
was proved for hybrid systems.

4.

Fig. 8 – Measured and simulated data for a tracer experimentof line 2 (Toscano et al., in press).

of 0.5 m (about 0.06–4 mm; d10 = 0.13 mm; d60 = 2.0 mm) at the topof the filter bed and coarse volcanic gravel as drainage 472 layer(about 16–32 mm, 25 cm depth). The average influent flow rate
for each line was about 40± 5 L h− 1.
Effluent flow rates measurements have been used tocalibrate the flow model. Fig. 8 shows measured and simulated data for a tracer experiment of line 2 (unplanted). The
simulated breakthrough curves matched the measured datawell. For reactive transport the parameter set as proposed byLangergraber (2001) was used. Mean values for organic matterand nitrogen parameters could be simulated for both lines byconsidering nutrient uptake by plants and oxygen release by
the plant roots. The influent concentrations for the simulations have been 23 mg COD L− 1 and 22.1 mg NH4–N L− 1 fortertiary treatment and 64 mg COD L− 1 and 8.6 mg NH4–N L− 1 forsecondary treatment, respectively.
Since no significant changes of the influent concentrations
over time occurred, a further experiment was conductedsimulating influent ammonia variations and checking theanswer of the model. This was done by adding ammonium
chloride (NH4Cl) as pulse tracer in the influent of lines 1 and 2
functioning as secondary treatment. In Fig. 9 measured andsimulated NH4–N concentrations are shown. The values forthe parameters describing plant uptake and oxygen releaseobtained in the previous simulations have been used for thesimulations. CW2D seems to be able to predict the dynamics

Discussion

Most of the currently available literature aimed at modelling of
CWs makes use of simplified equations. The underlying assumptions of such models are often that water flow can beconsidered homogeneous and mono-dimensional, or evensimplified to plug-flow with or without hydrodynamic dispersion. Tracer experiments for HF CWs have been carried out by
Sandoval-Cobo and Peña (submitted for publication) and Menaet al. (2007) and have been modelled using tanks-in-series(TIS), detention time gamma distribution (DTGD) and dispersed plug-flow (DPF) models. Mena et al. (2007) showed thatfor all their tracer experiments the simulations with HYDRUS(Šimůnek et al., 2006a) matched the data best and there was apoor match for DTGD and TIS models thus indicating that amechanistic model such as HYDRUS also performs best for
single-solute transport simulations in HF CWs.
Giraldi et al. (in press) found that mono-dimensional plugflow with dispersion can be suitable to simulate bothsaturated and unsaturated conditions in VF CWs. The hydrodynamic analysis of Giraldi et al. (in press) also pointed out thecorrelation between dispersivity and water content for VF
CWs; the dependency between dispersivity and saturation
degree seems to be linearly decreasing. The same authors
pointed out the need to introduce this dependency in the flow
model for a successful modelling of VF CWs.
All hydraulic retention times obtained from experimental
data have been found to be different than the calculated
nominal retention times (Giraldi et al., in press; SandovalCobo and Peña, submitted for publication; Mena et al., 2007).
This can be explained because by the fact that homogeneous
velocity profiles in HF CWs that are an underlying assumptionof the simplified models never occur. Key factors such as thepresence of dead zones, porosity change due to biofilm

Fig. 9 – Measured influent concentrations and measured and simulated effluent concentrations of NH4–N after the addition of
NH4Cl solution in lines 1 (left, planted) and 2 (right, unplanted) (Toscano et al., in press).

3940

S CI EN CE OF T H E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

accumulation, plant colonization and evapotranspiration are
not taken into account.
As for the biological cycling of nutrients and organic carbon, this models use simple first-order decay laws. SandovalCobo and Peña (submitted for publication) tested some ofthese methods on HF CWs. For pollutant removal using firstorder decay models Sandoval-Cobo and Peña (submitted for
publication) concluded that k–C⁎ models reproduced well theperformance of their HF CW and that k–C⁎ models may be used
for design purposes provided that overall rate constants andcharacteristic hydraulic parameters are properly derived fromexperimental data. It has to be mentioned that the adequacyof the k–C⁎ model to describe pollutant removal in HF CWs
may lay not on the typical background concentration interpretation but rather in the saturation-like kinetics shown byC⁎ value.
The k–C⁎ first-order degradation model is mostly applied to
HF CWs to describe removal of organic matter (e.g. Rousseau
et al., 2004; Mena, 2008). Other applications of the k–C⁎ model
found in literature are e.g. batch fed systems for which Steinet al. (2006) estimated the temperature dependence of themodel parameters of the k–C⁎ model for COD removal in abatch fed CW system for four plant species and Chan et al.
(2008) modelled ammonia removal using the k–C⁎ model.
Different design models for HF CWs have been comparedby Rousseau et al. (2004). By showing the output variability ofthese models it has been shown that the k–C⁎ model is themost suitable for design purposes. Already Kadlec (2000) concluded that first-order models are inadequate for the design oftreatment wetlands. However, Rousseau et al. (2004) claimed
that k–C⁎ model can be used for design only true when model
parameters from systems operating under similar conditions(climatic conditions, wastewater composition, porous filter
material and plant species) are available.
Langergraber (2008) reviewed mechanistic models available to simulate removal processes in SSF CWs. PHWAT(Brovelli et al., 2007) belongs to the category of models that can
only describe saturated water flow and are therefore only
applicable to HF CWs. This is due to the fact that the underlying model for water flow is MODFLOW which can only be
used for saturated porous media. However, by using MODFLOW as a sub-model PHWAT uses the most advanced waterflow model in this category as other models that are able todescribe both carbon and nitrogen transformations undersaturated water flow conditions use a tanks-in-series approach (Wynn and Liehr, 2001; Rousseau, 2005).
The model developed by Rousseau (2005) uses a reaction
model that is presented in matrix notation based on themathematical structure of the ASMs similar to the PHWAT
reaction model structure. Rousseau's model describes aerobic,
anoxic and anaerobic processes and has been applied to pilotscale HF CWs to evaluate porosity decreases in experimental
HF CWs submitted to different pretreatments (García et al.,2007).
Besides the multi-component reactive transport module
CW2D (Langergraber, 2001) four models could be been found inliterature that consider vadoze zone processes and are therefore suitable to simulate both HF and VF CWs. Langergraber(2008) reported that all these models are in a very early stage ofdevelopment. Different approaches are used for modellingvariably-saturated water flow, two models, Ojeda et al. (2006)
and Wanko et al. (2006), using the Richards equation andtherefore a mechanistic approach such as HYDRUS (Šimůnek
et al., 2006a), whereas Freire et al. (2006) use a combination of
CSTRs and dead zones and McGechan et al. (2005a) different
layers through which water flows vertically trough the CW.
CW2D is able to describe the biochemical transformation anddegradation processes for organic matter, nitrogen and phosphorus in SSF CWs whereas Ojeda et al. (2006) describe processes affecting solids, organic matter, nitrogen and sulphur
and Wanko et al. (2006) describe organic matter removal andoxygen transport only.
Recently a new mechanistic model (FITOVERT) was presented (Giraldi et al., 2008) that was developed to simulate VF
CWs. The hydraulic sub-model of FITOVERT also uses theRichards equation to describe water flow in one-dimensionalvertical layers. FITOVERT is able to simulate the fate of organic
matter, nitrogen, particulates and oxygen transfer. First
hydrodynamic simulations proved to be satisfying (Giraldi
et al., in press) and the calibration of the biochemical parameters has already started.
Results using CW2D/HYDRUS (Langergraber and Šimůnek,2005) showed that water flow, single-solute transport andreactive transport simulations resulted in good agreementswith measured data for HF and VF CWs as well as for hybrid
systems (Korkusuz et al., 2007; Mena et al., 2007; Toscano et al.,in press). The practical applications have shown that simulation results match the measured data when the hydraulicbehaviour of the system can be described well. A good matchof experimental data to reactive transport simulations canthen be obtained for CWs treating municipal wastewater using
the values for the CW2D model parameters as given by Langergraber and Šimůnek (2005). Therefore it is advisable tomeasure at least the porosity and saturated hydraulic conductivity of the filter material to obtain reasonable simulation
results for water flow (Langergraber, 2008).
For HF CWs the influence of plant species on the removal
efficiency was shown by Mena et al. (2007) and Toscano et al.
(in press). The removal performance has been influenced byplants as plants can increase dissolved oxygen concentrationand take-up nutrients. Differences between the efficiency ofnutrient and organic matter removal have been found not only
because of the presence of plants (Toscano et al., in press) but
also related to different plant species (Mena et al., 2007).
Meyer et al. (2007) used their experience from detailed
simulations with CW2D (Dittmer et al., 2005; Meyer et al., 2006)
to develop a simplified reaction model for CWs treating combined sewer overflow (CSO). The simplified reaction model
models the hydraulics as a sequence of stirred reactors withvariable water contents and considers only the most relevant
processes (filtration, sorption and degradation). The differences between CWs for treating municipal wastewater andCSO have also been pointed out by Henrichs et al. (2007). They
also used CW2D for simulating CWs for CSO treatment andfound limitations such as that the degradation processes indry periods should be different than those modelled in CW2D
for organic matter degradation and nitrification. Theselimitations can be overcome by the RSF_Sim module asproposed by Meyer et al. (2007) that is planned to be used asa design tool.

S CI EN C E OF TH E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

To use mechanistic models to gain insights into complex
systems such as “constructed wetlands” and to develop simplified design models based on these numerical simulations
could be a way toward bringing more reliability into the design
rules and to make modelling tools more acceptable for design
purposes in practice. These simplified design models should
comprise only the most relevant processes for the specificapplication and have therefore only few parameters. This can beof great importance for a better acceptance of models in thedesign of CWs. Even if a user-friendly operating environment isprovided for mechanistic models making their use rather simple
there are a lot of parameters to consider and it takes quite sometime to gain experience with the using of these models.

5.

Conclusion

Current developments in numerical modelling of CWs show
two well-differentiated but related objectives. Firstly, mechanistic models aimed at gaining insight into wetlands dynamics
and functioning, and secondly, simplified but robust andreliable models for design purposes. Both objectives needscientifically sound approaches and most likely the latter are anatural result of the former. Thus, joint work and proper
communication between these two areas of research have tobe pursued.

Acknowledgements
This paper summarizes the findings on numerical modellingof subsurface flow CWs that have been presented during the
2nd International Symposium on “Wetland Pollutant
Dynamics and Control WETPOL 2007” held from 16–20 September 2007 in Tartu, Estonia. The authors of the paper havebeen present and presented their work at WETPOL 2007 andespecially thank their co-workers for their contributions.
The work of David Giraldi was funded by the Italian Ministryof the Environment and Territory within Project 6 of the Italian–
Israeli Cooperation on Environmental Technologies.
The development of long-term simulation models for RSFs
by Daniel Meyer is funded by a PhD-scholarship of the German
Friedrich-Ebert-Stiftung.
Miguel Peña especially acknowledges Mr. J.J. SandovalCobo, M.Sc. student at the Universidad del Valle, Cali,
Colombia. Their research was funded by the Colombian
Institute of Science and Technology (COLCIENCIAS) under
the research contract No 1106-12-13835.
Alessandro Brovelli acknowledges his colleagues Serge
Baechler, Luca Rossi and D. Andrew Barry from EPFL, Lausanne,
Switzerland.
The research of Attilio Toscano was funded, within the
AQUATEC Project (PON 2000–2006), by the Italian Ministry ofInstruction, University and Research (MIUR).
The research of Ms. E. Asuman Korkusuz was funded by theEuropean Commission “Marie Curie Intra European Individual
Fellowship” project ONUREM (Project No. 515515: Optimization of Nutrient Removal in Constructed Wetlands Using
Special Substrates and Numerical Simulation, duration 09/
2005–08/2007).

3941

The experiments for the work presented by Günter Langergraber have been carried out within the project “Optimizationof subsurface vertical flow constructed wetlands” funded by
the Austrian Ministry for Agriculture, Forestry, Environment
and Water Management.

REFERENCES
Brix H, Johansen NH. Guidelines for vertical flow constructedwetland systems up to 30 PE. Økologisk Byfornyelse og
Spildevandsrensning No.52. Copenhagen, Denmark:Miljøstyrelsen, Miljøministeriet; 2004.
Brovelli A, Baechler S, Rossi L, Langergraber G, Barry DA. Coupledflow and hydro-geochemical modelling for design andoptimization of horizontal flow constructed wetlands. In:
Mander Ü, Kóiv M, Vohla C, editors. 2nd International
Symposium on “Wetland pollutant dynamics and control
WETPOL 2007” — extended abstracts, 16–20 September 2007,
Tartu, Estonia, vol. II. ; 2007. p. 393–5.
Chan SY, Tsang YF, Cui LH, Chua H. Domestic wastewater
treatment using batch-fed constructed wetland and predictivemodel development for NH3–N removal. Process Biochem
2008;43(3):297–305.
Chazarenc F, Merlin G, Gonthier Y. Hydrodynamics of horizontal
subsurface flow constructed wetlands. Ecol Eng
2003;21(2–3):165–73.
Dittmer U, Meyer D, Langergraber G. (2005): Simulation of a
subsurface vertical flow constructed wetland for CSO
treatment. Water Sci Technol 2005;51(9):225–32.
DWA-A 262. Grundsätze für Bemessung, Bau und Betrieb vonPflanzenkläranlagen mit bepflanzten Bodenfiltern zur
biologischen Reinigung kommunalen Abwassers. Hennef,
Germany: DWA — Deutsche Vereinigung für Wasserwirtschaft,
Abwasser und Abfall e.V., 2006. [in German].
DWA-M 178. Empfehlungen für Planung, Konstruktion und Betrieb
von Retentionsbodenfilteranlagen zur weitergehenden
Regenwasserbehandlung im Misch-und Trennsystem. Hennef,
Germany: DWA — Deutsche Vereinigung für Wasserwirtschaft,
Abwasser und Abfall e.V., 2005. [in German].
Freire FG, Davies LC, Vacas A, Pedro I, Novais JM, Martins-Dias S.
Continuous and intermittent loading of a vertical flow
constructed wetland for an azo dye treatment. In: IWA, editor.
Proceedings of the 10th IWA Specialized Group Conference on“Wetland systems for water pollution control”, 23–29
September 2006, Lisbon, Portugal, 3. ; 2006. p. 1501–9.
García J, Rousseau D, Caselles-Osorio A, Story A, De Pauw N,
Vanrolleghem P. Impact of prior physico-chemical treatment
on the clogging process of subsurface flow constructedwetlands: model-based evaluation. Water Air Soil Poll2007;185:101–9.
Giraldi, D., de'Michieli Vitturi, M., Zaramella, M., Marion, A., Iannelli,R.,. Hydrodynamics of vertical subsurface flow constructedwetlands: tracer tests with rhodamine WT and numericalmodelling. Ecol Eng in press. doi:10.1016/j.ecoleng.2008.06.004.
Giraldi D, de'Michieli Vitturi M, Iannelli R. FITOVERT: a dynamic
numerical model of subsurface vertical flow constructedwetlands. Proceedings of the “International Symposium on
Sanitary and Environment (SIDISA.08)”, 24–27 June 2008, Florence,
Italy; 2008. paper 336 (CD-ROM).
Haberl R, Grego S, Langergraber G, Kadlec RH, Cicalini AR,
Martins Dias S, et al. Constructed wetlands for thetreatment of organic pollutants. JSS — J Soils & Sediments
2003;3(2):109–24.
Henrichs M, Langergraber G, Uhl M. Modelling of organic matter
degradation in constructed wetlands for treatment of combined
sewer overflow. Sci Total Environ 2007;380(1–3):196–209.

3942

S CI EN CE OF T H E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

Henze M, Gujer W, Mino T, van Loosdrecht MCM. Activated sludge
models ASM1, ASM2, ASM2D and ASM3. IWA scientific andtechnical report 9. London, UK: IWA Publishing; 2000.
Kadlec RH. The inadequacy of first-order treatment kinetic
models. Ecol Eng 2000;15:105–19.
Kadlec RH. Effects of pollutant specification in treatment wetlands
design. Ecol Eng 2003;20(1):1–16.
Kadlec RH, Knight RL. Treatment wetlands. New York, USA: Lewis
Publishers; 1996.
Kadlec RH, Knight RL, Vymazal J, Brix H, Cooper P, Haberl R.
Constructed wetlands for pollution control. IWA scientific andtechnical report 8. London, UK: IWA Publishing; 2000.
Kaufmann I, Schmitt TG. Modelling constructed wetlands for CSO
treatment in long-term pollution load simulation. Proceedings
10th International Conference on Urban Drainage, 21–26 August
2005, Copenhagen, Denmark; 2005. CD-ROM.
Korkusuz EA, Meyer D, Langergraber G. CW2D simulation results oflab-scale vertical flow filters filled with special media and loadedwith municipal wastewater. In: Mander Ü, Kóiv M, Vohla C,
editors. 2nd International Symposium on “Wetland pollutant
dynamics and control WETPOL 2007” — extended abstracts,
16–20 September 2007, Tartu, Estonia, vol. II. ; 2007. p. 448–50.
Langergraber G. Development of a simulation tool for subsurface
flow constructed wetlands. Wiener Mitteilungen 2001;169
Vienna, 207 pp.
Langergraber G. Simulation of the treatment performance of
outdoor subsurface flow constructed wetlands in temperate
climates. Sci Total Environ 2007a;380(1–3):210–9.
Langergraber G. Simulation of the treatment performance of a two-stage
CW system in temperate climates. In: Mander Ü, Kóiv M, Vohla C,
editors. 2nd International Symposium on “ Wetland pollutant
dynamics and control WETPOL 2007” — extended abstracts, 16–20
September 2007, Tartu, Estonia, vol. I. ; 2007b. p. 180–2.
Langergraber G. Modeling of processes in subsurface flow constructedwetlands — a review. Vadoze Zone J 2008;7(2):830–42.
Langergraber G, Leroch K, Pressl A, Sleytr K, Rohrhofer R, Haberl R.
High-rate nitrogen removal in a two-stage subsurface verticalflow constructed wetland. Desalination in press.
Langergraber G, Šimůnek J. Modeling variably saturated water flow
and multicomponent reactive transport in constructedwetlands. Vadoze Zone J 2005;4(4):924–38.
Langergraber G, Šimůnek J. The multi-component reactive
transport module CW2D for constructed wetlands for the
HYDRUS Software Package. Riverside, CA, USA: Department ofEnvironmental Sciences, University of California Riverside;
2006. Hydrus Software Series 2, 72pp.
Langmuir D. Aqueous Environmental Chemistry. New Jersey, USA:
Upper Saddle River; 1997. Prentice Hall.
Levenspiel O. Chemical reaction engineering. New York, USA: 3rd
ed. John Wiley and Sons, Inc.; 1999
MacDonald MG, Harbaugh AW. A modular three-dimensional
finite-difference ground-water flow model. U.S. Geological
Survey techniques of water-resources investigations, report99-4259. U.S. Geological Survey, USA; 1988.
Mao X, Prommer H, Barry DA, Langevin CD, Panteleit B, Li L. Three
dimensional model for multi-component reactive transportwith variable density groundwater flow. Environ Modell Softw
2006;21(5):615–28.
Marsilli-Libelli S, Checchi N. Identification of dynamic models forhorizontal subsurface flow constructed wetlands. Ecol Model
2005;187(2–3):201–18.
McGechan MB, Moir SE, Castle K, Smit IPJ. Modelling oxygen
transport in a reedbed constructed wetland purification
system for dilute effluents. Biosyst Eng 2005;91(2):191–200.
Mena J. Tratamiento de aguas residuales urbanas y vínicas mediante
humedales artificiales de flujo subsuperficial. PhD thesis, Universityof Castilla-La Mancha, Ciudad Real, Spain; 2008. [in Spanish].
Mena J, Villaseñor J, Gómez R, De Lucas A, Langergraber G.
Simulation of the hydraulics and the removal efficiency ofpilot-scale horizontal subsurface flow constructed wetlands
with different plants (Phragmites australis and Lythrum salicaria).
In: Mander Ü, Kóiv M, Vohla C, editors. 2nd International
Symposium on “Wetland pollutant dynamics and control
WETPOL 2007” — extended abstracts, 16–20 September 2007,
Tartu, Estonia, vol. I. ; 2007. p. 212–4.
Meyer D, Sommer T, Thomas M, Schmitt TG, Hagen H.
Development of a long-term pollution-load model to simulate
CWs for CSO treatment. In: Mander Ü, Kóiv M, Vohla C, editors.
2nd International Symposium On “Wetland pollutant
dynamics and control WETPOL 2007” — extended abstracts,
16–20 September 2007, Tartu, Estonia, vol. I. ; 2007. p. 217–9.
Meyer D, Langergraber G, Dittmer U. Simulation of sorption
processes in subsurface vertical flow constructed wetlands for
CSO treatment. Proceedings of the 10th IWA SpecializedGroup Conference on “Wetland systems for water pollution
control”, 23–29 September 2006, Lisbon, Portugal; 2006. p.
599–609.
Ogata A. Theory of dispersion in a granular medium. US Geol Surv
Prof Pap 1970; 411-I: 1–34.
Ojeda E, Caldentey J, García J. 2D simulation model for evaluating
biogeochemical pathways involved in organic matter removal
in horizontal subsurface flow constructed wetlands.
Proceedings of the 10th IWA Specialized Group Conference on“Wetland systems for water pollution control”, 23–29
September 2006, Lisbon, Portugal; 2006. p. 1405–13.
ÖNORM B 2505. Bepflanzte Bodenfilter (Pflanzenkläranlagen) —Anwendung, Bemessung, Bau und Betrieb (Subsurface-flow
constructed wetlands — Application, dimensioning,
installation and operation). Vienna, Austria: Österreichisches
Normungsinstitut; 2005 [in German].
ÖNORM EN 12556-3. Kleinkläranlagen für bis zu 50 EW. Teil 3:
Vorgefertigte und/oder vor Ort montierte Anlagen zur Behandlung
von häuslichem Schmutzwasser (Small wastewater treatment
systems up to 50 PT. Part 3: Packaged and/or site assembled
domestic wastewater treatment plants). Vienna,
Austria: Österreichisches Normungsinstitut; 2001. in German.
Parkhurst DL, Appelo CAJ. User's guide to PHREEQC — a computer
program for speciation, reaction-path, 1D-transport, andinverse geochemical calculations. Technical reports. USA: U.S.
Geological Survey; 1999. p. 99–4529.
Pastor R, Benqlilou C, Paz D, Cardenas G, Espuña A, Puigjaner L.
Design optimisation of constructed wetlands for wastewater
treatment. Resour Conserv Recy 2003;37(3):193–204.
Persson J, Somes NLG, Wong THF. Hydraulic efficiency ofconstructed wetlands and ponds. Water Sci Technol
1999;40(3):291–300.
Rousseau DPL. Performance of constructed treatment wetlands:
modelbased evaluation and impact of operation andmaintenance. PhD thesis, Ghent University, Belgium; 2005.
Rousseau DPL, Vanrolleghem PA, De Pauw N. Model-based designof horizontal subsurface flow constructed treatment wetlands:a review. Water Res 2004;38(6):1484–93.
Sandoval-Cobo JJ, Peña MR. Performance analysis of SSF wetlands
in tropical regions based on first-order kinetic models for theremoval of organic matter. Ecol Eng submitted for publication.
Šimůnek J, Šejna M, van Genuchten MTh. The HYDRUS software
package for simulating the two- and three-dimensional
movement of water, heat, and multiple solutes invariably-saturated media — user manual, version 1.0.
PC-Progress, Prague, Czech Republic; 2006a. 160p.
Šimůnek J, Šejna M, van Genuchten MTh. The HYDRUS software
package for simulating the two- and three-dimensional
movement of water, heat, and multiple solutes invariably-saturated media — technical manual, version 1.0.
PC-Progress, Prague, Czech Republic; 2006b. 241p.
Stein OR, Biederman JA, Hook PB, Allen WC. Plant species andtemperature effects on the k–C⁎ first-order model for COD
removal in batch-loaded SSF wetlands. Ecol Eng 2006;26(2):100–12.

S CI EN C E OF TH E T OTAL EN V I RO N M EN T 4 0 7 ( 2 0 09 ) 39 3 1–3 94 3

Straub WA, Lynch DR. Models of landfill leaching: moisture flow
and inorganic strength. J Environ Eng-ASCE 1982;108(2):231–50.
Tomenko V, Ahmed S, Popov V. Modelling constructed wetland
treatment system performance. Ecol Model
2007;205(3–4):355–64.
Toscano A, Langergraber G, Consoli S, Cirelli GL. Modelling pollutant
removal in a pilot-scale two-stage subsurface flow constructedwetlands. Ecol Eng in press. doi:10.1016/j.ecoleng.2008.07.011.
Uhl M, Dittmer U. Constructed wetlands for CSO treatment — an
overview of practice and research in Germany. Water Sci
Technol 2005;51(9):23–30.
Vymazal J, Brix H, Cooper PF, Green MB, Haberl R, editors.
Constructed wetlands for wastewater treatment in Europe.
Leiden, The Netherlands: Backhuys Publishers; 1998.

3943

Wanko A, Mose R, Carrayrou J, Sadowski AG. Simulation of
biodegradation in infiltration seepage — model development
and hydrodynamic calibration. Water Air Soil Poll
2006;177:19–43.
Wynn M, Liehr SK. Development of a constructed subsurface-flowwetland simulation model. Ecol Eng 2001;16:519–36.
Zheng C, Wang PP. A modular three-dimensional multi-species
transport model for simulation of advection, dispersion andchemical reactions of contaminants in groundwater systems;
documentation and user's guide. U.S. Army Engineer Researchand Development Center Contract Report SERDP-99-1,
Vicksburg, Mississippi, USA; 1999.