Ecological Engineering 68 (2014) 209–213

Contents lists available at ScienceDirect

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

Short communication

Modelling the response of laboratory horizontal ﬂow constructedwetlands to unsteady organic loads with HYDRUS-CWM1
A. Rizzo a,∗ , G. Langergraber b , A. Galvão c , F. Boano a , R. Revelli a , L. Ridolﬁ aaDepartment of Environment, Land and Infrastructure Engineering (DIATI), Politecnico di Torino, Turin, Italy
Institute of Sanitary Engineering and Water Pollution Control, University of Natural Resources and Applied Life Sciences, Vienna, AustriacInstituto Superior Técnico, Technical University of Lisbon, Lisbon, Portugalb

a r t i c l e

i n f oArticle history:
Received 22 August 2013
Received in revised form 6 February 2014
Accepted 31 March 2014
Keywords:
Constructed wetland
Horizontal subsurface ﬂow
Unsteady organic loads
Process-based model
HYDRUSCWM1

a b s t r a c t
Although horizontal ﬂow constructed wetlands (HF CWs) are usually subjected to unsteady loads in realapplication, the modelling of HF CW response to time-variable loads has been scarcely studied in literatureyet. The aim of this study is to test the capability of HYDRUS-CWM1 to simulate the behavior of HF CWs
subjected to unsteady loads. Hence, we applied HYDRUS-CWM1 to simulate laboratory results of HF CWs
subjected to variable COD inﬂows. The modelling results adequately ﬁt the experimental data, with analmost perfect agreement in global COD removal efﬁciencies (67 and 68% from laboratory experiments
and simulations, respectively), and mean percent error equal to 20 and 31% for efﬂuent COD and NH+
4concentrations, respectively. The obtained results suggest that HYDRUS-CWM1 can be a powerful tool tosimulate the response of HF CWs under time-variable loads. Additionally, more detailed data are shown
to be crucial in order to better exploit process-based model tools.
© 2014 Elsevier B.V. All rights reserved.

1. Introduction
Constructed wetlands (CWs) are a wastewater treatment
technology that exploits the physical, chemical, and biological processes occurring in soils to improve water quality (Kadlec andWallace, 2009). The low operation and maintenance costs allow
for a widespread use of CWs, especially for treating wastewaterfrom small communities (Haberl et al., 2003; Kadlec and Wallace,2009).
Horizontal subsurface ﬂow (HF) is a common CW type to treat
organic loads (Kadlec and Wallace, 2009; Abou-Elela et al., 2013).
In HF CWs, wastewater ﬂows horizontally beneath the soil surface
and from the inlet to the outlet of a gravel bed planted with wetland vegetation (Haberl et al., 2003; Kadlec and Wallace, 2009). Dueto the reductive state developed, HF CWs are usually preferred toother CW types when anaerobic processes for treating wastewater
are sufﬁcient (Kadlec and Wallace, 2009). Moreover, HF CWs arealso preferred to surface ﬂow CWs for concerns about human contact with untreated wastewater, mosquito and odor control, and∗ Corresponding author. Tel.: +39 0110905684.
E-mail addresses: anacleto.rizzo@polito.it (A. Rizzo),
guenter.langergraber@boku.ac.at (G. Langergraber), fulvio.boano@polito.it(F. Boano), roberto.revelli@polito.it (R. Revelli), luca.ridolﬁ@polito.it (L. Ridolﬁ).
http://dx.doi.org/10.1016/j.ecoleng.2014.03.073
0925-8574/© 2014 Elsevier B.V. All rights reserved.

minimization of wildlife interactions within the wetland (Kadlec
and Wallace, 2009).
Process-based models are powerful tool to investigate thebehavior of HF CWs, since they are able to separately simulate theeffect of the different organic removal processes (Langergraber,2008). However, process-based models were typically used toinvestigate HF CWs subjected to steady state loads (e.g., Llorens
et al., 2011; Samso and Garcia, 2013a) and only recently they werestarted to be used for unsteady long-term simulations (e.g., Samsoand Garcia, 2013b). Generally, little is known about the ability ofprocess-based model to simulate the performance of HF CWs under
unsteady loads.
The aim of the work is to test the capability of HYDRUS-CWM1(Langergraber and Šimůnek, 20012) to model the HF CW response
to time-varying loads simulating the laboratory data collected byGalvão and Matos (2012).

2. Methods
In this section the main features of experimental dataset andof the used process-based model are brieﬂy exposed. For moredetails about data from laboratory HF CWs see Galvão and Matos(2012), while about HYDRUS-CWM1 see Langergraber et al. (2009),
Šimůnek et al. (2011a), and Langergraber and Šimůnek (2012).

210

A. Rizzo et al. / Ecological Engineering 68 (2014) 209–213

CODin [mg/l]

1500

1000

500

0

0

20

40

60
t [d]

80

100

130

Fig. 1. Unsteady COD inﬂow concentrations. The squares and circles represent themeasured data from Galvão and Matos (2012) for group A and B, respectively. The
continuous and dashed lines ﬁgure the inﬂow trend assumed in simulations for A
and B, respectively.

2.1. Experimental data
The data used to test the HYDRUS-CWM1 model were collectedfrom laboratory experiments on six HF CWs set up at the Technical University of Lisbon in 2010 (Galvão and Matos, 2012). EachHF CW consisted of a gravel bed with length 1.1 m, width 0.76 m,
and height 0.3 m. HF CWs were fed with an average hydraulic loadof 10 l/d as follows: from Tuesday to Thursday 10 l/d; on Mondays
and Fridays 20 l/d, to mimic the weekends (no feeding during Saturday and Sunday). A concentrated synthetic wastewater (i.e., nobacteria and no solids) was used with theoretical concentrationsof COD and total nitrogen (TN) equal to 39 and 6 g/l, respectively.
The synthetic wasterwater was composed with a mixture of urea,
acetate, peptone, starch, powdered milk, and soy oil (for concentrations see Galvão and Matos (2012)). This solution was diluted toobtain the desired concentration to feed each bed. After an inoculation period, the six beds were divided in two groups, A and B, fed
for 127 days with higher (A) and lower (B) average COD concentrations, respectively, as shown in Fig. 1. Each group was composedof three beds with different vegetation conﬁgurations: no plants, P.
australis, and Scirpus. Water level was constantly maintained 5 cm
below the surface. Since the HF CWs were situated inside laboratory, the resulting evaporation rates were very low (see Galvão andMatos (2012)) and are neglected in the present analysis.
The results from Galvão and Matos (2012) have shown no signiﬁcant role of plants, therefore the data from the three different
beds per group are averaged in the present work. In this way, wewith corresponding
obtain a single efﬂuent value of COD and NH+
4standard deviation per group, which is compared with results frommodelling simulations.
2.2. HYDRUS-CWM1
The HYDRUS Wetland Module is an extension of the HYDRUS
family codes, which allows for a detailed simulation of subsurface ﬂow CWs (Langergraber and Šimůnek, 20012). HYDRUS solves
systems of partial differential equations in three dimensions inorder to simulate (Šimůnek et al., 2011a): (i) variably-saturated
water ﬂow (Richards’ equation); (ii) transport of constituents (massbalance equations); (iii) inﬂuence of plants (water and nutrient
uptake, radial oxygen loss) and (iv) water temperature (heat transport equation). On this framework, the HYDRUS Wetland Module
adds two different biokinetic models: CW2D for aerobic and anoxic
processes, and CWM1 for aerobic, anoxic, and anaerobic transformations. Since HF CWs are principally anaerobic systems (Kadlec
and Wallace, 2009), CWM1 was adopted in our simulations. For
sake of simplicity, we refer to this option of HYDRUS Wetland
Module as HYDRUS-CWM1.

The biokinetic model CWM1 (Langergraber et al., 2009) simulates the transformation of organic matter (expressed as chemical
oxygen demand – COD), nitrogen (N), and sulfur (S) via biochemical transformation under aerobic, anoxic, and anaerobic conditions.
The organic matter is considered via both soluble and particulate COD fractions, while microorganisms are considered as ﬁxed
particulate components. The soluble and particulate groups arereferred as Sx and Xx , respectively, where x states the component
acronym. A list of the components included in CWM1 and of theseventeen biogeochemical transformations considered in CWM1
are reported in Supplementary online material.

2.3. Model set-up
The HF CWs studied by Galvão and Matos (2012) are modeled
via a vertical 2D rectangular domain 1.1 m long, and 0.30 m high.
In the soil hydraulic model, the van Genuchten–Mualem formulation is assumed for soil water retention and hydraulic conductivity(Šimůnek et al., 2011a). For the solute transport model, the Millington and Quirck formulation is set for soil tortuosity, while aLangmuir law is considered for NH+
adsorption (Šimůnek et al.,4
2011a).
Since the experiments of Galvão and Matos (2012) did not show
relevant difference between vegetated and unvegetated CWs, theplant effects (i.e. water and nutrient root uptake, radial oxygen loss)
are not considered. The scarce role of plant can be explained interm of evapotraspiration rate. Differently from outdoor locations
in which evapotranspiration plays an important role (Pedescoll
et al., 2013), very low evapotranspiration rates are reported for theindoor laboratory HF CWs studied by Galvão and Matos (2012). This
conﬁrms the scarce role of root water uptake in the hydrodynamics
of Galvão and Matos (2012)’s HF CWs. An additionally explanation can be related to root development. Indeed, Galvão and Matos(2012) have also mentioned no signiﬁcant removal rates insideeach group related to root compartment. This probably because
the plants were very young, so the root system was not expected
to be fully developed. For sake of simplicity, the inﬂuent O2 is set
equal to zero, in accordance with other studies dealing with modelling HF CWs (e.g., Samso and Garcia, 2013a,b). The sulfur cycle is
not modelled for lack of measured data.
Since the inﬂuent wastewater is synthetic, typical COD fractionations for domestic wastewater proposed in literature (e.g., Henzeet al., 2000) are not usable. Hence, the COD inﬂow concentration is
fractionated in accordance with inﬂuent load composition as follows: 62% for SF , 10% for SA , 3% for SI , 20% for XS , and 5% for XI .
Since NH+
inﬂow was not monitored, an average ratio between
4
and
COD
inﬂow, NH+/COD = 0.02, is used to estimate NH+
NH+
4
4
4
inﬂow concentration.
Initial conditions are obtained simulating 10 days of inoculation,with low initial component concentration and COD inﬂow concentrations equal to 1108 and 510 mg/l for group A and B, respectively.
The adopted boundary conditions have the following features:the inﬂow is modeled as distributed on the whole water depth;
outﬂow is located in a single node at 1 cm from the bottom andset as constant pressure head in order to maintain the water table
5 cm below the top; the upper boundary is not covered, therefore the atmospheric O2 is allowed to diffuse within CW; all otherboundaries are assumed impermeable and with no solute ﬂuxes.
A graphical representation of boundary conditions is reported in
Supplemental online material.
Only few parameters are modiﬁed from the default values
reported in Langergraber et al. (2009); their values are reported
together with other set parameters in Table 1. Note that the fractions of organic nitrogen are calibrated in order to guarantee theA. Rizzo et al. / Ecological Engineering 68 (2014) 209–213
Table 1
List of the parameters set to simulate experimental HF CWs of Galvão and Matos(2012).
CWM1 biokinetic model parameters
Rate of constant XFB lysis
Fraction of XI generated in biomass lysis
N content of SF
N content of SI
N content of XS
N content of XI

0.5a
0.01a
0.137b
0.117b
0.147b0.137bSoil hydraulic model parameters
Porosity
Saturated hydraulic conductivity
Residual soil water content (gravel–sand)Parameter ˛ in the soil water retention function(gravel–sand)Parameter n in the soil water retention function
Tortuosity parameter in the conductivity function(gravel–sand)

0.3c
72 m/da0.045d
14.5 m−1 d4a0.5dAdsorption isotherm model parameters for SNH
Langmuir coefﬁcient ks (gravel)
Langmuir coefﬁcient  (gravel)

0.411 m3 /ge
5.35 × 10−4 m3 /ge

Inﬂow composition parameteron COD ratio
Average inﬂow NH+
4

0.02a

a

Fit parameter.
b
Inﬂow total nitrogen balance guaranteed.
cGalvão and Matos (2012).
d
Šimůnek et al. (2011b).
e
Zhu et al. (2011).

211

results in term of global removal efﬁciency. Indeed, the average
global removal efﬁciencies of simulated COD efﬂuent (68 and 69%
for group A and B, respectively) result strongly similar to the values
reported in Galvão and Matos (2012)’s experiments (69 and 66% forgroup A and B, respectively), with a maximum discrepancy of only
3%. Also the simulated behavior of daily COD efﬂuent concentrations adequately
 follows the measured
 data, with a mean percent
error (MPE – (CODm − CODs )/CODm × 100 – subscript m and s
represent measured and simulated, respectively) equal to 22.0 and
17.5% for group A and B, respectively. The trend shows a weekly
pattern with sudden efﬂuent decreases driven by the absence ofloads during Saturday and Sunday.
efﬂuent concentration are reported in Fig. 2,
Simulated NH+
4(right panels), showing also in this case a good ﬁt with experimental data. MPE is satisfactorily low even if slightly higher than COD
efﬂuent (37.0 and 24.3% for group A and B, respectively), justiﬁedby higher uncertainty in NH+
inﬂow concentration driven by the
4/COD during the whole simulations.
need to assume a constant NH+
4
In order to test if the water temperature variations are responsible to differences between simulated results and experimental
data, additional simulations were run using the temperature data
recorded by Galvão and Matos (2012). The results indicate that
including temperature variations does not lead to appreciablechange in simulated efﬂuent concentrations from Galvão andMatos (2012)’s HF CWs (data not shown).
4. Discussion

total nitrogen balance in the inﬂow. The simulation are run assuming a constant temperature of 20 ◦ C.
The domain is discretized via triangular mesh elements withcatheti parallel to the HF CW length and height. The triangular meshelements are 1 cm high and 1–5 cm in length to have a ﬁner meshnear the inﬂow. The maximum time step is set equal to 0.001 d.
A modiﬁcation of the HYDRUS time step adjustment which takes
into account model component concentrations was implemented
in order to avoid numerical problems. More in detail, HYDRUS code
is now implemented with an automatically reduction of the time
step in case of a high error in total mass balance of any simulated
chemical compounds (e.g., when the O2 oxygen consumption is toofast). This new feature will be added into future HYDRUS version.
3. Results
The simulated COD efﬂuent concentrations are reported in Fig. 2(left panels). The simulations follow very well the experimental

The MPE on COD efﬂuent of 20% can be viewed as a good starting
point in the study of unsteady COD efﬂuent from HF CWs, since
up to now the validation of process-based model results is usuallyconcentrated only on the bulk system response (e.g., Llorens et al.,2011), or via a qualitative check with graphs similar to Fig. 2 (e.g.,Samso and Garcia, 2013b).
The good simulation results are obtained with a relatively simple manual calibration, changing only few parameters from thedefault values. A more reﬁned calibration (e.g., based on optimization algorithms) would not be very useful due to the lack of moredetailed input data, e.g. inﬂuent NH+
concentrations, tracer exper4
iments (e.g., Ranieri et al., 2013). Despite this lack of information,
HYDRUS-CWM1 shows a satisfactory ﬂexibility in calibration of HF
CWs with scarcely detailed data.
The great advantage of process-based models is the possibility to isolate the effect of single processes on the whole behaviorof the system. An example is given by the simulated COD efﬂuent
fractionation exposed in Fig. 3, in which acetate (SA ) results the

700

150
100

350

40

60

80

100

130

350

175

0
0

[mg/l]

20

4,out

0
0

NH

COD [mg/l]

50
0
0

20

40

60

80

100

20

40

60
t [d]

80

100

60

40

20

40

60 80
t [d]

100

130

20
0

Fig. 2. Comparison between COD (left panels) and NH+(right panels) efﬂuent concentrations simulated (continuous line) and measured (circle with standard deviation bar)4for group A (upper panels) and B (lower panels).

212

A. Rizzo et al. / Ecological Engineering 68 (2014) 209–213

500
COD

400

tot

300

COD

out

[mg/l]

200
100
0
0

(100 %)SF (2 %)SA (38.2 %)SI (8.78 %)XS (33.4 %)XI (17.6 %)
20

40

60

80

100

120

200
CODtot (100 %)150

SF (4.6 %)SA (37.8 %)SI (9.09 %)XS (29.7 %)XI (18.7 %)

100
50
0
0

20

40

60

80

100

120

t [d]
Fig. 3. Simulated fractionation of COD efﬂuent concentration for group A (upper panel) and B (lower panel): total COD (thick continuous line); fermentable, readily biodegradable soluble COD (SF – thin continuous line); fermentation products as acetate(SA – dashed line); inert soluble COD (SI – dash-dotted line); slowly biodegradable particulate
COD (XS – dotted line); inert particulate COD (XI – continuous line with circles). In the legend are reported the contribution of each component to the total efﬂuent COD
expressed in percentages.

most abundant component. This theoretical result is in accordancewith other modelling (Samso and Garcia, 2013b) and experimental(Huang et al., 2005) studies. Fig. 3 also shows that the highest COD
efﬂuent peak (around 20th and 30th days) is principally due to SA ,suggesting a lower ability of acetotrophic methanogenic bacteria(XAMB ) to degrade SA at suddenly high COD concentration. However, available experimental data are suitable to verify only thecapability of the model to simulate the global performance, i.e.
COD concentrations, and not detailed COD efﬂuent fractionation.
In the future, more accurate experiments (e.g., measurement ofmethane emissions (Corbella and Puigagut, 2013)) will be neededin order to test the ability of process-based model to simulate all
the included processes and not only the global response of CWs.
Generally, more detailed data would strengthen the usefulness ofprocess-based models, conﬁrming theoretical results and allowing
to build more detailed models.

Consequently, we would provide new insights to improve the HF
CW design procedure.

Acknowledgment
The Ph.D grant “Fondo Giovani MUR/09” of A. Rizzo is funded by
the the Italian Ministry of University and Research (MIUR).

Appendix A. Supplementary data
Supplementary data associated with this article can befound, in the online version, at http://dx.doi.org/10.1016/j.ecoleng.
2014.03.073.

References
5. Conclusions
HYDRUS-CWM1 models very well the global response of Galvão
and Matos (2012)’s HF CWs to variable loads, with an average
global removal efﬁciency of simulated COD efﬂuent (68%) similar
to the value registered in the experiments (67%). Moreover, modelling results show that HYDRUS-CWM1 adequately simulates also
HF CW response to unsteady loads at daily scale, with an acceptable error of 20% on COD efﬂuent concentrations. Even the dailyNH+
efﬂuent concentrations are satisfactorly captured under more
4
uncertain NH+
inﬂow information, with a MPE of 31%.
4
The results of this study suggest that HYDRUS-CWM1 can beused to study HF CWs under unsteady ﬂows, opening new possible future applications. Among these, the stochastic behavior of
HF CW outﬂows is of great interest. Indeed, the randomness of HF
CW efﬂuent is generally neglected during the typical design procedure, and the HF CW size should be highly overestimated to facereal outﬂow variability (Kadlec and Wallace, 2009). In this context, HYDRUS-CWM1 can be a powerful tool to better investigate
HF CW efﬂuent stochasticity. A future improvement could be themodelling of stochastic inﬂuent loads modeled as noisy signals(Ridolﬁ et al., 2011). In this way, we could investigate if the randomness of efﬂuent loads is driven by the variability of inﬂow loads.

Abou-Elela, S.I., Golinielli, G., Abou-Taleb, E.M., Hellal, M.S., 2013. Municipal wastewater treatment in horizontal and vertical ﬂows constructed wetlands. Ecol.
Eng. 61 (Part A), 460–468.
Corbella, C., Puigagut, J., 2013. Improving the reliability of closed chamber methodologies for methane emissions measurement in treatment wetlands. Water Sci.
Technol. 68 (9), 2097–2102.
Galvão, A., Matos, J., 2012. Response of horizontal sub-surface ﬂow constructedwetlands to sudden organic load changes. Ecol. Eng. 49, 123–129.
Haberl, R., Grego, S., Langergraber, G., Kadlec, R., Cicalini, A.-R., Dias, S., Novais, J.,
Aubert, S., Gerth, A., Thomas, H., Hebner, A., 2003. Constructed wetlands for thetreatment of organic pollutants. J. Soils Sediments 3 (2), 109–124.
Henze, M., Gujer, W., Mino, T., Van Loosdrecht, M., Henze, M., 2000. Activated Sludge
Models ASM1, ASM2, ASM2d and ASM3. IWA Publishing, Eighteenth century
collections online.
Huang, Y., Ortiz, L., Aguirre, P., Garcia, J., Mujeriego, R., Bayona, J., 2005. Effect ofdesign parameters in horizontal ﬂow constructed wetland on the behaviour ofvolatile fatty acids and volatile alkylsulﬁdes. Chemosphere 59 (6), 769–777.
Kadlec, R., Wallace, S., 2009. Treatment Wetlands, 2nd ed. CRC press, Boca Raton, FL.
Langergraber, G., 2008. Modeling of processes in subsurface ﬂow constructed wetlands: a review. Vadose Zone J. 7 (2), 830–842.
Langergraber, G., Šimůnek, J., 2012. Reactive transport modeling of subsurface ﬂow
constructed wetlands using the HYDRUS wetland module. Vadose Zone J. 11 (2),Special Issue “Reactive Transport Modeling”.
Langergraber, G., Rousseau, D.P.L., Garcia, J., Mena, J., 2009. CWM1: a general model
to describe biokinetic processes in subsurface ﬂow constructed wetlands. Water
Sci. Technol. 59 (9), 1687–1697.
Llorens, E., Saaltink, M.W., Poch, M., Garcia, J., 2011. Bacterial transformation andbiodegradation processes simulation in horizontal subsurface ﬂow constructedwetlands using CWM1-RETRASO. Bioresour. Technol. 102 (2), 928–936.

A. Rizzo et al. / Ecological Engineering 68 (2014) 209–213
Pedescoll, A., Sidrach-Cardona, R., Snchez, J., Bares, E., 2013. Evapotranspiration
affecting redox conditions in horizontal constructed wetlands under mediterranean climate: inﬂuence of plant species. Ecol. Eng. 58, 335–343.
Ranieri, E., Gorgoglione, A., Solimeno, A., 2013. A comparison between model andexperimental hydraulic performances in a pilot-scale horizontal subsurface ﬂow
constructed wetland. Ecol. Eng. 60, 45–49.
Ridolﬁ, L., D’Odorico, P., Laio, F., 2011. Noise-Induced Phenomena in the Environmental Sciences. Cambridge University Press, Cambridge.
Samso, R., Garcia, J., 2013a. Bacteria distribution and dynamics in constructed wetlands based on modelling results. Sci. Total Environ. 461, 430–440.
Samso, R., Garcia, J., 2013b. BIO PORE, a mathematical model to simulate bioﬁlm
growth and water quality improvement in porous media: application and calibration for constructed wetlands. Ecol. Eng. 54, 116–127.

213

Šimůnek, J., van Genuchten, M.T., Šejna, M., 2011a. The HYDRUS Software Package forSimulating Two- and Three-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media, Technical Manual, Version 2.0. PC
Progress, Prague, Czech Republic, 258 pp.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2011b. The HYDRUS Software Package
for Simulating Two- and Three-Dimensional Movement of Water, Heat, andMultiple Solutes in Variably-Saturated Media, User Manual, Version 2.0. PC
Progress, Prague, Czech Republic, 286 pp.
Zhu, W.-L., Cui, L.-H., Ouyang, Y., Long, C.-F., Tang, X.-D., 2011. Kinetic adsorptionof ammonium nitrogen by substrate materials for constructed wetlands. Pedosphere 21 (4), 454–463.