Ecological Engineering 68 (2014) 105–115

Contents lists available at ScienceDirect

Ecological Engineering
journal homepage: www.elsevier.com/locate/ecoleng

The veriﬁcation of the Constructed Wetland Model No. 1
implementation in HYDRUS using column experiment data
T.G. Pálfy ∗ , G. Langergraber ∗Institute of Sanitary Engineering and Water Pollution Control, University of Natural Resources and Applied Life Sciences (BOKU), Vienna, Austriaa r t i c l e

i n f oArticle history:
Received 5 November 2013
Received in revised form 9 January 2014
Accepted 24 March 2014
Keywords:
Batch operation
Constructed wetland
Column experiment
CWM1
Process-based modela b s t r a c t
The Constructed Wetland Model No. 1 (CWM1) is a biokinetic model describing microbial dynamics andtransformation and degradation processes in subsurface ﬂow constructed wetlands. The implementationof CWM1 in the HYDRUS program package (© PC Progress s.r.o.) was veriﬁed using data from controlledenvironment column experiments. These 20-day long batch experiments used synthetic wastewater andwere run at different temperatures. Columns were either planted with Carex rostrata Stokes, Schoenoplectus acutus Muhl. Ex Bigelow and with Typha latifolia L. or without plants. The minimal adjustment ofmodel parameters between the sixteen setups was targeted and (1) initial bacteria concentrations, (2)initial adsorbed ammonium nitrogen concentrations, and (3) root oxygen loss rates inevitably needed to
be set for each columns separately. For all other parameters – including the CWM1 model parameters –identical values have been used. Some biokinetic parameters had to be iteratively adjusted during calibration to allow anaerobic, anoxic and aerobic processes to occur in parallel in order to match measured data.
This can be justiﬁed by the local effect of root zone re-aeration. The simulation results were evaluatedby a new goodness of ﬁt analysis method, deﬂection analysis. Overall, simulated contaminant concentrations had a very good ﬁt to measured NH4 –N and SO4 –S concentrations, and a good ﬁt to measured COD
concentrations.
© 2014 Elsevier B.V. All rights reserved.

1. Introduction
Constructed wetlands (CWs) are state-of-the-art extensive
solutions for wastewater treatment (Campbell and Ogden, 1999).
They are capable of tackling a wide range of pollutants at low
maintenance and energy costs. The quality of the efﬂuent is comparable to secondary or tertiary efﬂuents of other technologies(Campbell and Ogden, 1999; Kadlec and Wallace, 2009). Subsurface ﬂow (SSF) CWs have the common feature that the water levelis kept below the surface of the ﬁlter bed. This wet media supports bacterial and bioﬁlm growth and also physical and chemicalcleansing processes on the surface of the soil particles. Degradation
processes are supported by plant uptake, root oxygen and exudate
transfer.
In the recent years, CWs have become increasingly important in several regions of the world. More efﬁcient designs have

∗ Corresponding authors. Tel.: +43 1 47654 5814; fax: +43 1 47654 5803.
E-mail addresses: palfy.tamas@gmail.com (T.G. Pálfy),
guenter.langergraber@boku.ac.at (G. Langergraber).
http://dx.doi.org/10.1016/j.ecoleng.2014.03.016
0925-8574/© 2014 Elsevier B.V. All rights reserved.

been achieved and the increasing knowledge and understanding
has led to numerical models similar to Activated Sludge Models (ASMs) (Henze et al., 2000) designed to simulate biologicaltreatment processes in SSF CWs (Langergraber, 2008). Numericalmodelling can provide better understanding of CW processes andoperation needs through facilitating insight in the functionality anddynamics of these complex systems, similarly to physical models
but, in many cases as a more ﬂexible and cost and time effective
option.
The implementation of the biokinetic model Constructed Wetland Model No. 1 (CWM1) (Langergraber et al., 2009) in the HYDRUS
Wetland Module was selected for veriﬁcation. HYDRUS-CWM1 is
one of the most comprehensive process-based models for SSF CWs
that are available (Langergraber and Šimůnek, 2012). A series ofphysical column experiments carried out at Montana State University (Allen et al., 2002) was compared with the results of theirnumerical counterparts simulated in CWM1. As the physical column experiments were well-known in details and carried out atvarious temperatures and with different vegetation settings, theseprovided the proper material to investigate the functionality andissues around HYDRUS-CWM1.

106

T.G. Pálfy, G. Langergraber / Ecological Engineering 68 (2014) 105–115

2. Methods
2.1. The physical column experiments
The source data and the model setup was based on a paper about
the physical column experiments carried out at Montana State University (Allen et al., 2002), more detailed results of measurements
were obtained from Mburu (2013). The greenhouse experimentused column-shaped constructed wetland microcosms. The experiment took place in Bozeman, Montana (N 45◦ 40 W 111◦ 03 ,at 1490 m elev.) and lasted about four years. A series of 20-day
long batch operations was conducted using synthetic wastewaterwith known constituents. For model veriﬁcation only the columns
with temperatures ranging from 12 to 24 ◦ C were modelled. No
artiﬁcial light sources were used, the cumulative daily net solar
radiation ranged from 1 to 8 MJ m−2 day−1 representing about
25% of the locally measured radiation because a fabric light ﬁlter
was employed to improve temperature control. Relative humidity
ranged from 30 to 70% with no seasonal pattern. Plant growth was
robust and plant activity showed seasonal cycles.
Eight replicates were built for each of the unplanted and planted
columns. Columns were planted with beaked sedge (Carex rostrata Stokes), hardstem bulrush (Schoenoplectus acutus Muhl. Ex
Bigelow) and broadleaf cattail (Typha latifolia L.) In this paper thecolumns are referred to as unplanted, Carex, Schoenoplectus andTypha. Each of them was made of a polyvinyl chloride pipe with60 cm height and 20 cm diameter ﬁlled to a depth of 50 cm. Alluvial pea gravel (D = 0.3–1.3 mm) of igneous and metamorphic origin
was used as ﬁlter material with a porosity of 0.27 and a porevolume of 4.3 L per column. Three access tubes were installed
but only the one ending at 15 cm below the ﬁlter surface wasused because no depth-correlated concentration differences were
measured. Evaporation losses were replenished with greenhousetemperature de-chlorinated tap water added to the bottom of
the columns. Evaporative losses were 0.7 L day−1 per column at
24 ◦ C and about 0.4 L day−1 per column at 4 ◦ C. The difference ofevapotranspiration losses between different plant species was notmeasured. Measured experiments began after plants had grown
to reasonable size on nutrient solutions and later on synthetic
wastewater.
A synthetic wastewater imitating secondary domestic efﬂuent
was used. Mean inﬂuent concentrations were: COD 470 mg L−1 ,Total N 44 mg L−1 , Amino-N 27 mg L−1 , NH4 –N 17 mg L−1 , SO4 –S
14 mg L−1 , PO4 –P 8 mg L−1 . Solution concentration data was available from directly after ﬁlling and from day 1, 3, 6, 9, 14 and 20 andwere used as a basis of comparison to simulation results. For moredetailed information about the physical experiments we refer toAllen et al. (2002).
2.2. HYDRUS-CWM1
The implementation of CWM1 comes with the HYDRUS Wetland Module which is an extension of the standard HYDRUS
features allowing detailed simulation of SSF CWs (Langergraber andŠimůnek, 2012). HYDRUS numerically solves systems of partial differential equations in three dimensions to simulate (Šimůnek et al.,2011): (1) variable-saturated water ﬂow as described by Richard’s
equation, (2): transport of constituents by solving mass balance
equations; (3) inﬂuence of plants via root water and nutrient uptake
and oxygen loss; and (4) heat transport.
CWM1 is one of the two biokinetic models which the HYDRUS
Wetland Module adds to the above described HYDRUS framework. It is capable to model bacteria dynamics and several
anaerobic, anoxic and aerobic transformation processes involvingorganic matter (expressed as COD), nitrogen and sulphur. Table 1

Table 1
The components of CWM1 and the phases they are deﬁned in L: liquid phase and S:
solid phase (Langergraber and Šimůnek, 2012).

1.
2.
3.
4.
5.
6.
7.
8.

9.
10.
11.
12.
13.
14.
15.
16.

SymbolSoluble components

Phase

SO
SF
SA
SI
SNH
SNO
SSO4
SH2S

Dissolved oxygen, O2
Fermentable, readily biodegradable soluble COD
Fermentation products as acetate
Inert soluble COD
Ammonium and ammonia nitrogen
Nitrate and nitrite nitrogen
Sulphate sulphur
Dihydrogensulphide sulphur

L
L+S
L+S
L+S
L+S
L
L
L

SymbolParticulate components

Phase

XS
XI
XH
XA
XFB
XAMB
XASRB
XSOB

Slowly biodegradable particulate COD
Inert particulate COD
Heterotrophic bacteria
Autotrophic nitrifying bacteria
Fermenting bacteria
Acetotrophic methanogenic bacteria
Acetotrophic sulphate reducing bacteria
Sulphide oxidizing bacteria

L+S
L+S
S
S
S
S
S
S

enumerates the model components categorized to soluble and particulate groups and describes in which phase these can be present
in the model. Bacteria are considered to be immobile and are therefore only deﬁned in the solid phase. The seventeen transformation
processes considered in CWM1 are resumed in Table 2. For more
information about CWM1 and the HYDRUS Wetland Module werefer to Langergraber et al. (2009) and Langergraber and Šimůnek(2011, 2012), respectively.
2.3. Model setup
The model was parameterized based on the descriptions of Allenet al. (2002) and, for initial runs before the iterative adjustments
which took place, based on some of the results and suggestionsof Mburu et al. (2012) who implemented CWM1 in the AQUASIM
software (Reichert, 1998). Where no data was available estimations
were made based on similarities to other experiments or materials

Table 2Modelled biokinetic processes in CWM1 (Langergraber and Šimůnek, 2012).
Heterotrophic bacteria
Hydrolysis: conversion of XS to SF
1.
Aerobic growth of XH on SF (mineralization of organic matter)2.
Aerobic growth of XH on SA (mineralization of organic matter)3.
Anoxic growth of XH on SF (denitriﬁcation)
4.
Anoxic growth of XH on SA (denitriﬁcation)
5.
6.
Lysis of XH
Autotrophic bacteria
Aerobic growth of XA on SNH (nitriﬁcation)
7.
Lysis of XA
8.
Fermenting bacteria
Growth of XFB (fermentation)
9.
Lysis of XFB
10.
Acetotrophic methanogenic bacteria
11.
Growth of XAMB (anaerobic growth on SA)
Lysis of XAMB
12.
Acetotrophic sulphate reducing bacteria
Growth of ASRB (anerobic)
13.
Lysis of XASRB
14.
Sulphide oxidizing bacteria
Aerobic growth of XSOB on SH2 S (oxidation to SSO4 )15.
Anoxic growth of XSOB on SH2 S (oxidation to SSO4 )16.
Lysis of XSOB
17.

T.G. Pálfy, G. Langergraber / Ecological Engineering 68 (2014) 105–115

Fig. 1. Comparison of the simulated domain mesh and its boundary conditions withthe setup of physical experiments (Allen et al., 2002). Boundaries: (a) atmosphericboundary, (b) constant head boundary, all other: no ﬂux boundary.

discussed in Henze et al. (2000), Langergraber (2001), Salvato et al.
(2003) and Kadlec and Wallace (2009).
Columns were modelled in a 10 cm × 50 cm 2D axisymmetricvertical ﬂow domain representing the radius and the ﬁlling depthof the physical column, respectively. Simulated processes included
water ﬂow, solute and heat transport and root processes. Maximum time step was set to 10−3 day. The domain was represented
via a ﬁnite element mesh consisting of 31 rows and 21 columns.
The top boundary was set to be atmospheric for modelling evapotranspiration and the two nodes at the bottom of the axis were set
to constant head boundary, responsible for the replenishment of
evapotranspiration losses. Other boundaries were set to no ﬂux.
Fig. 1 compares the simulated 2D domain mesh and its boundary
conditions with the setup of the physical experiments.
Evapotranspiration losses were linearly interpolated for all
planted columns using approximated values at 4 and 24 ◦ C
described by Allen et al. (2002). For unplanted columns correlations
were adopted from Kadlec and Wallace (2009: 110) to estimate theorder of magnitude of evaporation losses from unplanted media.
The porosity of the media was set to the measured value of
0.27. No hydraulic conductivity data was found so a value of
172,800 cm day−1 was selected after Salvato et al. (2003) whichis representative to pea gravel. This value has signiﬁcance eventhough the columns were batch operated because the evapotranspiration losses represent up to 8.26 cm depth of water in the porous
media per column per day. Filter material and root distribution
were set homogenously in all cases and observation nodes were set
to 5, 10 and 30 cm which did not show any difference in concentration via depth just like in the case of the physical experiments andthus data from the 10 cm deep observation nodes were selected
for further analysis. To ensure proper mixing, a high dispersion
coefﬁcient had to be selected.
The simulated domain was water saturated. Initial sorbed
ammonium concentrations were set after Mburu et al. (2012) but
were modiﬁed later in all instances. Initial dissolved oxygen concentration was set to 1 mg L−1 . Although the distilled water used
for producing synthetic wastewater should not contain any oxygen the mixing process of the synthetic contaminants probably

107

involved exposure to air. The low initial DO concentration was represented well by the initial measured sulphate dynamics. Nitrous
nitrogen, dihydrogen-sulphide sulphur and inert particulate COD
concentrations were set to zero; other initial contaminant concentrations were set after measured data provided by Mburu (2013).
Initial bacteria concentrations ranged on a 10−2 –100 ␮g g−1 order.
Initial adsorbed SNH was adjusted until a good ﬁt to measured
ammonium-nitrogen concentrations was achieved. Values of initial conditions used for the simulations are included in the Annex
in Table 14.
The initial COD concentration was proportionated between the
different fractions of SA, SF, SI and XS using the recommendations ofMburu et al. (2012) which is a slightly modiﬁed version of the ASM
models as no particulate organic matter was present in the synthetic wastewater. XI was set to zero, SF was estimated to be 25%,
SA 10%, SI 4% and XS 61% of the total COD. As Mburu et al. (2012) we
did not consider adsorption of COD fractions. Organic nitrogen isconsidered in CWM1 as part of the COD and was allocated to different COD fractions. Composition parameters for SF and SI were usedas given by Langergraber et al. (2009), which is 0.03 and 0.01 g N g−1
COD, representing about 3.5 mg L−1 and 0.2 mg L−1 N, respectively.
The rest of the organic nitrogen was allocated to XS by increasing
its nitrogen content to about 0.084 g N g−1 COD.
2.4. Calibration procedure
Biokinetic parameters were left as described by Langergraberand Šimůnek (2011) for early simulations except those which werechanged by Mburu et al. (2012). The unplanted columns were calibrated ﬁrst because they had simpler concentration dynamics asplanted ones. The second group to continue calibration with wascolumns planted with Carex as they had the most complex dynamics of measured sulphate concentrations. It was expected that it
will need the most effort to reproduce the results of these columns
and once the simulation results ﬁt the measured data well no further adjustment of the parameters will be required to simulate
Schoenoplectus and Typha columns.
Simulation results of unplanted columns were in good agreement with measured data with unchanged biokinetic parameters.
The measured values of NH4 –N and SO4 –S showed these columns
were dominated by anaerobic processes with no denitriﬁcation andcomplete sulphate reduction. Contrary to this, columns plantedwith Carex showed signs of hosting both anaerobic and aerobic
processes by simultaneous nitriﬁcation and sulphate reduction.
To match measured concentrations in the Carex columns, several inhibition/half-saturation concentration parameters had to bechanged to allow these processes happen parallel in the modelled
domains as well. After these changes were done, only (1) initialbacteria concentrations, (2) root oxygen release, (3) initial adsorbed
SNH concentration had to be set up for simulated Schoenoplectus
and Typha domains and to ﬁne-tune previous results.
1. Values for initial bacteria concentrations given by Mburu et al.
(2012) were reduced by a factor of 10 to avoid immediatebacteria die-off after the start of the simulation. Initial bacteria
concentrations were also used as a ﬁnal ﬁne-tuning parameterat the end of the calibration procedure.
2. Root oxygen release was calibrated based on the SSO4 and SNH
concentration dynamics. A small value would have already led toa minor SO concentration which would have inhibited anaerobic
processes in the system if the standard biokinetic parameters
had been used.
3. Initial adsorbed SNH concentrations had to be set to be able tomodel both the initial concentration increase of dissolved SNH
in the case of unplanted and the decrease of dissolved SNH in

108

T.G. Pálfy, G. Langergraber / Ecological Engineering 68 (2014) 105–115

Table 3
Values of biokinetic model parameters used for the simulations compared to values given by Langergraber et al. (2009) and used by Mburu et al. (2012).
Parameter name [unit]−1

Hydrolysis rate constant [day ]Half-saturation coefﬁcient of autotrophic bacteria for SO [mg L−1 ]Maximum growth rate for autotrophic bacteria [day−1 ]Maximum growth rate for fermenting bacteria [day−1 ]Inhibition coefﬁcient of acetotrophic methanogenic bacteria for SO [mg L−1 ]Inhibition coefﬁcient of acetotrophic methanogenic bacteria for SNO [mg L−1 ]Maximum growth rate for acetotrophic sulphate reducing bacteria [day−1 ]Inhibition coefﬁcient of acetotrophic sulphate reducing bacteria for SO [mg L−1 ]Inhibition coefﬁcient of acetotrophic sulphate reducing bacteria for SNO [mg −1 ]Half-saturation coefﬁcient of acetotrophic sulphate reducing bacteria for SSO4 [mg L−1 ]Yield coefﬁcient for acetotrophic methanogenic bacteria
Yield coefﬁcient for acetotrophic sulphate reducing bacteria

the case of the planted columns based on a realistic equilibrium
isotherm and bacteria concentrations.
The calibration process is described in detail in Pálfy (2013).
2.5. Data evaluation
The simulation results were evaluated using a standardized
MS Excel© workbook which enabled the rapid evaluation of thehigh number of different outcomes from the partly trial-and-error
approach and because the modelling tool handles COD in different
fractions and does not give any graphical or summarized numericalresults for it. The ﬁnal worksheets contained calculus and additional diagrams for different goodness of ﬁt analysis approaches.
A new goodness of ﬁt analysis method, deﬂection analysis (Pálfy,2013) was used to evaluate simulation results because the generally used coefﬁcient of determination (R2 ) and coefﬁcient of
efﬁciency (E) as suggested e.g. by Ahnert et al. (2007) have beenshown to be inadequate for evaluating the ﬁt of the simulations
to measured data given by mean and standard deviation (SD).
Deﬂection analysis is comparing the proportion of the area set
by the deviation of simulation results from the region deﬁned bymean ± SD to the area determined by mean ± SD values (deﬂection, DEF) and to the area below measurement means (magnitudeof deﬂection, MAG), using the trapezoidal rule of numerical integration. Fig. 2 shows the visual interpretation of DEF and MAG.
The deviation (DEVcum ) which is represented on both sides ofFig. 2 by light grey colour can be calculated using Eq. (1):
DEVcum ≈n−1
(ti+1 − ti )i=1

 x̄ + SD + x̄s + si+1 
 i
i
i+1
− i


4
4





SDi,i+1
 x̄ − SDi + x̄i+1 − SDi+1 si + si+1 
+ i
−
 − 2
4
4

(1)

The area enclosed by mean ± SD (SDcum ), marked by dark grey
colour on the left part of Fig. 2 is calculated using Eq. (2):SDcum ≈n−1
(ti+1 − ti )(SDi + SDi+1 )

(2)

i=1



Langergraber et al. (2009)

Mburu et al. (2012)

This work

Kh
KOA
␮A
␮FB
KOAMB
KNOAMB␮ASRB
KOASRBKNOASRB
KSOASRB
YAMB
YASRB

3.0
1.0
1.0
3.0
0.00020.00050.18
0.00020.0005
19.0
0.032
0.05

0.58

0.58
0.0026
1.50
3.77
0.02
0.05
0.50
0.05
0.50
5.00
0.04
0.04

3.77

0.31

0.04
0.04

whilst the magnitude of deﬂection as
MAG =DEVcum
.
MEScum(5)

Simulations with DEF > 50% and MAG > 25% were considered as
weak results. These threshold values have been derived from theevaluation of our simulation results and might need to be revised
after more practical experience with using this method is available.
3. Results and discussion
3.1. Model calibration
Table 3 shows those biokinetic parameters at 20 ◦ C for which different values were used compared to those given by Langergraber
et al. (2009) and by Mburu et al. (2012). Tables of all adjusted model
parameters and initial concentrations can be found in the Annex.
The changes were the most signiﬁcant in the half saturationand inhibition coefﬁcients for dissolved oxygen and nitrate (SO and
SNO). These were increased or decreased by about two logs making
possible simultaneous nitriﬁcation and sulphate reduction which
can occur within a single ﬁnite element (Langergraber, 2011). Rootsurfaces can act as a second type environment in the ﬁlter media
favouring for anoxic or aerobic processes. In the model there isno spatial segregation of root surface zones thus adjusting these
coefﬁcients was the way to include the effects of roots on microbialactivity. These settings promoted microbial activity of XA becauseof much lower half-saturation coefﬁcient for SO and of XAMB and
XASRB because of much higher inhibition concentrations for SO
and SNO. Without these settings the simulated activity of thesemicroorganisms would have been close to zero. The processes
linked to these bacteria groups were necessary to follow the SO4 –S
and NH4 –N dynamics of the real columns. Fig. 3 shows the effectof a higher inhibition coefﬁcient for SO (0.0002 → 0.05 mg L−1 ) andSNO (0.0005 → 0.50 mg L−1 ) and the lower half-saturation coefﬁcient for the substrate SSO4 (19.0 → 5.0 mg L−1 ) of ASRB on thesimulation results of Carex 24 ◦ C SO4 –S.
3.2. Simulation results

The integral below measurement means (MEScum ), marked by darkgrey colour on the right part of Fig. 2 can be estimated by Eq. (3):n−1

MEScum ≈Symbol(ti+1 − ti )i=1

x̄i + x̄j+1
2(3)

As deﬁned in Pálfy (2013) Deﬂection can be simply calculated as
DEVcum
DEF =
SDcum

(4)

Simulated COD concentrations are shown in Fig. 4, the results ofdeﬂection analysis in Table 4. Most of the simulated data ﬁts well tothe measured data; however deviations occur in the case of several
columns. These can be categorized in two groups:
1. Simulated concentrations do not follow the dynamics between
day one and six which can be described as a slight increase or
slowdown in the decrease of concentrations. Columns where thisT.G. Pálfy, G. Langergraber / Ecological Engineering 68 (2014) 105–115

109

Fig. 2. Visual explanation of deﬂection (DEF) and magnitude (MAG). DEF is the proportion of the light grey ﬁeld and the diagonally hatched area (measurement means ± SD, left).
MAG is the proportion of the light grey ﬁeld and the horizontally hatched area (area below measurement means, right). Polygon markers are representing the measurementmeans, the dark grey line the simulation results.

SO4-S

concentraon[mg L-1 ]20

SO4-S
Measured concentraons
Model results

20

15

15

10

10

5

5

0
0

10

5

15

20
me [days]

Measured concentraons
Model results

0
0

10

5

15

20
me [days]

Fig. 3. The combined effect of the ﬁnal, higher inhibition coefﬁcient for SO and SNO and the ﬁnal, lower half-saturation coefﬁcient for the substrate SO4 –S of XASRB on thesimulated SO4 –S concentrations of Carex 24 ◦ C during calibration. Left: original parameter values from Langergraber et al. (2009), right: parameter values adjusted.

problem manifests most pronounced are Unplanted 12 ◦ C andTypha 12 ◦ C.
2. Simulated concentrations underestimate COD removal from daynine or through the whole batch, seemingly caused by an underestimated initial drop of the concentrations within day one.
Columns where this problem is most notable are Carex 16, 20
and 24 ◦ C.
The results of deﬂection analyses for Carex 16, 20 and 24 ◦ C andSchoenoplectus 16 ◦ C have values of DEF > 50% and MAG > 25% thus
the simulation outcomes does not match measured data properly.
The weak ﬁt can be visually seen on Fig. 4 as the ﬁnal COD concentrations are overestimated.
Simulated and measured NH4 –N concentrations for all columns
are shown in Fig. 5, the results of deﬂection analysis in Table 5.
Table 4
The results of deﬂection analysis (in %) on simulated COD concentrations. The four
weakest results of all simulations were all COD (highlighted in grey).
Unplanted

Carex

Schoenopect.

Typha

Model results ﬁt very well to the measured data. The relatively
simple dynamics of ammonium nitrogen concentrations could be
reproduced well by changing initial adsorbed NH4 –N concentrations to scale the initial decrease or rise. Later changes could be
simulated by very slow nitriﬁcation and root nutrient uptake.
Table 5 shows the results of deﬂection analysis for NH4 –N,
marking Carex 20 ◦ C as column with the weakest ﬁt, with 49% DEF
and 40% MAG. In this case the low SD and measurement meansincreased the sensitivity of the method to any deviation from measurements. Schoenoplectus 24 ◦ C has a bad simulated ﬁt on day one
but because of the averagely high SD this is not considered lesssigniﬁcant by deﬂection analysis.
Fig. 6 shows simulated and measured SO4 –S concentrations,
Table 6 the results of deﬂection analysis. All simulation results
describe well sulphate dynamics in anaerobic columns witha quick transformation into hydrogen sulphide. Columns with
Table 5
The results of deﬂection analysis (in %) on simulated NH4 –N concentrations.
Unplanted

Carex

Schoenoplect.

Typha

DEF
MAG

12 ◦ C

16
8

12
8

19
10

19
9

DEF
MAG

12 ◦ C

0
0

4
8

1
1

1
1

DEF
MAG

16 ◦ C

0
0

169
93

75
27

9
4

DEF
MAG

16 ◦ C

0
0

11
15

0
0

1
1

DEF
MAG

20 ◦ C

0
0

68
44

62
20

15
4

DEF
MAG

20 ◦ C

0
0

49
40

0
0

0
1

DEF
MAG

24 ◦ C

0
0

53
26

45
14

2
1

DEF
MAG

24 ◦ C

0
0

9
10

2
2

0
0

110

T.G. Pálfy, G. Langergraber / Ecological Engineering 68 (2014) 105–115

Fig. 4. Simulation results for COD using CWM1 in HYDRUS. Continuous line: simulation results, dots: measurement means, intervals: standard deviation.

better root re-aeration include Schoenoplectus 12 ◦ C and Carex 12 ◦ C
and 16 ◦ C and are showing higher SO4 –S concentrations. The latter two had a complex dynamics which was simulated well until
about day thirteen. After this point values still stay in the range
set by SD but it is questionable if the processes are described properly because the simulated concentrations are decreasing while themeans of the measurements are increasing. In the model available
SNH is depleted and this results in bacterial inactivity and die-offof sulphide oxidizers, whilst in reality adsorbed nitrogen might bedirectly available in the bioﬁlm.
Table 6
The results of deﬂection analysis (in %) on simulated SO4 –S concentrations.

Table 6 shows simulated results match well measured ones. The
weakest result is the unplanted 12 ◦ C column with 37% DEF and 14%
MAG. Considering the low SD values, measurement means and thevisual interpretation this can still be considered as a very good ﬁt.
Visual goodness of ﬁt evaluation shows in general a very good
representativeness of the simulation results of the CWM1 implementation in HYDRUS. This is in accordance with the average ofthe results of deﬂection analyses carried out. Table 7 summarizes
these results by contaminant species for all, for unplanted and forplanted columns and includes values for an overall ﬁt (column“All”). The COD of the planted columns has the worst goodness ofTable 7
Average goodness of ﬁt (in %) by contaminant species for all, for unplanted and forplanted columns.

Unplanted

Carex

Schoenopl.

Typha

12 ◦ C

37
14

0
0

0
0

13
14

COD

NH4 –N

SO4 –S

All

DEF
MAG

16 ◦ C

1
2

0
0

2
2

0
1

All columns

DEF
MAG

35
17

5
5

5
5

15
9

DEF
MAG

20 ◦ C

5
7

9
15

5
10

2
2

Unplanted columns

DEF
MAG

4
2

0
0

11
6

5
3

DEF
MAG

24 ◦ C

2
3

2
4

0
0

0
0

Planted column

DEF
MAG

46
22

7
7

3
4

18
11

DEF
MAG

T.G. Pálfy, G. Langergraber / Ecological Engineering 68 (2014) 105–115

111

Fig. 5. Simulation results of NH4 –N using CWM1 in HYDRUS. Continuous line: simulation results, dots: measurement means, intervals: standard deviation.

ﬁt between measured and simulated values, close to the threshold of 50% DEF and 25% MAG. Unplanted columns had simpler
contaminant dynamics compared to planted ones and their simulated results for all three contaminant group ﬁt very well measured
values.

3.3. Issues of model setup and suggestions for improvement
Columns represent a physical model type where there is excessive growing space for aboveground biomass compared to plot or
full-scale experiments. The humidity can be also different compared to the one in the microclimate formed by large-scale stands.
In addition, radiation was only about 25% of locally measured (Allenet al., 2002). These can change plant physiology and have an implication on root nutrient uptake and oxygen and exudate supply.
To minimize these effects, planted columns should be arranged
as close to each other as possible and surrounded by a darkened
foil as high as the plants are to mimic microclimate and full-scale
habitat.
The reliability of calibrated root oxygen supply is impoverished because the water used in the physical experiments forevapotranspiration loss replenishment was not de-oxygenated.
Table 8 shows the potential oxygen input through this pathway andcompares it to the simulated cumulated ﬂux supplied by roots inthe case of the 24 ◦ C Carex model setup. Evapotranspiration losses
were not measured by plant species which means oxygen entering
the microcosms via water replenishment could have had different effect on columns with different plants. De-oxygenated water
should be preferred for water loss compensation.
There is a large uncertainty about the three parameters calibrated separately for each domain setup in CWM1: (1) initialbacteria concentrations, (2) initial adsorbed ammonium nitrogenconcentrations and (3) root re-aeration. It would be of high value for
CWM1-based modelling if these parameters could have been measured, just as knowing the bulk density and hydraulic permeabilityof the used materials.

Table 8Potential oxygen input through evaporative loss compensation.
Assumed DO
concentration in
replenishing water[mg L−1 ]Cumulative DO ﬂux
during batch [mg]Cumulative DO ﬂux
expressed in % of thesimulated root oxygen
release (%)

0.1
0.5
1.0
2.0
3.4

1.4
7.0
14.0
28.0
47.6

3
15
29
59
100

112

T.G. Pálfy, G. Langergraber / Ecological Engineering 68 (2014) 105–115

Fig. 6. Simulation results of SO4 –S using CWM1 in HYDRUS. Continuous line: simulation results, dots: measurement means, intervals: standard deviation.

4. Summary and conclusions
Experiments with batch-operated constructed wetland microcosms were carried out with different plant species at known
temperatures and synthetic wastewater (Allen et al., 2002). Mburu
et al. (2012) simulated data from these experiments using animplementation of the Constructed Wetland Model 1 (CWM1,
Langergraber et al., 2009) in the AQUASIM software (Reichert,1998) which has a model parameter-estimation function. The work
presented in this paper has veriﬁed the implementation of CWM1
in the HYDRUS Wetland Module (Langergraber and Šimůnek, 2012)using the same experimental data as Mburu et al. (2012). These
results show that parameters of CWM1 can be calibrated manuallyby starting with extremities using a trial-and-error approach. Firstthe simplest column, the anaerobic unplanted was parameterized,then the one with the most complex SO4 –S dynamics, Carex 12 ◦ C.
After further minor adjustments to unplanted columns and to all
Carex columns no changes were made to ﬁt any of the Schoenoplectus or Typha columns. This means the ﬁnal parameter set worked
well in describing various anaerobic and anoxic columns with different plant species and different temperatures.
After deﬁning a parameter set the goal was to adjust as few
parameters as possible between the simulations of each column
to achieve the desired goodness of ﬁt. In contrast to Mburu et al.
(2012) no adsorption parameters for SNH or slowly biodegradable
COD have been adjusted for different columns. It was satisfactory
to adjust only initial SNH and bacteria concentrations in the solid
phase and root oxygen loss parameter.
The deﬁned parameters allow aerobic and anaerobic processes to occur in parallel after the adjustment of some of thehalf-saturation and inhibition coefﬁcients. This can be justiﬁed
considering root re-aeration and aerobic conditions at root surfaces
which were not spatially discretized but modelled with anaerobic
processes in the elements of the mesh.
Deﬂection analysis has been introduced as a goodness of ﬁt
analysis method capable to take into account the standard deviation of measured data. It was shown by using this methodthat the HYDRUS implementation of CWM1 was working ﬁne
and provided an overall good ﬁt that followed the contaminant
dynamics of the measured concentrations in the batch-operated
columns. The weakest match have been the simulations for COD
for the Carex 16 ◦ C, 20 ◦ C and 24 ◦ C and Schoenoplectus 20 ◦ Cmicrocosms.
Model results for root oxygen release cannot be considered asreliable due to some aspects of the physical experiment. The oxygen
content of the water used for the compensation of evapotranspiration losses was unknown as well as the exact evapotranspiration
rates by column type and temperature. Simulated root oxygen
release includes also the oxygen intake by water used for thecompensation of evapotranspiration losses. The use of distilled orde-oxygenated water is recommended for water replenishment inphysical column experiments.

Acknowledgements
Tamás Pálfy thanks Zoltán Gribovszki for guiding into theﬁeld of constructed wetlands. The authors thank Otto R. Stein,Montana State University, for providing additional detailed information regarding the column experiments.

T.G. Pálfy, G. Langergraber / Ecological Engineering 68 (2014) 105–115
Table A.2
Adsorption parameters.

Annex. Adjusted model parameters
Tables A.1–A.5.

Table A.1
Domain, water ﬂow and solute transport parameters.
Unit

Domain geometry
Type of geometry
2D-domain options

RadiusHeightHorizontal discretization in x

cm
cm

Vertical discretization in yTime information
Initial time step
Minimum time step
Maximum time stepdddSoil hydraulic model
Hydraulic model

–
–
cm−1
–
cm d−1–

Solute transport
Time weighting scheme

2D – simple
2D –axysymmetricalvertical ﬂow
1050
21 (at every
0.50 cm)
31 (at every
1.66 cm)
1.00E−05
3.00E−06
1.00E−03

1
Pea gravel
0.045
0.27
0.145
2.68
172,800
0.5
No temperaturedependence
Crank-Nicolson
schemeGalerkin ﬁnite
elements
Millington & Quirk
CWM1Space weighting scheme
Use tortuosity
Wetland module
Iteration criteria for nonlinear
adsorption
Absolute concentration tolerance
Relative concentration tolerance
Maximum number of iterations
Soil speciﬁc parameters
Bulk. D.
Disp. L.
Disp. T.
Solute speciﬁc parameters
Diffusion of oxygen in water
Diffusion of organics in water
Diffusion of nutrients in water
Diffusion of oxygen in gaseous
phase
Diffusion of organics in gaseous
phase
Diffusion of nutrients in gaseous
phase

Value ordescription,
respectively

vanGenuchten–Mualem,
w/o air-entry valueof −2 cm
No Hysteresis

Hysteresis
Material properties for water ﬂow
Number of materials
Name
Residual water, Qr
Porosity, Qs
Alphan
Saturated hydraulic conductivity, Ks
I
Temperature dependence

113

0.001
0.001
10
g cm−3
cm
cm

1.79
500
400

cm2 d−1
cm2 d−1
cm2 d−1
cm2 d−1

1.73041.0944
1.922418456

cm2 d−10

cm2 d−10

Adsorption parameters – ammonium and ammonia nitrogen (NH4 )Kd
cm3 ␮g−1
cm3 ␮g−1
Nu
–
Beta
d−1
Alpha

0.9
0.003
1.2
2

Adsorption parameters – slowly biodegradable particulate COD
Kd
cm3 ␮g−1
–
Beta
d−1
Alpha

0.575
1
2.4

Table A.3
Constructed Wetland Model No.1 parameters.
Unit

Parameter

Hydrolysis
Hydrolysis rate constant

d−1

0.58

Autotrophic bacteria
Max. growth rate XA
Sat./Inh. coef. for SO

d−1
mg O2 L−1

1.5
0.0026

Fermenting bacteria
Max. growth rate XFB

d−1

3.77

Acetotrophic methanogenic bacteria
Sat./Inh. coef. for SO
Sat./Inh. coef. for SNO

mg O2 L−1
mg N L−1

0.02
0.05

Acetotrophic sulphate reducing bacteria
d−1Max. Growth Rate XASRB
Sat./Inh. coef. for SO
mg O2 L−1
mg N L−1
Sat./Inh. coef. for SNO
Sat./Inh. coef. for SSO4
mg S L−1
Yield coefﬁcient for XAMBYield coefﬁcient for XASRB

0.5
0.05
0.5
5
0.04
0.04

Composition parameters
N Content of XS
N Content of XIg N g CODXS −1
g N g CODXI −1

0.079
0

OxygenRate O2

d−1

0

114

T.G. Pálfy, G. Langergraber / Ecological Engineering 68 (2014) 105–115

Table A.4
Root model and atmospheric boundary condition parameters for different temperatures.
12 ◦ C

Root water and solute uptake

Dissolved oxygen (SO)␮g cm−3cRootAmmonium and ammonia nitrogen (SNH)cRoot␮g cm−3
Water uptake reduction modelCritical stress index
Solute stress model
S-shape parameters
P50
cm
P3
–
PW
cm

16 ◦ C20 ◦ C24 ◦ C

U.

C.

T.

Sch.

U.

C.

T.

Sch.

U.

C.

T.

Sch.

U.

C.

T.

Sch.

−4.2

−4.9

−5.6

−5.1

−27.6

−4.4

−3.4

−3.0

−3.9

−0.8

−1.2

−0.7

−3.8

−0.9

−0.8

−0.8

4.7
S-Shaped1
No solute stress−800
3−1.00E+10
12 ◦ C

Atmospheric boundary condition

16 ◦ C20 ◦ C−1

Precip.
Evap.
Transp.
hCritA
Surface area associated with transpiration

cm d
cm d−1cm d−1
cm
cm2

1.66

0
0
−0.01
314

1.85

24 ◦ C

2.04

2.23

U. = unplanted, C. = Carex, T. = Typha, Sch. = Schoenoplectus.

Table A.5
Initial concentrations.
12 ◦ C

Initial concentrationsU.

16 ◦ C

C.

T.

Sch.

U.

20 ◦ C

C.

T.

−3

Sch.

U.

24 ◦ C

C.

T.

Sch.

U.

C.

SO
SF
SA
SI

␮g cm
␮g cm−3␮g cm−3␮g cm−3

113.0
45.2
18.1

116.9
46.7
18.7

111.8
44.7
17.9

120.8
48.3
19.3

SNH (in liquid
phase)
SNH (in solid
phase)SNO

␮g cm−3

18.2

18.7

18.6

16.4

SSO4
SH2 S

␮g cm−3␮g cm−3

13.6

XS (in liquid
phase)XS (in solid
phase)XI

␮g cm−3

275.7

␮g g−1

XH
XA
XFB
XAMB
XASRB
XSOB

␮g g−1
␮g g−1
␮g g−1
␮g g−1
␮g g−1
␮g g−1

␮g g−1

T.

Sch.

1

32.7

1.0

2.8

2.4

38.5

1.5

3.0

−3

␮g cm

3.5

40.0

1.8

5.0

6.0

38.5

2.0

5.0

4.2

0.69
0.16
0.76
0.93
0.77
0.04

1.03
0.15
0.84
1.66
0.77
0.16

0.0
13.0

16.0

12.8

272.7

294.6

0.0
285.1

␮g g−1

0.0
0.0
5.48
1.10
0.02
0.89
0.90
0.04

4.47
0.40
0.69
3.05
0.46
0.37

2.24
0.80
0.08
0.77
0.84
0.12

4.49
0.81
1.53
2.52
0.99
0.15

0.20
1.20
0.84
0.88
1.01
0.03

5.00
0.60
1.70
2.80
1.01
0.10

0.12
1.20
0.84
1.40
1.12
0.07

2.08
1.20
0.84
1.98
0.80
0.12

1.26
0.57
0.56
1.08
1.06
0.04

0.68
0.58
1.92
2.04
1.36
0.07

1.30
0.63
0.66
1.55
1.01
0.04

0.73
0.58
0.84
1.42
1.01
0.25

5.48
1.01
0.76
1.44
2.17
0.04

0.68
0.63
0.76
1.44
1.00
0.10

U. = unplanted, C. = Carex, T. = Typha, Sch. = Schoenoplectus.

References
Ahnert, M., Blumensaat, F., Langergraber, G., Alex, J., Woerner, D., Frehmann,T., Halft, N., Hobus, I., Plattes, M., Spering, V., Winkler, S., 2007. Goodnessof-ﬁt measures for numerical modelling in urban water management –a review to support practical applications. In: IWA (Ed.), Proceedings ofthe 10th IWA Specialised Conference on “Design, Operation and Economicsof Large Wastewater Treatment Plants”, 9–13 September. Vienna, Austria,
pp. 69–72.
Allen, W.C., Hook, P.B., Biederman, J.A., Stein, O.R., 2002. Temperature and wetland plant species effects on wastewater treatment and root zone oxidation. J.
Environ. Qual. 31, 1010–1016.
Campbell, C.S., Ogden, M., 1999. Constructed Wetlands in the Sustainable Landscape.
Wiley & Sons, New York, pp. 45, 48–50, 54–55.
Henze, M., Gujer, W., Mino, T., Loosdrecht, M. (Eds.), 2000. Activated sludge models
ASM1, ASM2, ASM2d and ASM3. IWA Scientiﬁc and Technical Report No. 9. IWA
Publishing, London, UK.

Kadlec, R.H., Wallace, S.D., 2009. Treatment Wetlands, 2nd ed. CRC Press, Boca Raton,FL.
Langergraber, G., (PhD thesis) 2001. Development of a simulation tool for subsurface
ﬂow constructed wetlands. Institute of Sanitary Engineering and Water Pollution Control, University of Natural Resources and Life Sciences (BOKU), Vienna,
Austria.
Langergraber, G., 2008. Modeling of processes in subsurface ﬂow constructed wetlands – a review. Vadoze Zone J. 7 (2), 830–842.
Langergraber, G., 2011. Numerical modelling: a tool for better constructed wetland
design? Water Sci. Technol. 64 (1), 14–21.
Langergraber, G., Šimůnek, J., 2011. HYDRUS Wetland Module Manual, Version 2.
Hydrus Software Series 4. Department of Environmental Sciences, University ofCalifornia Riverside, Riverside, CA, 46 pp.
Langergraber, G., Šimůnek, J., 2012. Reactive transport modelling of subsurface ﬂow
constructed wetlands using the HYDRUS wetland module. Vadose Zone J. 11(2), http://dx.doi.org/10.2136/vzj2011.0104 (special issue “Reactive Transport
Modelling”).

T.G. Pálfy, G. Langergraber / Ecological Engineering 68 (2014) 105–115
Langergraber, G., Rousseau, D., García, J., Mena, J., 2009. CWM1 – a general model todescribe biokinetic processes in subsurface ﬂow constructed wetlands. Water
Sci. Technol. 59 (9), 1687–1697.
Mburu, N., 2013. Personal communication.
Mburu, N., Sanchez-Ramos, D., Rousseau, D.P.L., van Bruggen, J.J.A., Thumbi, G., Stein,O.R., Hook, P.B., Lens, P.N.L., 2012. Simulation of carbon, nitrogen and sulphur
conversion in a batch-operated experimental wetland mesocosm. Ecol. Eng. 42,
304–315.
Pálfy, T.G., (Master thesis) 2013. Veriﬁcation of the implementation of CWM1 in the
HYDRUS Wetland Module. Institute of Sanitary Engineering and Water Pollution Control, University of Natural Resources and Life Sciences (BOKU), Vienna,
Austria, 66 pp.

115

Reichert, P., 1998. AQUASIM 2.0 – User Manual. Swiss Federal Institute for Environmental Science and Technology (EAWAG), Dübendorf, Switzerland, ISBN:3-906484-16-5.
Salvato, J.E., Nemerow, N.L., Agardy, F.J. (Eds.), 2003. Environmental Engineering.
Wiley and Sons, New Jersey, p. 269.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2011. The HYDRUS software package
for simulating the two- and three-dimensional movement of water, heat, andmultiple solutes in variably-saturated media. Technical Manual, version 2.0. PC
Progress, Prague, Czech Republic.