FITOVERT: A dynamic numerical model of subsurface vertical flow constructed wetlands

of 22 P 612 9 March 2009 Accepted 18 May 2009 Available online 16 June 2009 Keywords: Program languages: Matlab 6.0� language Program size: 1480 KB Availability and costs: contact developers even distribution of fed wastewater, oxygen diffusion, adequate retention time and limited clogging at the same time. Further research is still needed to optimize this technology and to make it completely reliable. Numerical models can help to obtain a better understanding of the processes occurring in CWs for the improvement of existing design and operation criteria (Langer- graber, 2007). Moreover, some novel configurations of VSSF-CW have been recently proposed and only partially investigated (Babatunde et al., * Corresponding author. Tel.: þ39 0502217718; fax: þ39 0502217730. E-mail addresses: [email protected] (D. Giraldi), [email protected] Contents lists availab Environmental Mod journal homepage: www.els Environmental Modelling & Software 25 (2010) 633–640 (R. Iannelli). Name of the software: FITOVERT – Version 0.1 Developers: R. Iannelli, D. Giraldi, M. de Michieli Vitturi Contact address: R. Iannelli (first contact for software request): University of Pisa – Department of Civil Engineering, via C.F. Gabba 22 – 56122 Pisa (PI), Italy, Tel.: þ39 050 2217.718 fax: þ39 050 2217.730, e-mail: r.iannelli@ing. unipi.it Year first available: 2007 PC hardware required: the same as for Matlab 6.0� Software required: Matlab 6.0� Different types of constructed wetlands (CWs) are currently used to treat wastewater; basically they can be classified, according to the path of wastewater flux, in: surface flow, horizontal subsurface flow, and vertical subsurface flow (VSSF-CW) systems (Kadlec and Knight, 1996; USEPA, 2000; Kadlec et al., 2000). Increasingly VSSF-CWs are appreciated for their remarkable capacity to oxidize both organic matter and nitrogen compounds. However, their operation is more complex than the other CW configurations: usually they are intermittently fed and sometimes they present several vertically stacked layers in order to provide Constructed wetlands Hydrodynamics Modelling Reactive transport Vertical subsurface flow Unsaturated flow Software availability 1364-8152/$ – see front matter � 2009 Elsevier Ltd. doi:10.1016/j.envsoft.2009.05.007 opment of FITOVERT was to keep the complexity of the model to an acceptable level, so as to provide a practical tool for design and operation optimization. The dynamic formulation of the model allows to simulate the typical non stationary feeding-emptying operation of VSSF-CWs. FITOVERT is able to describe the water flow through porous media in unsaturated conditions, combined with evapotrans- piration; its biochemical module describes the degradation of both organic matter and nitrogen; the transport in the liquid phase is implemented for both dissolved and particulate components; the oxygen transport in the gaseous phase of the soil and its exchange with the liquid phase are also considered. As a main advantage, compared to the few currently available dedicated numerical models, FITOVERT is able to handle the porosity reduction due to bacteria growth and accumulation of particulate components, so that the clogging process is also simulated as an effect of the pore size reduction on the hydraulic conductivity of the simulated system. The performance of the model was firstly analyzed by comparison with hydrodynamic tests recorded in an experimental VSSF-CW pilot plant: tracer test were carried out in three different saturation conditions (fully saturated, partially saturated, and completely drained). FITOVERT proved to accurately simulate the hydraulic behaviour of VSSF-CWs in both saturated and unsaturated conditions. The needs for model improvements and further calibration are finally discussed. � 2009 Elsevier Ltd. All rights reserved. 1. Introduction Received 14 October 2008 Received in revised form This paper introduces a mathematical model (FITOVERT) specifically developed to simulate the behav- iour of vertical subsurface flow constructed wetlands (VSSF-CWs). One of the main goals of the devel- Article history: FITOVERT: A dynamic numerical model flow constructed wetlands D. Giraldi a, M. de Michieli Vitturi b, R. Iannelli a,* aDepartment of Civil Engineering, University of Pisa, via Carlo Francesco Gabba 22, 561 b Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Pisa, Via della Faggiola, 32, 5 a r t i c l e i n f o a b s t r a c t All rights reserved. subsurface vertical isa, Italy 6 Pisa, Italy le at ScienceDirect elling & Software evier .com/locate/envsoft Developed in the Matlab� environment, FITOVERT is completed with a user-friendly interface for quick access to the main input parameters and for a clear graphic display of the simulation results. Since numerous and complex physical, chemical, and biolog- ical processes are involved in VSSF-CWs, a preliminary selection of the most relevant of them was carried out basing on an extended literature analysis. The selected processes were then implemented in the model by way of a number of dedicated interlinked modules. The single modules of FITOVERT are hereafter briefly described. Fur further details one can refer to Iannelli et al. (2008). 2.1. Topological configuration FITOVERT has amono-dimensional vertical schemewith a freely selectable number of horizontally homogeneous layers (Fig. 1). The total depth of the bed and the main characteristics of each layer (thickness and hydraulic parameters) can be defined by the user. Each layer of the bed is discretized by stacked numerical elements whose number and thickness (equal for all the elements) can be defined as well. All the equations in FITOVERT refer to single numerical elements connected each other by the transport equations. A wastewater flux runs through the system from the top to the bottom of the bed. Horizontal flow is neglected. delling & Software 25 (2010) 633–640 2008); among others we mention ‘‘Tidal Flow’’ (Green et al., 1997; Sun et al., 1999, 2005, 2006; Zhao et al., 2004a; Austin, 2006; Lee and Scholz, 2006), ‘‘Anti-sized’’ (Zhao et al., 2004b) and ‘‘Recipro- cating’’ (Behrends et al., 2001) CW systems. Numerical models can act as an efficient platform to test these new configurations and to compare themwith the traditional ones, thus reducing the required efforts for experimental studies and evaluations. In the past, only simple models were developed, ranging from regressions (black box stochastic models) to deterministic models based on first-order or Monod-type equations (e.g. Rousseau et al., 2004). Simple models aim to offer basic tools for the design of CWs, but they provide only a limited understanding of the system: optimization of the facility and insight into the treatment process are not main objectives for this class of models. The mechanistic approach for modelling CWs has been adopted only recently, since it requires a significant effort for numerical implementation. Relatively few numerical models specifically developed to simulate CW systems have been reported (Brovelli et al., 2007). Most of them can simulate horizontal subsurface flow CWs. The numerical modelling of VSSF-CWs requires the implementation of reactive transport models for variably saturated conditions that present a higher complexity compared to models for saturated conditions (Langergraber et al., 2009). In his recent review, Lan- gergraber (2008) reported only four models that are able to simulate also VSSF-CWs, besides the most popular CW2D/HYDRUS (Langergraber and Simunek, 2005). However, even high complexity models can be inadequate in fully forecasting the behaviour of CWs and in providing reliable estimations of their parameters (Marsili- Libelli and Checchi, 2005). In this study we present a new numerical model (FITOVERT) specifically developed to simulate VSSF-CWs. The complexity of the model was maintained as far as possible to an acceptable level for its effective calibration and practical use. After a brief introduction of the code, the single modules are presented and discussed. A summary of the numerical imple- mentation of the complex system of equations the model is based on is also provided. The model has been firstly analyzed in comparison with hydrodynamic tests recorded in an experimental VSSF-CW: the results are reported and the efficiency of the hydraulic module is discussed. Finally we discuss the needs for the calibration of the whole model and some possible improvements of the code. 2. Model description The code is a dynamic multi-component reactive transport model for variably saturated conditions with volume filtration of particulates. It was expressly designed to simulate subsurface vertical flow CWs and therefore is referred as FITOVERT, a short- hand merge of the Italian words standing for CWs (FITOdepur- azione) with vertical (VERTicale) flow. The dynamic formulation of the model allows to simulate the typical non stationary behaviour of VSSF-CWs. FITOVERT is able to describe the water flow through porous media in unsaturated conditions combined with evapo- transpiration; its biochemical module describes the degradation of both organic matter and nitrogen; the transport in the liquid phase is implemented for both dissolved and particulate components; the oxygen transport in gaseous phase and its exchange with the liquid phase are considered as well. As a main advantage, compared to the few currently available models, FITOVERT is able to handle the porosity reduction due to bacteria growth and accumulation of particulate components, so that the clogging process is also simulated as an effect of pore size D. Giraldi et al. / Environmental Mo634 reduction on the hydraulic conductivity of the simulated system. Boundary conditions for the inlet allow to define time-varying hydraulic and organic loads. If the system becomes saturated, FITOVERT is able to automatically handle the ponding on the surface of the bed by properly changing the hydraulic boundary conditions. Boundary conditions for the outlet allow choosing between closed and open outlet (valve function). For the open outlet configuration it is possible to define the hydraulic head on the outlet. Boundary conditions for the outlet also allow defining aeration from the bottom of the bed to simulate aeration pipes. output input numerical elementlayer valve hydraulic head on the outlet Fig. 1. Topological scheme of FITOVER. for this difference by including the diffusion of dissolved compo- Table 1 Pilot VSSF-CW: layers and van Genuchten–Mualem parameters (layers are listed from the top to the bottom of the bed) (adapted from Giraldi et al., 2009). Layer Thickness, cm Particle size, mm qs qr a, cm �1 m Ks, cm s �1 1st 8 20–30 0.480 0.035 0.145 0.627 2 2nd 9 2–3 0.447 0.070 0.145 0.627 0.169 3rd 8 5–10 0.483 0.036 0.145 0.627 0.398 4th 8 8 0.456 0.035 0.145 0.627 0.495 5th 15 3–6 0.442 0.050 0.145 0.627 0.390 0 D. Giraldi et al. / Environmental Modelling & Software 25 (2010) 633–640 635 2.2. Hydraulic module The vertical water flow through porous media in unsaturated conditions isdescribedusing theRichardsequation (e.g. Raats, 2001): vq vt ¼ þ v vz � K � vh vz � 1 �� (1) where q [–] is the volumetric water content, t [T] is the time, z [L] is the spatial coordinate (positive downwards), K [LT�1] is the unsaturated hydraulic conductivity, h [L] is the matrix potential. The constitutive relationships between pressure head, hydraulic conductivity, and water content are handled using van Genuchten– Mualem functions and the related empirical parameters a [L�1] and m [–] for unsaturated conditions (van Genuchten, 1980): q ¼ qr þ qs � qr� 1þ jahj 11�m �m; K ¼ Ks � q� qr qs � qr �0:5241� 1� � q� qr qs � qr �1 m !m352 (2) where qr [–] and qs [–] are the residual and saturated volumetric water content, respectively, and Ks [LT �1] is the saturated hydraulic conductivity. Typical values for van Genuchten–Mualem parameters are available in literature for standard soil textures (sand, loam, silt, clay) and their combinations (e.g. Meyer et al., 1997). However, in VSSF-CWs, coarse gravel can be also found as medium material; moreover plant roots and biomass can modify the hydraulic parameters of the original medium. Therefore, dedicated analyses are often required to calibrate these parameters. We set up and carried out experimental procedures to define van Genuchten– Mualem parameters for the evaluation of the hydraulic module on a VSSF-CW pilot plant (see next section). The hydraulic module can run as a stand-alone model with significant saving of computation time if only hydraulic simulations are required. 2.3. Biochemical module The treatment processes for biodegradable organic matter and 6th 20 50–60 0.48 nitrogen compounds are described following the Activated Sludge Model 1 (ASM1) introduced by Henze et al. (1987). Thirteen components are taken into account by this model, seven of which Table 2 Operation conditions for tracer tests and results of the model analysis (adapted from Gir Test label Saturation condition Hydraulic load, m A Partial 1047 B Complete 2094 C Complete 1047 D Draining 1047 E Draining 1047 nents through the bio-film: the kinetics of diffusion can limit the global reaction velocity and therefore the overall treatment performance of the biological system. The kinetics of diffusion can become a limiting factor especially in presence of thick bio-films, while for thin bio-films it can be neglected. Since in VSSF-CWs the bio-film is supposed to be always thin due to the small size of the supporting medium, in FITOVERT the bio-film diffusion is consid- ered by a proper calibration of the kinetic parameters, with particular reference to the half saturation coefficients. The uptake of nutrients and metals by the plants is neglected in the current version of FITOVERT. This is not completely true in the short period; however the expected maximum removal rates of nitrogen, phosphorus and metals by direct plant uptake and har- vesting are small compared to typical loadings (e.g. USEPA, 2000). Adsorption is not considered as well: persistent pollutants and phosphorus removal are not simulated in this release. The modular coding of the model allows the substitution of ASM1 scheme with the updated ASM2 or ASM3 versions, able to simulate also the fate of phosphorus, and the implementation of saturation absorption equations as well; however the significantly increased number of calibration parameters suggested to delay these implementations to future releases of the model. 2.4. Transport of dissolved components The advection–dispersion transport in the liquid phase for soluble components is described according to Bresler’s equation (Bresler, 1973): q$ vc vt ¼ v vz � d$ vc vz � � q$vc vz þ R$q (3) are soluble (dissolved components) and six in suspension (partic- ulate components). The kinetics of biological degradation is based on a system of eight Monod-type equations with a wide set of stoichiometric, bio-kinetic, and mass balance parameters, for the calibration of which an extended set of literature work is available. ASM1 is a standard reference for modelling activated sludge (i.e. suspended growth biomass), while in VSSF-CWs, the active biomass is mainly attached to the rhizosphere and to the granular material forming the medium. Theoretically, one should account 0.035 0.145 0.627 2 which slightly differs from the other commonly used formulation (e.g. van Genuchten, 1982) in the definition of the dispersion coefficient. In Eq. (3), c [ML�3] is the concentration of single soluble aldi et al., 2009). m/d Dispersivity, cm Efficiency of the model 10 0.990 4.5 0.979 4.5 0.989 14 0.972 14 0.986 0 10 20 30 40 50 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 Time (hours) Time (hours) Time (hours) Time (hours) Tr ac er c on ce nt ra tio n (m icr og ra ms /l) 0 10 20 30 40 50 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 0 10 20 30 40 50 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Tr ac er c on ce nt ra tio n (m icr og ra ms /l) Tr ac er c on ce nt ra tio n (m icr og ra ms /l) Tr ac er c on ce nt ra tio n (m icr og ra ms /l) D A C E 0 10 20 30 40 50 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 satu D. Giraldi et al. / Environmental Modelling & Software 25 (2010) 633–640636 component in water, d [L2T�1] is the dispersion coefficient, q [L3L�2T�1] is the specific flow rate, and R [ML�3T�1] is the reaction term. The dispersion coefficient theoretically accounts both for the diffusion and the mechanical dispersion, but for the liquid phase the effect of diffusion compared to dispersion can be neglected. Under saturated and steady flow conditions, mechanical dispersion may be considered as proportional to average flow velocity by means of the dispersivity l [L] (Ogata, 1970). d ¼ l ���q q ��� (4) Fig. 2. Experimental data (dots) and numerical simulation (line) for Test A (partial Following Bresler (1973), we also assumed that the same rela- tionship is valid for unsaturated transient conditions. The reaction 0 5 10 15 20 25 266 266.5 267 267.5 268 268.5 269 269.5 270 Depth [cm] O 2 ga s [ M (O 2)L - 3 ] Depth [cm] 0 5 10 15 20 25 3.5 3.6 3.7 3.8 3.9 4 4.1 4.2 4.3 SN H [M (N )L - 3 ] Fig. 3. Time variations of gaseous oxygen (O2 gas), dissolved oxygen (SO), ammonia nitr a filling–empting cycle of a CW plant. term is determined by the biochemical module, as well as the diffusive exchange for dissolved oxygen. 2.5. Transport of particulate components and clogging process The transport and filtration of particulate components are described with an original scheme based on the work of Iwasaki (1937), later developed by Ives (1969) for the numerical analysis of the sand filtration process in saturated condition: ration), C (complete saturation), D, E (draining) (adapted from Giraldi et al., 2009). �q vTm vz ¼ þq vTm vt þ q$f $Tm (5) 0 5 10 15 20 25 2 3 4 5 6 7 8 9 Depth [cm] Depth [cm] SO [M (C OD )L - 3 ] day 28 time 00:00:00 day 28 time 00:25:54 day 28 time 00:51:55 day 28 time 01:17:41 day 28 time 01:43:27 day 28 time 02:08:45 day 28 time 02:35:35 day 28 time 03:00:00 0 5 10 15 20 25 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 SN O [M (N )L - 3 ] ogen (SNH) and nitrite–nitrate nitrogen (SNO) during a demonstrative simulation of where Tm [ML ] is the concentration of each single particulate component in water, and f [L�1] is the filter coefficient. the t delli FITOVERT is able to handle the porosity reduction due to bacteria growth and filtration of particulate components, by mean of a continuous update of the total volumetric specific deposit Dvtot [L3L�3]. The effect of pore size reduction on the hydraulic conductivity is also considered by the Carman–Kozeny’s equation (Boller and Kavanaugh, 1995), modified according to O’Melia and Ali (1978). The final formulation leads to define a correction term for the hydraulic conductivity: K ¼ K0h� 1þ p Dvtot30 �x� 1� Dvtot30 �yi (6) where K0 [LT �1] is the hydraulic conductivity for the clean filter, 30 [L3L�3] is the porosity of the clean filter, and p [–], x [–], and y [–] are empirical parameters. The transport of particulate components, jointly with the bio- logical module applied to the six considered particulate compo- nents, completely govern the fate of the biomass present in the bed in terms of biological activity, bed composition, dynamically time- varying porosity and hydraulic conductivity, allowing to simulate also the time-evolution of the clogging process. This is an innova- tive aspect of this model, since while new approaches are being developed to assess the build-up of bio-fouling in porous media (Brovelli et al., 2009), no other available specific CW models currently include this estimation. �3 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0.46 0.462 0.464 0.466 0.468 0.47 0.472 0.474 0.476 0.478 0.48 Depth [cm] Ep sil on day 0 time 00:00:00 day 4 time 06:51:51 day 8 time 13:43:31 day 12 time 20:35:11 day 17 time 03:26:51 day 21 time 10:18:31 day 25 time 17:10:11 day 30 time 00:00:00 Fig. 4. Porosity (left) and hydraulic conductivity (right) decrease in D. Giraldi et al. / Environmental Mo 2.6. Oxygen transfer The oxygen transport in liquid phase is described by means of the same transport equation used for the other dissolved compo- nents. However, in this case the reaction term accounts for both the biochemical reactions and the diffusive exchanges with the gaseous phase. The exchanges of oxygen between the liquid and the gaseous phases are modelled according to the Fick’s law. FITOVERT adopts a global coefficient of oxygen transport KLa [T �1], which can vary to account for the reduction of interchange surface due to the increase of water content: KLa ¼ � 1� q qs � KLa0 (7) where KLa0 [T �1] is referred to drained conditions. The oxygen transport in gaseous phase is described by means of a mass balance equation, with the following assumptions: pressure not specifically implemented, since prevalent literature considers the oxygen leaking from plant roots as secondary compared to other oxygen supply mechanisms, especially for VSSF-CWs (e.g. USEPA, 2000). Nevertheless, oxygen transfer by plants can be simulated by a specific calibration of the oxygen transfer parameters. 2.7. Evapotranspiration The potential evapotranspiration is assumed among the input data, since it can be easily estimated by well-established models (e.g. Priestley and Taylor, 1972; Allen et al., 1998). The evapotrans- piration is then considered as the sum of the evaporation from the free surface and the transpiration from plants in VSSF-CWs. The two terms are identified according to the leaf area index (Varado et al., 2006; Liu et al., 2005). The evaporation losses on the surface are covered by the water flux from the layers below, according to the suction head and the hydraulic conductivity of the system itself. gradient is neglected (i.e. atmospheric pressure all over the bed) (e.g. Massmann and Farrier, 1992); diffusion is the main transport mechanism for oxygen gas in steady water content conditions (Jury et al., 1991); the diffusion coefficient for oxygen in free airDa [L 2T�1] is corrected by means of a tortuosity factor (e.g. Patwardhan et al., 1988), according to the formulation of Marshall (1959): D ¼ ð3� qÞ32Da (8) where D [L2T�1] is the diffusion coefficient for oxygen in the bed and 3 [L3L�3] is the current porosity. Oxygen transfer by plants from the atmosphere to their roots is 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.5 1 Depth [cm] k [cm /s] day 0 time 00:00:00 day 4 time 06:51:51 day 8 time 13:43:31 day 12 time 20:35:11 day 17 time 03:26:51 day 21 time 10:18:31 day 25 time 17:10:11 day 30 time 00:00:00 op 5 cm of a VSSF-CW during a demonstrative 30 days simulation. 1.5 x 10-3 ng & Software 25 (2010) 633–640 637 The actual evaporation losses are equal to the potential ones when the medium is saturated, that is the water table is within the capillary fringe. The actual evaporation rate is defined according to Lappala et al. (1987). The transpiration losses depend on the root water uptake, and hence affect the deeper layers according to the root development. To describe the root water uptake, we adopted a microscopic approach (Varado et al., 2006), consisting in a sink term added to the Richard’s equation. The sink term depends on the water potential gradient between the medium and the plant, the hydraulic conductivity, the root density distribution and the atmospheric demand (Chabot et al., 2002). Many models have been proposed for the sink term, however the canonical form (Lai and Katul, 2000) was adopted. 2.8. Numerical implementation The mixed form of the Richards equation and the transport equation for the dissolved components were discretized in space delli using a Galerkin finite element (FE) method (e.g. Strang and Fix, 1973; Ern and Guermond, 2004). The domain was partitioned with a uniform grid into a set of finite elements connected at fixed nodal points. Next, through the variational form of the equation, the solution was approximated in the functional space of the contin- uous piecewise linear functions on the elements of the domain. The discretization of the mixed form of the Richards equation, compared to the h-based form, guarantees mass conservation. Regarding temporal discretization, the time derivative in the Richards equation was approximated with a finite difference scheme. Due to the non-linearity of the equation, the Picard itera- tive method (e.g. Dunham, 1921; Kangle et al., 1998) was adopted for the simplicity of the code and the computational efficiency. As result of the spatial and temporal discretization, a tridiagonal linear system with the nodal values of the pressure head as unknowns was obtained. This system was solved by means of the Thomas algorithm (e.g. Conte and de Boor, 1972). As the characteristic time steps for the solution of the biochemical system and the transport model are different, the chemical computations have been decoupled from the solution of the transport equation. Thus, at the beginning of every global time step, the non-linear first-order system of differential equations describing the biochemical process is numerically integrated, explicitly in time, with a finite difference scheme at each nodal point of the uniform grid. The transport equation was solved adopting the Galerkin finite element scheme for the discretization in space, and an implicit Cranck–Nicolson scheme for the discretization in time (e.g. Crank and Nicolson, 1947; Paniconi et al., 1991). The Cranck–Nicolson scheme is second-order accurate in time and unconditionally stable, i.e. it does not impose strict conditions on the time step. Nevertheless, without corrections, it does not guarantee that the solution will be greater than zero. The flux correction procedure (FEM-FCT) developed by Kuzmin and Turek (2002), which is suit- able for implicit time discretizations, was adopted to preserve the sign of the solution, and to provide a sharp resolution of steep gradients. The mass conservation equations for the particulate compo- nents are discretized in spacewith a Petrov–Galerkin finite element scheme, well suited for advection-dominated flows (Brooks and Hughes, 1982). The time derivative is discretized with a finite difference scheme, and the obtained non-linear discrete system is then solved with the same scheme used for the Richards equation. 3. Model evaluation 3.1. Validation of the hydraulic module A proper hydraulic calibration of mechanistic models describing CWs is considered the basis to obtain reliable overall simulations of VSSF-CW systems (e.g. Langergraber, 2002). Therefore we firstly evaluated FITOVERT by calibrating and testing its hydraulic module (Giraldi et al., 2009). Tests were carried out by comparing model outputs for inert components with experimental breakthrough curves of a pilot plant. The pilot plant was a 33 m2 VSSF-CW with six gravel layers for a total depth of 68 cm (Table 1), and plantedwith different plant species (Cyperus papyrus, Canna sp., Iris pseudoacorus, Phragmites australis and Juncus ensifolius). We analyzed three different satu- ration conditions under several constant flux regimes: complete saturation, partial saturation with the free water table 20 cm over the bottom of the bed, complete drainage (Table 2). The hydraulic parameters of the model were experimentally determined during the same measurement campaign on the pilot D. Giraldi et al. / Environmental Mo638 plant itself (Table 1), with the exception of the dispersivity coefficient, which was assessed using a trial and error procedure over experimental data from tracer tests. Therefore, the inputs of the model were in the number of 38 (five van Genuchten–Mualem parameters and thickness for each of the 6 layers of the bed; hydraulic load; inlet tracer concentration) with only one degree of freedom (the dispersivity). In the model, no sources or sinks of water and tracer were considered. This means that the mass recovery was theoretically 100%, and the tracer was perfectly conservative. Since this assumption was not matched in the tests, the simulations were performed assuming that the inlet tracer mass was equal to the mass of tracer actually leaving the system. This assumption implies that the effect of evapotranspiration is not considered significant for the hydrodynamics of CWs; this is not always true (Chazarenc et al., 2004), but in this case the hydraulic load was so high that the incidence of evapotranspiration can be neglected. The goodness of fit between simulated and measured data was evaluated using the Nash–Sutcliffe efficiency criteria (Nash and Sutcliffe, 1970). FITOVERT obtained satisfying tracer test simulations both for saturated and unsaturated conditions (Fig. 2 and Table 2). The adopted numerical scheme appeared to be very robust, and the numerical diffusion was avoided, as demonstrated by simulation tests not shown here. The hydraulic calibration of the model pointed out that the dispersivity parameter decreased as the saturation of the system increased (Table 2). This finding is in agreement with the literature for soil column experiments (Padilla, 1999) and was confirmed by a mechanical interpretation of the mixing phenomena (Giraldi et al., 2009). Since the water content typically varies dynamically during the standard operation of VSSF-CWs, the relation between dispersivity and saturation degree should be included in the numerical models, like FITOVERT, aimed to simulate the VSSF-CW behaviour. In our study, this relationship seemed to be linearly decreasing, but more data on different unsaturated conditions need to be analyzed to confirm this statement or to validate any other non-linear rela- tionships. Recently developed methods to estimate such relation- ship in porous media (Catania and Paladino, 2009) could also be used for VSSF-CW assessment. 3.2. Evaluation of the biochemical and transport modules Our tests demonstrated that a proper calibration of the hydraulic module is a mandatory prerequisite for the stable running of the complete model. Thus, the hydraulic module has been calibrated and validated in a number of situations, as shown in the previous section. Conversely, the calibration of the biochemical and transport modules is currently being performed with reference to experi- mental data, with the aim to validate the model for full-scale operation simulations and forecasts. Since calibrated simulations are not available yet, Fig. 3 presents some results of demonstrative simulations showing time variations of gaseous oxygen, dissolved oxygen, ammonia nitrogen and nitrite–nitrate nitrogen during a simulated filling–empting cycle of a CW plant. The calibration of the biochemical kinetics is particularly important, since literature values for parameters included in the implemented ASM1 pattern refer to suspended growth biomass, while VSSF-CWs typically support attached growth biomass. As already stated, the bio-film diffusion must be considered by way of a dedicated calibration of the biochemical module. Respirometry, which is an investigation technique of biological systems based on the dynamic measurement of the oxygen consumption, has been identified as the best tool to assess the biochemical dynamics of ng & Software 25 (2010) 633–640 VSSF-CWs (Giraldi and Iannelli, 2003), and therefore to calibrate particular mechanism is required too. This is being carried out on provide useful insights into the behaviour of this kind of waste- delli water treatment system. Acknowledgements The present work was partially funded by the Italian Ministry for the Environment, Land and Sea within the Italian-Israeli Cooperation on Environmental Technologies – Project 6. FITOVERT – Version 0.1 has been tested with the contribution of the following colleagues and friends: F. Castrogiovanni, A. Del Genovese, D. Fasano, G. Guido, E. Perrone, S. Tempestini. References Allen, R.G., Pereira, L.R., Raes, D., Smith, M., 1998. Crop Evapotranspiration – Guidelines for Computing Crop Water Requirements. FAO Irrigation and Drainage Paper, 56, Rome, Italy. Andreottola, G., Oliveira, E., Foladori, P., Peterlini, R., Ziglio, G., 2007. Respirometric techniques for assessment of biological kinetics in constructed wetlands. Water Sci. Technol. 56 (3), 255–261. Austin, D., 2006. Influence of cation exchange capacity (CEC) in a tidal flow, flood and drain wastewater treatment wetland. Ecol. Eng. 28 (1), 35–43. Babatunde, A.O., Zhao, Y.Q., O’Neill, M., O’Sullivan, B., 2008. Constructed wetlands for environmental pollution control: a review of developments, research and practice in Ireland. Environ. Int. 34, 116–126. Behrends, L., Houke, L., Bailey, E., Jansen, P., Brown, D., 2001. Reciprocating con- Crine et al., 1992; Se´guret et al., 2000), since these systems have several characteristics in common with VSSF-CWs. 4. Conclusions We have herein presented FITOVERT, a numerical model specifically developed to simulate and forecast the behaviour and treatment properties of VSSF-CWs. The first calibration and validation of the hydraulic modulewere satisfying for both saturated and unsaturated conditions, but highlighted the need for a better mathematical description of the relationship between dispersivity and saturation degree. Further experiments are required to fully calibrate the other modules that constitute themodel. Most of them have already been started. Once the calibration of FITOVERT is completed, the model could be used to optimize both design and operation of VSSF-CWs and to test new operation policies. The calibration process itself can the basis of experimental data of water content profiles over time measured on CW pilot plants byway of a capacitance probe (Giraldi and Iannelli, 2009). Since the calibration-validation process is still in progress, Fig. 4 shows a simulated porosity and hydraulic conductivity decrease during a demonstrative long-time run of the model. Finally, possible improvements of the model and procedures of calibration can be derived from the experience gained in trickling filter modelling (e.g. Kshirsagar et al., 1972; Sant’Anna et al., 1982; the model. Nevertheless this technique, whilst widely used to calibrate activated sludge models (e.g. Spanjers et al., 1998), still presents very limited examples of application to VSSF-CWs (Giraldo and Zarate, 2001; Andreottola et al., 2007), thus requiring further efforts for its extensive use as a standard calibration tool. Respiro- metric tests were carried out on the same VSSF-CW pilot plant used for the hydraulic calibration; the collected data are being used to calibrate also the other modules of FITOVERT. As FITOVERT is able to simulate the porosity and conductivity reduction due to biomass growth and particle transport and deposition throughout the filtering media, a calibration of this D. Giraldi et al. / Environmental Mo structed wetlands for treating industrial, municipal and agricultural waste- water. Water Sci. Technol. 44 (11–12), 399–405. Boller, M.A., Kavanaugh, M.C., 1995. Particle characteristics and headloss increase in granular media filtration. Water Res. 29 (4), 1139–1149. Bresler, E., 1973. Simultaneous transport of solutes and water under transient unsaturated flow conditions. Water Resour. Res. 9 (4), 975–986. Brooks, A.N., Hughes, T.J.R., 1982. Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompress- ible Navier–Stokes equations. Comput. Method Appl. M. 32 (1–3), 199–259. Brovelli, A., Baechler, S., Rossi, L., Langergraber, G., Barry, D.A., 2007. Coupled flow and hydro-geochemical modelling for design and optimization of horizontal flow constructed wetlands. In: Mander, U¨, Ko´iv, M., Vohla, C. (Eds.), 2nd Inter- national Symposium on ‘‘Wetland Pollutant Dynamics and Control WETPOL 2007’’ – Extended Abstracts, 16–20 September 2007, Tartu, Estonia, 2007, vol. II, pp. 393–395. Brovelli, A., Malaguerra, F., Barry, D.A., 2009. Bioclogging in porous media: model development and sensitivity to initial conditions. Environ. Model. Soft. 24 (5), 611–626. Catania, F., Paladino, O., 2009. Optimal sampling for the estimation of dispersion parameters in soil columns using an iterative genetic algorithm. Environ. Model. Soft. 24 (1), 115–123. Chabot, R., Bouarfa, S., Zimmer, D., Chaumont, C., Duprez, C., 2002. Sugarcane transpiration with shallow water table: sap flow measurements and modelling. Agric. Water. Manag. 54, 17–36. Chazarenc, F., Naylor, S., Brisson, J., Merlin, G., Comeau, Y., 2004. Effect of evapo- transpiration on hydrodynamics and performance of constructed wetlands. In: Proc. 9th International Conference on Wetlands system for Water Pollution Control, Avignon, 26th September–1st October 2004. Conte, S.D., de Boor, C., 1972. Elementary Numerical Analysis: An Algorithmic Approach, second ed. McGraw-Hill, New York, USA. Crank, J., Nicolson, P., 1947. A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type. Proc. Cambridge Philos. Soc. 43, 50–64. Crine, M., Marchot, P., L’Homme, G., 1992. Statistical hydrodynamics in trickle flow columns. Am. Inst. Chem. Engineers J. 38 (1), 136–147. Dunham, J., 1921. Note on the Picard method of successive approximations. Ann. Math. 23 (2), 75–77. Ern, A., Guermond, J.L., 2004. Theory and Practice of Finite Elements. Springer, New York, USA. Giraldi, D., de’Michieli Vitturi, M., Zaramella, M., Marion, A., Iannelli, R., 2009. Hydrodynamics of vertical subsurface flow constructed wetlands: tracer tests with rhodamine WT and numerical modelling. Ecol. Eng. 35 (2), 265–273. Giraldi, D., Iannelli, R., 2003. Tecniche respirometriche per lo studio dei processi biologici nei terreni e nelle biomasse adese (Respirometric techniques for the study of biological process in soils and attached biomass systems). In: Proc. 24� Corso di aggiornamento in Tecniche Per La Difesa Dall’inquinamento, Guardia Piemontese Terme (CS), June 18–21, 2003 (in Italian). Giraldi, D., Iannelli, R., 2009. Measurements of water content distribution in vertical subsurface flow constructed wetlands using a capacitance probe: benefits and limitations. Desalination 243, 182–194. Giraldo, E., Zarate, E., 2001. Development of a conceptual model for vertical flow wetland metabolism. Water Sci. Technol. 44 (11–12), 273–280. Green, M., Friedler, E., Ruskol, Y., Safrai, I., 1997. Investigation of alternative method for nitrification in constructed wetlands. Water Sci. Technol. 35 (5), 63–70. Henze, M., Gardy Jr., C.P.L., Gujer, W., Marais, G.v.R., Matsuo, T., 1987. Activated Sludge Model No. 1. Scientific and Technical Report No. 1. IAWPRC, London, UK. Iannelli, R., Giraldi, D., de’Michieli Vitturi, M., Perrone, E., Tempestini, S., 2008. Polishing Municipal Secondary Effluent for Stream Rehabilitation – Italian Israeli Cooperation on Environmental Technologies: Project 6. Final Report, January 2008. Ives, K.J., 1969. Theory of Filtration. In: Proceeding of the ‘‘International Water Supply Association, Eight Congress’’, Vienna, vol. I, pp. K1–K28. Iwasaki, I., 1937. Some notes on sand filtration. J. Am. Water Works Assoc. 29, 1591– 1602. Jury, W.A., Gardner, W.R., Gardner, W.H., 1991. Soil Physics. JohnWiley and Sons Inc., New York, USA. Kadlec, R.H., Knight, R.L., Vymazal, J., Brix, H., Cooper, P., Haberl, R., 2000. Con- structed Wetlands for Pollution Control. IWA Scientific and Technical Report 8. IWA Publishing, London, UK. Kadlec, R., Knight, R., 1996. Treatment Wetlands. CRC Press. Kangle, H., Mohanty, B.P., Leij, F.J., van Genuchten, M.T., 1998. Solution of the nonlinear transport equation using modified Picard iteration. Adv. Water Resour. 21 (3), 237–249. Kshirsagar, S.R., Phadke, N.S., Tipnis, S.S., 1972. Detention time studies in trickling filter. Indian J. Environ. Health 14 (1), 95–104. Kuzmin, D., Turek, S., 2002. Flux correction tools for finite elements. J. Comput. Phys. 175, 525–558. Lai, C., Katul, G., 2000. The dynamic role of root water uptake in coupling potential to actual transpiration. Adv. Water Resour. 23, 427–439. Langergraber, G., 2002. Calibration of a simulation tool for subsurface flow con- structed wetlands for wastewater treatment. In: Hassanizadeh, S.M., Schotting, R.J., Gray, W.G., Pinder, G.F. (Eds.), Developments in Water Science 47: Computational Methods in Water Resources, vol. 1. Elsevier Science BV, Amsterdam, The Netherlands, pp. 663–670. Langergraber, G., 2007. Simulation of the treatment performance of outdoor ng & Software 25 (2010) 633–640 639 subsurface flow constructed wetlands in temperate climates. Sci. Total Environ. 380, 210–219. Langergraber, G., 2008. Modeling of processes in subsurface flow constructed wetlands – a review. Vadoze Zone J. 7, 830–842. Langergraber, G., Giraldi, D., Mena, J., Meyer, D., Pen˜a, M., Toscano, A., Brovelli, A., Korkusuz, E.A., 2009. Recent developments in numerical modelling of subsur- face flow constructed wetlands. Sci. Total Environ 407 (13), 3931–3943. Langergraber, G., Simunek, J., 2005. Modeling variably saturated water flow and multi- component reactive transport in constructed wetlands. Vadose Zone J. 4, 924–938. Lappala, E.G., Healy, R.W., Weeks, E.P., 1987. Documentation of Computer Program VS2D to Solve the Equation of Fluid Flow in Variably Saturated Porous Media. Water-Resources Investigations Report 83-4099. US Geological Survey, Denver, Colorado, USA. Lee, B.H., Scholz, M., 2006. Application of the self-organizing map (SOM) to assess the heavy metal removal performance in experimental constructed wetlands. Water Res. 40 (18), 3367–3374. Liu, S., Graham, W.D., Jacobs, J.M., 2005. Daily potential evapotranspiration and diurnal climate forcing: influence on the numerical modelling of soil water dynamics and evapotranspiration. J. Hydrol. 309, 39–52. Marshall, T.J., 1959. The diffusion of gas trough porous media. J. Soil Sci. 10, 79–82. Marsili-Libelli, S., Checchi, N., 2005. Identification of dynamic models for horizontal subsurface constructed wetlands. Ecol. Mod. 187, 201–218. Massmann, J., Farrier, D.F., 1992. Effects of atmospheric pressure on gas transport in the vadose zone. Water Resour. Res. 28 (3), 777–791. Meyer, P.D., Rockhold, M.L., Gee, G.W., 1997. Uncertainty Analyses of Infiltration and Subsurface Flow and Transport for SDMP Sites. NUREG/CR-6565. Pacific Northwest National Laboratory Prepared for U.S. Nuclear Regulatory Commis- sion, Washington, DC, USA. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models, part I – a discussion of principles. J. Hydrol. 10, 282–290. O’Melia, C.R., Ali, W., 1978. The role of retained particles in deep bed filtration. Prog. Water Technol. 10, 167–182. Ogata, A., 1970. Theory of dispersion in a granular medium. U.S. Geol. Surv. Prof. Pap., 411-I, pp. 1–34. Padilla, I.Y., Jim Yeh, T.-C., Conklin, M.H., 1999. The effect of water content on solute transport in unsaturated porous media. Water Resour. Res. 35, 3303–3313. Paniconi, C., Aldama, A.A., Wood, E.F., 1991. Numerical evaluation of iterative and noniterative methods for the solution of the nonlinear Richards equation. Water Resour. Res. 27, 1147–1163. Patwardhan, A.S., Nieber, J.L., Moore, I.D., 1988. Oxygen, carbon dioxide and water transfer in soils: mechanism and crop response. Trans. ASAE 31 (5), 1383– 1395. Priestley, C.H.B., Taylor, R.J., 1972. On the assessment of surface heat flux and evaporation using large scale parameters. Mon. Weather Rev. 100 (2), 81–92. Raats, P.A.C., 2001. Developments in soil–water physics since the mid 1960s. Geo- derma 100 (3–4), 355–387. Rousseau, D.P.L., Vanrolleghem, P.A., De Pauw, N., 2004. Model-based design of horizontal subsurface flow constructed treatment wetlands: a review. Water Res. 38 (6), 1484–1493. Sant’Anna, G., Roustan, M., Roques, H., Nyadziehe, K.T., Ben Aim, R., 1982. Hydro- dynamics of plastic media trickling filters. Environ. Technol. Lett. 3, 395–404. Se´guret, F., Racault, Y., Sardin, M., 2000. Hydrodynamic behaviour of full scale trickling filters. Water Res. 34 (5), 1551–1558. Spanjers, H., Vanrolleghem, P., Olsson, G., Dold, P., 1998. Respirometry in Control of the Activated Sludge Process: Principles. Scientific and Technical Report No. 7. IAWQ, London, UK. Strang, G., Fix, G.J., 1973. An Analysis of the Finite Element Method. Prentice-Hall, Englewood Cliffs, USA. Sun, G., Gray, K.R., Biddlestone, A.J., Cooper, D.J., 1999. Treatment of agricultural wastewater in a combined tidal flow-downflow reed bed system. Water Sci. Technol. 40 (3), 139–146. Sun, G., Zhao, Y.Q., Allen, S.J., 2005. Enhanced removal of organic matter and ammoniacal-nitrogen in a column experiment of tidal flow constructed wetland system. J. Biotechnol. 115, 189–197. Sun, G., Zhao, Y.Q., Allen, S.J., Cooper, D., 2006. Generating ‘‘tide’’ in pilot-scale constructed wetlands to enhance agricultural wastewater treatment. Eng. Life Sci. 6 (6), 560–565. USEPA, 2000 Constructed Wetlands Treatment of Municipal Wastewaters – Manual, EPA/625/R-99/010. Cincinnati, USA. van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soil. Soil Sci. Soc. Am. J. 44, 892–898. van Genuchten, M.T., 1982. A comparison of numerical solutions of the one- dimensional unsaturated-saturated flow and mass transport equation. Adv. Water Resour. 5, 47–55. Varado, N., Braud, I., Ross, P.J., 2006. Development and assessment of an efficient vadose zone module solving the 1D Richards’ equation and including root extraction by plants. J. Hydrol. 323, 257–258. Zhao, Y.Q., Sun, G., Allen, S.J., 2004a. Purification capacity of a highly loaded labo- ratory scale tidal flow reed bed system with effluent recirculation. Sci. Total Environ. 330 (1–3), 1–8. Zhao, Y.Q., Sun, G., Allen, S.J., 2004b. Anti-sized reed bed system for animal wastewater treatment: a comparative study. Water Res. 38 (12), 2907–2917. D. Giraldi et al. / Environmental Modelling & Software 25 (2010) 633–640640 FITOVERT: A dynamic numerical model of subsurface vertical flow constructed wetlands Introduction Model description Topological configuration Hydraulic module Biochemical module Transport of dissolved components Transport of particulate components and clogging process Oxygen transfer Evapotranspiration Numerical implementation Model evaluation Validation of the hydraulic module Evaluation of the biochemical and transport modules Conclusions Acknowledgements References