VZJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 13 September 2005
Published in Vadose Zone J 4:924-938 (2005)
DOI: 10.2136/vzj2004.0166
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (8)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Langergraber, G.
Right arrow Articles by Simunek, J.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Langergraber, G.
Right arrow Articles by Simunek, J.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Langergraber, G.
Right arrow Articles by Simunek, J.
Related Collections
Right arrow Wetlands and Aquatic Processes
Right arrow Water Pollution
Right arrow Software
Right arrow Vadose Zone Processes and Chemical Transport

ORIGINAL RESEARCH

Modeling Variably Saturated Water Flow and Multicomponent Reactive Transport in Constructed Wetlands

Günter Langergrabera,* and Jirka Simunekb

a Institute of Sanitary Engineering and Water Pollution Control, BOKU Univ. of Natural Resources and Applied Life Sciences, Vienna, Muthgasse 18, A-1190 Vienna, Austria
b Dep. of Environmental Sciences, Univ. of California, Riverside, CA 92521, USA
* Corresponding author (guenter.langergraber{at}boku.ac.at)

Received 25 November 2004.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 THEORY
 APPLICATIONS TO PILOT-SCALE...
 DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
Constructed wetlands (CWs) are becoming increasingly popular worldwide for removing organic matter (OM), nutrients, trace elements, pathogens, or other pollutants from wastewater and/or runoff water. We present a multicomponent reactive transport model CW2D (i.e., Constructed Wetlands 2D), as an extension of the HYDRUS-2D variably saturated water flow and solute transport software package. CW2D was developed to model the biochemical transformation and degradation processes in subsurface-flow CWs. Such wetlands involve a complex mixture of water, substrate, plants, litter, and a variety of microorganisms to provide optimal conditions for improving water quality. The water flow regime in subsurface-flow CWs can be highly dynamic and requires the use of a transient variably saturated flow model. The biochemical components defined in CW2D include dissolved oxygen (DO), three fractions of OM (readily and slowly biodegradable, and inert), four N compounds (ammonium, nitrite, nitrate, and dinitrogen), inorganic P, and heterotrophic and autotrophic microorganisms. Organic N and organic P were modeled as part of the OM. The biochemical degradation and transformation processes were based on Monod-type rate expressions. All process rates and diffusion coefficients were assumed to be temperature dependent. Heterotrophic bacteria were assumed to be responsible for hydrolysis, mineralization of OM (aerobic growth), and denitrification (anoxic growth). Autotrophic bacteria were assumed to be responsible for nitrification, which was modeled as a two-step process. Lysis was considered to be the sum of all decay and sink processes. We demonstrate the performance of the model for one- and two-stage subsurface vertical flow CWs. Model simulations of water flow, tracer transport, and selected biochemical compounds are compared with experimental observations. Limitations of the model are discussed, and needs for model improvements are summarized.

Abbreviations: ASM, Activated Sludge Model • COD, chemical oxygen demand • CW, constructed wetland • CW2D, Constructed Wetlands 2D • DO, dissolved oxygen • IP, inorganic phosphorus • OM, organic matter • OM CI, inert OM • OM CR, readily biodegradable OM • OM CS, slowly biodegradable OM • PSCW, pilot-scale subsurface vertical flow constructed wetland • SSP, Small-Scale Plot • TOC, total organic carbon • XANb, autotrophic microorganisms, Nitrobacter • XANs, autotrophic microorganisms, Nitrosomonas • XH, heterotrophic microorganisms


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
THEORY
 APPLICATIONS TO PILOT-SCALE...
 DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
OUR ABILITY TO MODEL flow and transport processes in the vadose zone between the soil surface and the groundwater table has increased enormously during the past several decades (van Genuchten and Simunek, 2004). Vadose zone models are now routinely used for a large number of applications, both in research of subsurface processes and management of our soil and water resources. Typical applications are evaluation of the performance and effectiveness of engineered covers to minimize infiltration into underlying waste, quantifying groundwater recharge, evaluating the impact of climate and land-use changes on subsurface flow, identifying important factors in controlling subsurface flow and/or contaminant transport, and/or evaluating regional and global water cycles (Scanlon et al., 2002).

Many or most vadose zone models consider the transport of only one solute and assume that the fate and transport of this solute is independent of all other species that may be present in the soil solution. In reality, the soil solution is always a mixture of many chemical species and microorganisms that may mutually interact, create complexed species, precipitate or dissolve, affect each others degradation and decay, and/or compete with each other for sorption sites (van Genuchten and Simunek, 2004). Many important environmental problems hence require a simultaneous analysis of prevailing flow and transport processes with multiple chemical and biological reactions and transformations. Examples involving mostly geochemical reactions as listed by Simunek and Valocchi (2002) include acid mine drainage (Walter et al., 1994; Lichtner, 1996), radionuclide transport (Viswanathan et al., 1998), and reactive permeable barriers for aquifer remediation (Fryar and Schwartz, 1994). Examples that additionally involve biochemical processes include degradation of NTA (Tebes-Stevens et al., 1998), fate and transport of metal–organic mixed wastes (Rittmann and VanBriesen, 1996; VanBriesen, 1998), the formation of redox zones in organic-contaminated aquifers (Abrams et al., 1998; Essaid et al., 1995), C and N cycles in soils (Parton et al., 1988), bacteria-induced changes in the soil hydraulic properties (Rockhold et al., 2002), and complex biogeochemical processes in CWs (Langergraber, 2001, 2003).

Constructed wetlands are engineering structures used worldwide to improve water quality (Kadlec et al., 2000; Langergraber and Haberl, 2001; Haberl et al., 2003). They involve a complex mixture of water, substrate, plants, litter, and a variety of microorganisms to produce the optimal conditions for removing OM, N, P, and for decreasing the concentrations of toxic trace metals, organic chemicals, and pathogens. Constructed wetlands include surface flow and subsurface flow CWs (Kadlec and Knight, 1996). Surface-flow CWs are generally densely vegetated and typically have water depths of <0.4 m. In subsurface-flow CWs no free water level is visible. Subsurface vertical flow CWs with intermittent feeding are now used widely due to their efficiency in removing ammonia N (e.g., Langergraber and Haberl, 2001). Although much experience exists in constructing and operating such systems, their design is still based mostly on empirical rules, such as requiring a certain area per person or using simple first-order kinetic rates (e.g., Kadlec and Knight, 1996). However, Kadlec (2000) summarizes evidence that first-order models are inadequate for the design of treatment wetlands.

Understanding how CWs function is difficult due to the large number of physical, chemical, and biological processes that are active at the same time and mutually influence each other. Numerical simulators can represent valuable tools for analyzing and improving our understanding of complex systems such as CWs. Only few simulation results describing the behavior of subsurface flow CWs have been published. Schwager and Boller (1997) simulated tracer experiments and O2 transport in intermittent sand filters using HYDRUS (an older version of HYDRUS-1D) and MOFAT programs, respectively. However, they did not consider reactive transport. Other models focused mainly on a description of seasonal trends (e.g., Wynn and Liehr, 2001) or used simplifications such as first-order reaction rates to describe the degradation of a substance along the flow path (Werner and Kadlec, 2000). A detailed reaction model using a similar approach as described in this paper was used by Tanner et al. (1999) to model reactions in surface flow CWs.

To our knowledge, CW2D is the first model that describes detailed reactions in subsurface flow CWs in combination with a flow model capable of describing variably saturated flow. Our main objectives in this article are (i) to present the development of a multicomponent reactive transport model, CW2D (i.e., Constructed Wetlands 2D) (Langergraber, 2001), which is an extension of the variably-saturated water flow and solute transport program HYDRUS-2D and which preserves all numerical and graphical capabilities of HYDRUS-2D (Simunek et al., 1999); (ii) to demonstrate the performance of the model for one- and two-stage subsurface flow CWs; and (iii) to compare model simulations of water flow, tracer transport, and selected biochemical compounds against experimental observations.


   THEORY
TOP
 ABSTRACT
 INTRODUCTION
 THEORY
APPLICATIONS TO PILOT-SCALE...
 DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
Variably Saturated Flow
The water flow regime in subsurface-flow CWs can be highly dynamic because of intermittent loadings of wastewater, and hence requires the use of a transient variably saturated flow model. The computational module of the HYDRUS-2D software package numerically simulates two-dimensional variably saturated water flow using the Richards equation (Simunek et al., 1999). Although HYDRUS-2D offers various models for the soil hydraulic properties, in this study we used only the formulation proposed by van Genuchten (1980). This model has six parameters, namely the residual and saturated water contents, {theta}r and {theta}s, respectively; the saturated hydraulic conductivity, Ks, the van Genuchten shape parameters {alpha} and n, and Mualem's pore connectivity parameter l. The HYDRUS-2D model also offers a wide range of both system-dependent and independent boundary conditions. For common applications to CWs, two additional boundary conditions had to be implemented. One concerns surface ponding during wastewater loadings that exceed the infiltration capacity of a vertical flow CW, which typically occurs during its operation. Surface ponding was implemented as an extension of the atmospheric boundary condition. For applications to CWs to treat combined sewer overflow, it was also necessary to implement a limited outflow boundary condition at the bottom of the CW to enable modeling of the reduction of peak flows. Reduction of peak flows is used to increase the treatment performance and to reduce the impact of combined sewer overflow on the receiving water (Dittmer et al., 2004).

Solute Transport
The governing equation for the macroscopic transport of component i can be written in the form

[1]
where i = 1,..., N (N is the number of components), ci is the concentration in the aqueous phase (M L–3), si is the concentration in the solid phase (M M–1), {theta} is the volumetric water content (L3 L–3), {rho} is the soil bulk density (M L–3), Di is the effective dispersion tensor (L2 T–1), q is the volumetric flux density (L3 L–2 T–1), S is the source–sink term (L3 L–3 T–1); cs,i is the concentration of the source–sink (M L–3), and ri is the reaction term (M L–3 T–1). Components of the effective dispersion tensor Di include molecular diffusion, and longitudinal and transverse dispersion. The reaction term ri in the multicomponent reactive transport model CW2D is described in detail below. Liquid- and solid-phase concentrations at equilibrium can be related using either linear, Frendlich, or Langmuir adsorption isotherms (Simunek et al., 1999).

The concept of two-region, dual-porosity type solute transport (van Genuchten and Wierenga, 1976) is implemented in HYDRUS-2D to permit consideration of physical nonequilibrium transport. The physical nonequilibrium transport model divides the liquid phase into mobile (flowing) and immobile (stagnant) regions. Solute exchange between the mobile and immobile regions is modeled as a first-order process.

The CW2D Multicomponent Reactive Transport Module
The need for a numerical model describing biochemical transformation and degradation processes in subsurface flow CWs motivated the development of the CW2D multicomponent reactive transport module (Langergraber, 2001) for incorporation into the HYDRUS-2D variably saturated water flow and solute transport program. CW2D is able to model the biochemical transformation and degradation processes for OM, N, and P. The reactive transport module CW2D considers 12 components and nine processes. The components include DO, OM (three fractions of different degradability), ammonium, nitrite, nitrate, and N2 gas, inorganic P, and heterotrophic and autotrophic microorganisms. Organic N and organic P were modeled as nutrient contents of the OM. The processes considered are hydrolysis, mineralization of OM, nitrification (modeled as a two-step process), denitrification, and a lysis process for the microorganisms.

The mathematical structure of CW2D was based on the structure of the Activated Sludge Models (ASMs) introduced by Henze et al. (2000). Major assumptions in all ASMs are (a) a constant value of pH, (b) constant coefficients in the rate equations, and (c) constant stoichiometric factors. The various ASMs are valid only for domestic wastewater and are applicable within the temperature range of 10 to 25°C (Henze et al., 2000). ASMs assume that nitrification occurs in one step, i.e., that ammonia is directly converted (nitrified) into nitrate. The two-step nitrification model used in CW2D (i.e., sequential nitrification of ammonium into nitrite and nitrate) is based on models presented by Nowak (1996) and Brouwer et al. (1998).

Components
The components defined in CW2D (Table 1) are:


View this table:
[in this window]
[in a new window]
  		Table 1. Components defined in the reactive transport module.

 
Organic N and organic P were modeled as N and P contents of the COD using the composition parameters as given below. We assumed that the OM was present only in the aqueous phase and that all reactions occurred only in the aqueous phase as well. Adsorption was considered for ammonium, NH4+, and IP. Adsorption was assumed to be a kinetic process.

Biofilm in CW2D was assumed to be ideal, composed of a homogeneous matrix of microorganisms, with all species uniformly distributed over the biofilm. All microorganisms were assumed to be immobile and hence were associated exclusively with the solid phase. Lysis was considered to be the sum of all decay and sink processes of all microorganisms involved, with the rate of lysis being independent of environmental conditions.

Reactions
The biochemical degradation and transformation processes in CW2D were modeled using Monod-type rate expressions as used in the ASMs (Henze et al., 2000). Equations [4] through [12] give the reaction rates rcj for the nine processes considered. Parameters in these equations are defined in Table 2. The various rate parameters in the kinetic growth and decay expressions were assumed to be temperature dependent. Values in Table 2 are given for a reference temperature of 20°C and for temperature dependent parameters also for 10°C. The temperature dependence of all process rates and diffusion coefficients was assumed to be described using the Arrhenius equation.


View this table:
[in this window]
[in a new window]
  		Table 2. Kinetic parameters (at 20°C, values in brackets at 10°C).{dagger}

 
Two factors are used in Eq. [4] through [12] as an abbreviation. The first factor describes the conversion of solid-phase concentrations into liquid-phase concentrations:

[2]
where the subscripts Y = H, ANs, and ANb are used for the heterotrophic and two autotrophic organisms, respectively. The second factor is the Monod term for nutrients:

[3]
where the subscript x for Het, DN, and ANb represents the processes of heterotrophic growth, denitrification and growth of Nitrobacter, respectively.

In this study we assumed heterotrophic bacteria to be responsible for hydrolysis, mineralization of OM (aerobic growth), and denitrification (anoxic growth).

Hydrolysis.
Hydrolysis describes the conversion of slowly biodegradable OM CS into readily biodegradable OM CR, with a small fraction being converted into inert OM CI. Ammonium (NH4+) and IP are released during this transformation process. We further assume that hydrolysis takes place independently of the O2 conditions, leading to the following rate equation for hydrolysis:

[4]

Aerobic Growth of Heterotrophic Bacteria.
This process consumes O2 and OM CR, while NH4+ and IP are incorporated in the biomass as follows:

[5]

Nitrate-Based Growth of Heterotrophs on Readily Biodegradable COD (denitrification).
This process consumes NO3– and OM CR. Nitrate is reduced to dinitrogen (N2). Again, NH4+ and IP are incorporated in the biomass as follows:

[6]
Inhibition due to NO3––N is modeled as a noncompetitive process using a soluble nonbiodegradable inhibitor (Nowak et al., 1995).

Nitrite-Based Growth of Heterotrophs on Readily Biodegradable COD (denitrification).
Similarly as for NO3–, this process consumes NO2–, OM CR, NH4+, and IP. Nitrite is reduced to N2. The inclusion of both denitrification processes is necessary because of the use of a two-step nitrification model.


[7]

Lysis of Heterotrophic Organisms.
Lysis produces OM (CS, CR, and CI), NH4+, and IP. Lysis is assumed to represent the sum of all decay and sink processes.


[8]
Autotrophic bacteria are responsible for nitrification, which is modeled as a two-step process. Nitrosomonas spp. consume NH4+ and produce NO2–, while Nitrobacter spp. transform NO2– to NO3–. Both steps of the nitrification process are strictly aerobic and can therefore only occur when DO is available.

Aerobic Growth of Nitrosomonas on Ammonium (First Step of Nitrification).
This process consumes NH4+ and O2 and produces NO2–. Inorganic P and a small portion of NH4+ are incorporated in the biomass.


[9]

Lysis of Nitrosomonas.
This lysis produces OM (CS, CR, and CI), NH4+, and IP as follows:

[10]

Aerobic growth of Nitrobacter on Nitrite (Second Step of Nitrification).
This process consumes NO2– and O2 and produces NO3–. Ammonium and IP are incorporated in the biomass as follows:

[11]

Lysis of Nitrobacter.
This process produces OM (CS, CR, and CI), NH4+, and IP as follows:

[12]

Stoichiometric Matrix
Table 3 represents the stoichiometric matrix for CW2D. The concept of yield coefficients is used to describe biomass growth. Yield coefficients represent the amount of substrate that is used for biomass growth. Stoichiometric coefficients (including the yield coefficients) that appear in Table 3 are defined, together with their default values, in Table 4 (for OM and microorganisms), Table 5 (for NH4+), and Table 6 (for IP).


View this table:
[in this window]
[in a new window]
  		Table 3. Stoichiometric matrix of reactions in CW2D.{dagger}

 

View this table:
[in this window]
[in a new window]
  		Table 4. Stoichiometric parameters for organic matter and microorganisms and composition parameters.

 

View this table:
[in this window]
[in a new window]
  		Table 5. Stoichiometric coefficients for NH+4–N.{dagger}

 

View this table:
[in this window]
[in a new window]
  		Table 6. Stoichiometric coefficients for inorganic phosphorus.{dagger}

 
Table 4 defines the stoichiometric and composition parameters used in Table 3. The composition parameters provide information about the N and P contents of the OM and biomass, and therefore are measures of the mass of organic N and organic P, respectively.

Stoichiometric parameters for NH4+–N, {upsilon}1,N through {upsilon}9,N, and IP, {upsilon}1,P through {upsilon}9,P, in Table 3 can be calculated from the mass balance equations. These mass balance equations can be written for each process j as follows:

[13]
where i = 1,..., N (N is the number of components), j = 1,..., R (R is the number of processes), {upsilon}j,i is the stoichiometric factor for component i and process j (mgi mgj–1), c represents the material to which the mass balance is applied (either N or P), and ic,i is a conversion factor that converts the unit of component i to the unit of material c. The resulting equations for the stoichiometric coefficients {upsilon}j,N, and {upsilon}j,P are given in Tables 5 and 6, respectively.

Calculation of the Reaction Term
It is assumed that reactions take place only in the aqueous phase. Therefore the overall reaction term ri for component i can be calculated as follows:

[14]
where as before i = 1,..., N, j = 1,..., R, and where ri is the reaction term for component i (M L–3 T–1), {theta} is the volumetric water content (L3 L–3), {upsilon}j,i is the stoichiometric factor for component i and process j (M M–1) (Table 3), and rcj is the zero-order reaction rate for process j in the aqueous phase (M L–3 T–1) as defined in Eq. [4] through [12].

Exchange of Oxygen
The exchange of O2 from the gas phase into the aqueous phase is described with a model commonly used in wastewater treatment (e.g., Gujer and Boller, 1990):

[15]
where a is the air content (L3 L–3), kaer,O2 is the O2 reaeration rate (T–1), and cO2,sat is the temperature dependent saturation concentration of DO (M L–3). This term must be added to the reaction rate constant for O2 since this is the only way to increase the concentration of DO (with the exception of water flow). Note that the exchange of O2 in Eq. [15] is proportional to the level of disequilibrium in the liquid phase and that the effect of the gas phase is included only by multiplying the rate with the air content.

Implementation of CW2D into HYDRUS-2D
The original HYDRUS-2D code considers consecutive first-order decay and degradation, involving the sequential degradation of one into another solute. This approach allows one to solve the transport equations for individual species of the decay chain sequentially since the second solute depends on the first one (and so on), but not vice versa. However, since various components in the new reactive CW2D module depend in a relatively complex manner on each other (e.g., degradation of one species produces multiple other species), the solute transport equations for all components must be evaluated simultaneously at each iteration step of the numerical solution. The iteration procedure was modified such that the iteration included all components simultaneously, rather than being performed at the level of each single component. If convergence was not reached for any single component, all components were recalculated during the next iteration. A self-adjusting time stepping scheme was required to provide a relatively smooth and efficient numerical solution. We note that Galerkin finite element method was still used to numerically solve both the variably saturated flow and solute transport equations.


   APPLICATIONS TO PILOT-SCALE CONSTRUCTED WETLANDS
TOP
 ABSTRACT
 INTRODUCTION
 THEORY
 APPLICATIONS TO PILOT-SCALE...
DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES
 
Pilot-Scale Subsurface Vertical Flow Constructed Wetland
System Description
Figure 1 shows a schematic of a pilot-scale subsurface vertical flow constructed wetland (PSCW). The PSCW is based on Austrian standards for vertical flow CWs for wastewater treatment (ÖNORM B 2505, 1997). The surface area of the wetland was 1 m2, and the height of the main layer of the filter bed 60 cm. The main layer consisted of sand (gravel size 0.06–4 mm or 1–4 mm in two PSCWs, respectively). The sandy material used for the main layer had a porosity of 0.30 and a saturated hydraulic conductivity, Ks, of 117 cm h–1. An intermediate layer of 10-cm thickness with a gravel size of 4 to 8 mm prevents fine particles from being washed out into the drainage layer (15 cm thick; gravel 16–32 cm) where the effluent was collected by means of tile drains (Langergraber, 2003).



View larger version (13K):
[in this window]
[in a new window]
  		Fig. 1. Schematic of the pilot-scale subsurface constructed wetlands (values are in centimeters).

 
The PSCW was loaded intermittently with wastewater. Following Austrian standards, the daily loading rate was 40 L at application intervals of 6 h (i.e., a single loading comprised 10 L).

Tracer experiments using NaCl were run for daily hydraulic loading rates of 40 and 60 L. The electrical conductivity of the effluent, which is directly related to the tracer concentration, was measured automatically. The measured influent conductivity of the tracer was 30 mS cm–1, while the background conductivity during the tracer experiment was 1.5 mS cm–1. Influent and effluent concentrations of various N species and total organic carbon (TOC) were measured periodically. Table 7 shows the median values and confidence intervals of the measured influent and effluent concentrations for NH4–N, NO3–N, and TOC.


View this table:
[in this window]
[in a new window]
  		Table 7. Median values and confidence intervals (in parenthesis) of the measured influent and effluent concentrations of NH+4, NO–3, and total organic carbon (TOC) (Langergraber, 2003).

 
The width of the transport domain in the numerical simulations was 1 m and its depth 0.6 m (only the main layer was simulated), while the transport domain itself was discretized into 11 columns and 40 rows. This results in a two-dimensional finite element mesh consisting of 440 nodes and 780 triangular finite elements. An atmospheric boundary condition was assigned to the top of the system, a constant pressure head boundary condition (constant head of –2 cm) to the bottom of the 0.06- to 4-mm main layer and a free drainage boundary condition to the bottom of the 1- to 4-mm main layer. Water flow and single-solute transport were evaluated using the standard HYDRUS-2D software, while multicomponent reactive transport was calculated using the CW2D/HYDRUS-2D software.

Flow Simulations
Table 8 shows various parameter sets of the van Genuchten–Mualem soil hydraulic property model used in flow simulations (Fig. 2 and 3) . The literature Parameter Set 1 was based on default parameters provided by HYDRUS-2D for sand using pedotransfer functions published by Carsel and Parrish (1988). The shape factors for sand were also used for Parameter Set 2 in Table 8, but now with the measured values of the saturated water content {theta}s and the saturated hydraulic conductivity Ks. Parameter Set 3 in Table 8 was obtained by calibrating the model against cumulative effluent data, while Parameter Set 4 represents parameters optimized using the effluent rates. Parameter Set 5 are those reported by Langergraber (2003), who in addition to the effluent rates, also used water content measurements at the 25-cm depth, while simultaneously optimizing the saturated hydraulic conductivity Ks. Finally the parameters used for the coarser gravel (1–4 mm) are also given.


View this table:
[in this window]
[in a new window]
  		Table 8. Literature and fitted soil hydraulic parameters of the van Genuchten–Mualem model for the 0.06- to 4-mm substrate.

 


View larger version (18K):
[in this window]
[in a new window]
  		Fig. 2. Measured and simulated cumulated effluent for a single application of 10 L (four loadings per day; daily hydraulic load of 40 L).

 


View larger version (21K):
[in this window]
[in a new window]
  		Fig. 3. Measured and simulated effluent rates for a single application of 10 L (four loadings per day; daily hydraulic load of 40 L).

 
Figures 2 and 3 show the measured and simulated cumulated effluents and effluent flow rates, respectively. Simulation results are presented for the parameter sets given in Table 8. The simulation using Parameter Set 3 showed a good match with the measured cumulated effluent (Fig. 2). However, peak effluent rates were significantly overpredicted (Fig. 3). Parameter Set 4 described the effluent rates much better. These results indicate that a good calibration of the flow model requires measurements of at least the saturated water content {theta}s and the saturated hydraulic conductivity Ks of the main layer, as well as the effluent flow rate between two successive loadings. However, the best fit was obtained (Fig. 4) with the parameter set of Langergraber (2003); this indicates that additional water content data and more degrees of freedom in the optimization (i.e., one additional parameter was fitted) can result in a better calibration of the flow model.



View larger version (13K):
[in this window]
[in a new window]
  		Fig. 4. Measured and simulated effluent rates (left) and cumulated effluent (right) for a single feeding of 10 L (four loadings per day, daily hydraulic load of 40 L).

 
Simulation of Single-Component Transport
Results for a single component solute transport were previously presented by Langergraber (2003) and are given here for completeness. Table 9 shows parameters that were used for the general equilibrium solute transport model, and Table 10 gives parameters for the physical nonequilibrium transport model. The equilibrium transport model assumes that all pore water is mobile, whereas the physical nonequilibrium transport model considers a fraction of the pore water to be immobile (stagnant). Figure 5 compares simulated and measured breakthrough curves. It can be seen that for both hydraulic loading rates the simulation results matched the measured data relatively well.


View this table:
[in this window]
[in a new window]
  		Table 9. Parameters of the general equilibrium solute transport model (Langergraber, 2003).

 

View this table:
[in this window]
[in a new window]
  		Table 10. Parameters for the physical nonequilibrium transport model (Langergraber, 2003).

 


View larger version (23K):
[in this window]
[in a new window]
  		Fig. 5. Measured and simulated effluent electrical conductivities for different daily hydraulic loads.

 
Simulation of Multicomponent Reactive Transport
Table 11 shows measured influent and effluent concentrations of selected components, as well as the influent concentrations used in the CW2D/HYDRUS-2D simulation, and the resulting simulated effluent concentrations. The simulated and measured effluent concentrations of NH4+ are lower for the 0.06- to 4-mm gravel than for the coarser material (gravel size of 1–4 mm). The higher flow rates in the main layer with the coarse material reduced the contact time between water and the bacteria. The effluent concentrations of NH4+–N and the OM CS were significantly higher for the coarser material. Effluent NO3– concentrations were also higher because of less readily degradable OM CR produced by hydrolysis that could be utilized for denitrification in deeper layers.


View this table:
[in this window]
[in a new window]
  		Table 11. Measured and simulated influent and effluent concentrations.

 
Results of the reactive transport simulations showed an overall good match with the measured data. Measured higher NH4+–N and NO3––N concentrations in the system with the 1- to 4-mm substrate could be well described by the model. However, further investigations regarding parameters of the biokinetic model are still required. The close match could be partly explained by the calibration of the flow model, which was made possible by having data that provided a good definition of the flow processes in the system.

Below we present several figures showing time series of the simulation results for two loadings at three locations (5-, 15-, and 50-cm depths) in the 60-cm main layer of the PSCW system filled with gravel size 0.06 to 4 mm. Figure 6 shows time series for the pressure head and DO. Oxygen depletion at depths of 5 and 15 cm can be observed to occur immediately after loading due to degradation. However, the O2 concentration at 50 cm stayed constant at a level close to the saturation concentration during the entire loading interval. After <2 h, the O2 concentration at depths of 5 and 15 cm reached saturation again.



View larger version (15K):
[in this window]
[in a new window]
  		Fig. 6. Pressure head h (left) and dissolved oxygen (DO) (right).

 
Time series and vertical concentration profiles at the end of the simulation (12 h) are shown in Fig. 7 for heterotrophic organisms XH. During two loading intervals the concentrations of the microorganisms remained constant, thus indicating quasi steady state conditions when growth and lysis processes are in balance (Fig. 7, left). Almost all microorganisms were present in the top 20 cm of the filter (Fig. 7, right), where most transformations and degradation processes took place. Autotrophic bacteria were distributed in a very similar manner. The maximum concentration of autotrophic bacteria was simulated to be lower than the maximum concentration of heterotrophic bacteria (data not shown).



View larger version (11K):
[in this window]
[in a new window]
  		Fig. 7. Simulated time series (left) and vertical profile (right) of heterotrophic organisms XH.

 
Figure 8 shows time series for NH4+–N and NO3––N, respectively. Effluent concentrations for NH4+–N and NO3––N were about 0.01 and 41 mg L–1, respectively. Vertical concentration profiles at simulation times of 1 min (end of loading), 6 min, 30 min, and 2 h (only for NO3–) are shown in Fig. 9 for NH4+–N and NO3––N. Notice that the concentration peaks for NH4+ at 5 and 15 cm coincide with minima in the NO3– simulations. By contrast, concentrations of NH4+ and NO3– remained very stable at 50 cm depth.



View larger version (16K):
[in this window]
[in a new window]
  		Fig. 8. Simulated time series for NH4+–N (left) and NO3––N (right).

 


View larger version (13K):
[in this window]
[in a new window]
  		Fig. 9. Vertical profiles of NH4+–N (left) and NO3––N (right) at t = 1 min (end of loading), 6 min, 30 min, and 2 h.

 
Two-Stage Subsurface Vertical Flow Constructed Wetland
System Description
The Small-Scale Plot (SSP) system described by Grosse et al. (1999) is a PSCW for the treatment of heavily polluted surface water. The total surface area of the SSP was 2 m2, divided into downflow and upflow chambers (each having a surface area of 1 m2). The inlet was situated on top of the downflow chamber. The effluent was collected in the upflow chamber by means of perforated pipes. A separating wall was built between the downflow and upflow chambers. Figure 10 shows a schematic of the SSP (Perfler et al., 1999). Heights of the downflow and upflow chambers were 85 and 60 cm, respectively. The drainage layer at the bottom of both chambers had a height of 15 cm. This layer was connected to the two chambers and contained two tile drains. The main 55- and 45-cm-thick layers in the downflow and upflow chambers, respectively, were filled with sand (gravel size 0.06–4 mm). The top layer of the downflow chamber had a thickness of 15 cm (gravel size 4–8 mm). Algae growth in the main layer was prevented by keeping the water table during application beneath the soil surface.



View larger version (18K):
[in this window]
[in a new window]
  		Fig. 10. Schematic of the two-stage vertical flow constructed wetlands (distances are in centimeters).

 
The unstructured two-dimensional finite element mesh used for the SSP simulations consisted of 1135 nodes and 2057 elements. The inlet distribution pipe at the top of the downflow chamber was modeled as a time dependent atmospheric boundary condition. Outlet (effluent) drainage pipes were modeled as circular openings in the transport domain, and subjected to seepage face boundary conditions. The soil hydraulic parameters of the van Genuchten–Mualem model used for the simulations are presented in Table 12. Values for the sand were the same as those in Table 8 (Parameter Sets 1 and 2), while the parameters for the top and drainage layer were again the same as those obtained by Langergraber (2001).


View this table:
[in this window]
[in a new window]
  		Table 12. van Genuchten–Mualem soil hydraulic parameters for the two-stage constructed wetlands system.

 
Since simulations of single-component reactive transport using the physical equilibrium and nonequilibrium (mobile–immobile water) transport models were presented by Langergraber (2003), we will focus below on multicomponent reactive transport. Langergraber (2003) showed that the presence of stagnant (immobile) regions of the pore water produced a much better match with the measured concentration data, compared with considering only mobile pore water.

Simulation of Multicomponent Reactive Transport
Table 13 shows influent concentrations for the multicomponent reactive transport simulations. Since this system was used to treat polluted surface water, the concentrations were much lower compared with wastewater. The hydraulic loading rate was 50 L every 2 h (i.e., 600 L d–1).


View this table:
[in this window]
[in a new window]
  		Table 13. Influent concentrations for the multicomponent reactive transport simulations.

 
Figure 11 shows the distribution of heterotrophic organisms XH. The maximum XH concentration was found to be in the upper part of the main layer of the downflow chamber. Bacteria concentrations decreased with depth. Concentrations were low in the upflow chamber because of a lack of substrate available.



View larger version (13K):
[in this window]
[in a new window]
  		Fig. 11. Simulated steady state distribution of heterotrophic organisms XH.

 
Figure 12 shows the distribution of Nitrosomonas XANs. The distribution of the second autotrophic organism, Nitrobacter XANb, was found to be very similar (not further shown here). Again, most autotrophic bacteria were located in the upper part of the downflow chamber. The maximum concentration of autotrophic bacteria was lower than the maximum concentration of the heterotrophic bacteria. The concentration decreased faster with depth for the autotrophic bacteria than for the heterotrophic bacteria. This is a direct consequence of nitrification being a strictly aerobic process, with O2 being available only in the upper parts of the downflow chamber (see also Fig. 13) .



View larger version (12K):
[in this window]
[in a new window]
  		Fig. 12. Calculated steady state distribution of Nitrosomonas XANs.

 


View larger version (22K):
[in this window]
[in a new window]
  		Fig. 13. Simulated time series of dissolved oxygen (DO).

 
Figure 13 shows time series of DO for two consecutive loadings to the main layer in the downflow chamber (0 cm) at two different depths (5 and 15 cm, located in the unsaturated and saturated zones, respectively). After loading, the DO concentration decreased at the surface and at the 5-cm depth. At the 15-cm depth, however, the DO concentration increased due to advective DO transport with the infiltrating water. Still, O2 decreased quickly (within 15 min after loading) due to the consumption of O2. No O2 was found at the 15-cm depth (i.e., in the saturated zone) during the remainder of the simulation. In the unsaturated zone (5-cm depth) the DO concentration increased again due to re-aeration, and reaching O2 saturation after about 1 h after loading.

Figure 14 shows time series of NH4+–N and NO3––N at the surface of the main layer (0 cm) and at three different depths (5 cm, located in the unsaturated zone, and 15 and 30 cm, located in the saturated zone). After loading, NH4+ initially increased whereas NO3– decreased. Later, NH4+ is nitrified, resulting in higher NO3– concentrations. At the 15-cm depth, changes in the concentration occurred mainly due to advection since O2 was already a limiting factor (Fig. 13). Notice, again, that only small changes in the concentrations occurred at the 30-cm depth. Figure 15 shows simulated and measured concentrations of NH4+–N and NO3––N in a cross section of the downflow chamber. The simulation results matched the measured data well.



View larger version (18K):
[in this window]
[in a new window]
  		Fig. 14. Simulated time series of NH4+–N (left) and NO3––N (right).

 


View larger version (18K):
[in this window]
[in a new window]
  		Fig. 15. Simulated and measured concentrations of NH4+–N and NO3––N in a cross section of the downflow chamber.

 

   DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 THEORY
 APPLICATIONS TO PILOT-SCALE...
 DISCUSSION
SUMMARY AND CONCLUSIONS
 REFERENCES
 
Our simulation results showed a good match with the measured data for both single-component (tracer experiments, e.g., Fig. 5) and multicomponent reactive transport (e.g., Table 11 and Fig. 14) when the hydraulic behavior of the system was accurately described (e.g., effluent flow rates in Fig. 2 and 4). The importance of the water flow model was also evident from a sensitivity analysis that showed that results of the reactive transport simulations were very sensitive to model parameters describing the unsaturated properties of the substrate. The good match of the experimental data with the reactive transport simulations could be obtained using literature values for the CW2D model parameters (Langergraber, 2001).

The CW2D module considers the main constituents of domestic wastewater. CW2D hence is only applicable to CWs that treat domestic wastewater or water with properties that can be described with processes and parameters currently incorporated into CW2D. If other components and/or waters are to be modeled, an extension to the present model is needed.

Further limitations of CW2D result from assumptions of the basic model, the Activated Sludge Model, used here. These include consideration of (i) a constant value of pH, (ii) constant coefficients in the rate equations, (iii) constant stoichiometric factors, and (iv) a limited temperature range from 10 to 25°C (Henze et al., 2000).

Although simulations showed a good match with the measured data, much research remains so as to produce a reliable tool for CW design. The main research needs regarding modeling of CW processes as listed by Langergraber (2003) are:

Detailed Hydraulic Investigations of Full-Scale CWs.
A good match between simulation results and measured data can be obtained when sufficient data are available to describe the hydraulic behavior of the system, and thus when the flow model can be successfully calibrated. For full-scale systems, the match between the model and data was thus far relatively poor, most probably due to a lack of data to accurately describe the hydraulic behavior of the system (Langergraber, 2001) and due to the inherent heterogeneity of large-scale systems. The heterogeneity of the substrate and uneven distribution of influent water on the surface are the main reasons for needing detailed hydraulic investigations of full-scale systems.

Numerical Simulations of Outdoor CW Systems.
Simulations of outdoor (full-scale) CW systems are also needed to validate the temperature dependencies of the various processes considered in CW2D.

Investigation of Plant Uptake Models.
Models describing nutrient (N and P) uptake by plants have to be considered to describe the contribution of plants to the overall treatment process. Langergraber (2004) showed the important influence of plant uptake on effluent concentrations and removal efficiency using plant uptake models that are currently available in HYDRUS-2D (Simunek et al., 1999). These models describe nutrient uptake associated with water uptake (i.e., active uptake). For high loaded systems (e.g., CWs that treat wastewater) the model gives good results, whereas for waters with lower nutrient contents the simulation results indicate that potential nutrient uptake will be substantially overestimated when using literature values to calculate the potential uptake rate, compared with simulations using models that couple nutrient uptake to water uptake.

Need for a Module Describing Substrate Clogging.
One limitation of CW2D is that up to now only dissolved solutes are being considered. Wastewater also contains suspended solids that may cause pore clogging. Clogging is considered to be the biggest operational problem affecting vertical flow in CWs. It causes a reduction in the infiltration capacity of the substrate surface, and therefore often a rapid decline in the treatment performance of the system (e.g., Platzer and Mauch, 1997; Winter and Goetz, 2003; Langergraber et al., 2003). The clogging module should be able to describe pore size reductions due to the settling of suspended solids and due to bacteria growth. Both particle transport and bacteria growth models are currently available (e.g., Mackie and Bai, 1993; Vandevivere et al., 1995; Kildsgaard and Engesgaard, 2001) and could be incorporated into CW2D/HYDRUS-2D. The effects of pore size reductions on the hydraulic properties of the substrate may also need to be taken into account (e.g., Taylor and Jaffé, 1990; Magnico, 2000).

Improved Experimental Methods for Estimating CW2D Model Parameters.
There is also a need for improving technologies to measure CW2D model parameters rather than always requiring calibration of the reactive transport part of the model. This is important before the model can be used as a design tool. A number of measurement techniques have been developed to characterize the parameters of Activated Sludge Models as applied to activated sludge wastewater treatment plants (Henze et al., 2000; Vanrolleghem et al., 1999). However, no experimental techniques are currently available to measure CW2D model parameters. Detailed investigations on microbial biocoenosis in particular must be performed, with a need to adopt existing measurement methods to soil and/or wastewater environments.


   SUMMARY AND CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 THEORY
 APPLICATIONS TO PILOT-SCALE...
 DISCUSSION
 SUMMARY AND CONCLUSIONS
REFERENCES
 
We presented a multicomponent reactive transport model, CW2D, which is an extension of the HYDRUS-2D variably saturated water flow and solute transport program, to include biochemical transformation and degradation processes in subsurface flow CWs. Constructed wetlands represent a very complex mixture of water, substrate, plants, litter, and a variety of microorganisms, in which a large number of simultaneous processes takes place. A complete description of such a complex system involving mutually dependent physical, chemical, and biological reactions and chemical and biological compounds is nearly impossible without the help of a numerical model that considers the majority of these interactions. The water flow regime in subsurface vertical flow CWs is also highly dynamic, adding to the complexity of the overall system by requiring a transient variably saturated flow model. By developing CW2D/HYDRUS-2D we have attempted to develop a numerical tool that can capture at least the main processes of such a complex system, and that can help in analyzing and studying the various transient processes operative in CWs.

In this paper we described the main processes and components considered in CW2D and demonstrated the performance of the model for one- and two-stage subsurface vertical flow CWs. Model simulations for water flow, tracer transport, and selected biochemical compounds were compared against experimental observations. We also discussed limitations of the model and needs for further improvements of the model. The research needs are largely those listed by Langergraber (2003), but were updated to reflect results of our recent investigations.

CW2D/HYDRUS-2D was developed to model processes in subsurface flow CWs intended for wastewater treatment. However, the model has already been successfully applied to CWs for the treatment of combined sewer overflows (Dittmer et al., 2004). Possible future applications include (i) CWs for other types of wastewater, such as industrial and agricultural wastewaters having compositions similar to municipal wastewater; (ii) infiltration of treated (waste)water into soil for groundwater recharge or other applications; (iii) C and N cycling in soils; (iv) processes in natural wetlands; and (v) processes in riparian areas.


   ACKNOWLEDGMENTS
 
The senior author was partly funded by the Austrian research project "Characterisation of Microbial Biocoenosis to Optimise Removal Efficiency and Design of Subsurface Flow Constructed Wetlands for Wastewater Treatment" (funded by the Austrian Science Fund, FWF, Project P16212-B06) and partly by the Agricultural Experiment Station of the University of California Riverside. The work of Dr. Simunek was supported in part by SAHRA (Sustainability of Semi-Arid Hydrology and Riparian Areas) under the STC Program of the National Science Foundation, Agreement EAR-9876800 and the Terrestrial Sciences Program of the Army Research Office (Terrestrial Processes and Landscape Dynamics and Terrestrial System Modeling and Model Integration).


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 THEORY
 APPLICATIONS TO PILOT-SCALE...
 DISCUSSION
 SUMMARY AND CONCLUSIONS
 REFERENCES