- Email: [email protected]
Diphasic modelling of vertical flow filter
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/ecoleng
Diphasic modelling of vertical flow filter Nicolas Forquet ∗ , Adrien Wanko, Robert Mosé, Antoine-Georges Sadowski Laboratoire Systèmes Hydrauliques Ubrains, Ecole Nationale du Génie de l’Eau et de l’Environnement de Strasbourg, 1 quai Koch, B.P. 61039, 67070 Strasbourg Cedex, France
a r t i c l e
i n f o
a b s t r a c t
Article history:
Wastewater treatment quality in vertical flow filter fed intermittently depends on the oxygen
Received 25 April 2008
renewal capacities. Convective transport by the air phase is one source. Diphasic modelling
Received in revised form
has been shown to be an efficient tool to describe complicated movement of air and water
14 August 2008
during a feeding sequence and thus to estimate the convective dioxygen income. However,
Accepted 11 September 2008
earlier numerical works do not describe the full complexity of the system, like the seepage front existing at the interface with the drain, and are limited to one scenario. In this paper, we introduce an adaptation of an existing diphasic numerical model to the specificity of ver-
Keywords:
tical flow filter with a particular emphasis on the choice of boundary and initial conditions.
Diphasic modelling
Completely different behaviour was observed as ponding does or does not occur. Taking into
Convection
account the influence of air also enables us to understand a phenomenon like early seepage
Dioxygen
at the bottom due to air pushed downward and reduction of infiltration speed linked to air
Hydraulic parameters
compression in the top layers. We investigate relationships between the feeding parameters and the quantity of dioxygen entering the porous media by convection. Hydraulic loading rate does not seem to affect the quantity of dioxygen entering by convection. Increase of the daily hydraulic augments the dioxygen income. However, this phenomenon seems limited as it tends towards an asymptotic value. The number of flushes per day has a contradictory impact: it reduces the income of dioxygen per flush but have a positive impact over a day. Finally we initiate a benchmark of simulation scenarios, which might be a basis for model comparison and experimental validation. © 2008 Elsevier B.V. All rights reserved.
1.
Introduction
Efficient wastewater treatment for small communities and on-site systems is an on-going challenge in Europe, especially in France. Eighty-seven percent of French communities have less than 2000 inhabitants and together represent 15 million inhabitants (IFEN, 2006). Some of them are located on the top of catchment areas where the river ecosystem is particularly sensitive to pollution by ammonium. About 11 million inhabitants make use of on-site wastewater treatment (IFEN, 2006). Badly designed and poorly operating, many of them exhibit poor treatment performance and are a threat for groundwater. Thus, two current challenges are: cutting
∗
Corresponding author. Tel.: +33 388248271; fax: +33 388248283. E-mail address: [email protected] (N. Forquet). 0925-8574/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.ecoleng.2008.09.003
ammonium discharge from small wastewater treatment plant and reducing risks of diffuse pollution from on-site treatment systems. Vertical flow sand filter (VFSF) and vertical flow constructed wetland (VFCW) fed intermittently are among the most promising technologies to reach these goals. The survey studies done over the 20 last years (Brix and Arias, 2005; Molle et al., 2005; Cooper and Green, 1995) and pilot-scale experiments (Prochaska et al., 2007; Wanko et al., 2005) exhibit good efficiency to degrade ammonium and constituents that create biochemical oxygen demand (BOD) mainly because aerobic conditions maintained within the porous media. These systems are robust and require little maintenance. Molle et al. (2006) found that VFCW can handle
48
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
hydraulic overloads during rain events. However, phenomena occurring within these systems are not well understood. Thus, design criteria remain mainly empirical. As highlighted before, good performance depends on aerobic conditions: dioxygen is needed for removal of organic matter (e.g. BOD). To achieve a high ammonium removal, the dioxygen condition in the filter needs to be even better. Then dioxygen renewal is a keystone of any attempt to model transport and degradation of pollutants. We distinguish between dioxygen carried by convection and diffusion. Based on experimental data, Platzer and Mauch (1997) estimate the diffusion transfer rate at 1 gO2 m−2 h−1 except during 1.5 h per day when the soil is too saturated (typically after a flush). They assume that an equivalent volume of air replaces the water volume passing through the filter in order to calculate the amount of dioxygen entering by convection. Kayser and Kunst (2005) compare experimental mass balances with the theoretical approach of Platzer and conclude that the latter overestimates the influence of convection on the total dioxygen transfer. Platzer’s hypothesis considers plug flow for both air and water phases. This is inappropriate during flush especially when ponding occurs on the top. Vachaud et al. (1974), as well as Touma and Vauclin (1986), experimentally describe the situation where air cannot escape at the top and the subsequent compression of the air phase in the porous media. It not only changes the direction of air flow but also influences water infiltration speed. Their work was confirmed by the numerical studies of Phuc and Morel-Seytroux (1972) and Celia and Binning (1992) who used two-phase flow models. In the field of VFSF and VFCW, Schwager and Boller (1997) were the first to introduce two-phase flow modelling in order to represent the air phase velocity during a loading sequence (the flush duration plus the time before the next one) and to estimate dioxygen transfer. They conclude that diffusion is the dominant phenomenon (70% of dioxygen income). The potentials of their approach are great but seem to have been largely unexploited. The aim of this article is to continue the descriptive work initiated by Schwager and Boller in a more complex and realistic configuration and to complete it by a study of the evolution of the quantity of air renewed at each loading sequence as a function of different operating parameters (daily load, number of flush per day and loading rate). Finally, we introduce benchmark scenarios which might be used in future work to test the capability of the model to represent reality by comparisons with experimental data and with simulations by other models.
Table 1 – List of model constants and their values. Name
Value
Specific storativity, Ssw (m−1 ) Water dynamic viscosity, w (kg m−1 s−1 ) Air dynamic viscosity, a (kg m−1 s−1 ) Air density at the reference pressure, 0a (kg m−3 ) Air reference pressure head, h0a (m) Water density, 0w (kg m−3 ) Air–water superficial tension, w (N m−1 ) Acceleration of gravity, g (m s−2 )
1 × 10−4 1 × 10−3 1.82.10−5 1.2 10.1325 1000 7.3 × 10−2 9.81
where w is the water content (m3 m−3 ), Ssw is the specific storativity (m−1 ), hw and ha are the pressure heads for water and air, respectively (mH2 O), is the porosity, krw and kra are the saturation-dependent relative conductivities for water and air, respectively, k is the intrinsic permeability(m2 ), w and a are the dynamic viscosities of water and air, respectively (kg m−1 s−1 ), 0a is the air density at the reference pressure (kg m−3 ), h0a is the reference pressure (mH2 O), a and 0w are the air and water densities, respectively (kg m−3 ). The specific storativity term represents the volume of water released from storage per unit decline in hydraulic head. This usually neglected term has been included here to ensure the validity of Eq. (1) in the saturated zone. The storage effect is the result of a slight compression of the water. To account for air compressibility, air density is expressed as a function of air pressure (Touma and Vauclin, 1986):
a = 0a 1 +
Model description
We used the one-dimensional two-phase flow numerical model developed by Binning (1994) for the simulations. Combining mass conservation and Darcy velocity equations for each phase leads to a set of two equations which have to be solved simultaneously:
⎧ ⎨
−
ws wr + wr 1 + (˛hc )n 1−1/n
(1)
(3)
where ws is the water content at saturation (m3 m−3 ), wr is the residual water content (m3 m−3 ) and ˛ (m−1 ) and n are fitting parameters. In association with the Mualem (1973) theory on unsaturated conductivity, the latter relationship gives:
1/m m
krw (w ) = e 1/2 1 − (1 − e kra (w ) = (1 − e )
krw (w ) k ∂hw w ∂hw ∂ ∂w =0 + Ssw − −1 ∂t ∂t ∂z w ∂z k k ( ) ∂h ∂ ∂h ∂ ra w 0a a w a a ⎩ ( − ) a =0 − a − − ∂t ∂z a ∂z 0w h0a ∂t
(2)
Table 1 summarises values of model constants. The model has been implemented in high-level language to take advantage of vectorization. It allowed us to use smaller time step than in the original model by Binning, which is of great interest while modelling quick variations at the beginning and end of flushes. Fig. 1 shows the modelled filter. It consists of 50 cm of coarse sand over 10 cm of gravel. Water is injected on the top and drains freely at the bottom. Set of Eq. (1) requires specifying three soil-dependent relationships: the pressuresaturation and the relative permeability-saturation for air and for water. The capillary pressure head hc = ha − hw is related to the water content by (Van Genuchten, 1980): w =
2.
ha h0a
1/2
2
)
1/m 2m
1 − e
(4)
where m = 1−1/n and e = ( w − wr )/( ws − wr ). For coarse sand, we used experimental data collected on an alluvial sand, which fits the French guidelines for
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
49
Particular attention was given to the bottom boundary condition. The drain exposes it to atmospheric pressure and thus disrupted the pressure gradient (Abdou and Flury, 2004). As a result, a water-saturated zone must form at the bottom before water can escape. This boundary condition is known as seepage face. Downward seepage occurs only when the gravity force compensates the capillary retention potential: hc,N+1 = 0
(5)
where hc,N+1 is the capillary pressure at the last node. Since gravel and coarse sand have show a large difference in their constitutive relations, the capillary pressure will be discontinuous at the interface. A phenomenon similar to the seepage front occurs and leads to the development of a saturated layer at the bottom of sand. Numerical modelling is also highly sensitive to the choice of initial conditions. We started from dry material with a uniform water pressure head of −15 cm and applied a succession of loading sequence to reach a so-called pseudo-steady-state where variation of water content during a loading sequence is similar to the one at the previous loading sequence (see Fig. 2).
3.
Fig. 1 – Layout of the modelled vertical flow filter.
VFSF and VFCW (Liénard et al., 2001). For this purpose two experiments were done: a sand box experiment to obtain capillary-saturation fitting parameters and a measurement of the saturated conductivity on a small laboratory column to obtain the intrinsic permeability. Results are presented in Table 2 along with the values taken in literature for the gravel.
Table 2 – Sand and gravel parameter values. Data for sand come from experiments done on an alluvial sand collected near Strasbourg. Gravel data have been taken from Meyer and Serne (1999). Sand ws (water content at saturation) (m3 m−3 ) wr (residual water content) (m3 m−3 ) ˛ (fitting parameter) (m−1 ) N (fitting parameter) K (intrinsic permeability) (m2 ) ˚ (porosity) d10 (grain size that is 10% finer by weight) (mm) d60 (grain size that is 60% finer by weight) (mm) Density
Gravel
0.3276
0.42
0.02
0.005
0.03705 3.371 3.3078 × 10−11 0.35 0.21
493 2.29 1 × 10−7 0.45 –
0.57
–
1.7
–
Effects of air entrapment
This section is devoted to the description of air and water phase movements during a loading sequence after the pseudo-steady-state has been reached. We compare a case (1) where ponding does not occur with a case (2) where it does. In case 1, the loading rate is 0.5 m3 m−2 h−1 and the duration of the flush is 180 s. The loading rate does not exceed the infiltration capacities of the sand and the surface remains free of water. The left part of Fig. 3 shows the evolution of water content along a vertical versus time. On the top of the domain, water content increases until feeding stops, which is symbolised on the figure by the thick line labelled “End of flush”. Isolines of saturation give information on the saturation front progression. The variation on water content is smoother as we go deeper in the porous media. At the bottom of the figure, there are two narrow dark areas at depth 0.1 and 0 m which corresponds to high water content. It is consistent with the description made in the previous section of the seepage front at the interface between coarse sand and gravel and at the bottom boundary condition. On the right side of Fig. 3 are plotted three different air Darcian velocity profiles taken at three remarkable times. t1 is located after the flush starts. Air flows in two opposite directions. Above 0.52 m high positive velocities illustrates a fast upward movement of the air phase, which is escaping the porous media at the top. Under 0.52 m, the velocities are slightly negative, which means that only a small proportion of the air volume is forced downward by the incoming water front. At time t2 the flush is over. Air velocity on the top is falling to zero while an upward convective peak is still in action deeper in the porous media. Within a few seconds, the velocities on the top become negative, which means that air is drawn up at the same time that saturation decreases. This is the active step of air renewal as illustrated at time t3 . Velocities are much smaller than during the flush but keeps on until the end of the resting period. An integration of
50
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Fig. 2 – Evolution of water content along a vertical versus simulation time. Hydraulic loading rate = 0.5 m3 m−2 h−1 ; number of load per day = 12; daily hydraulic load = 0.3 m m−2 .
Fig. 3 – On the left is plotted the evolution of the water content along a vertical versus time during a loading sequence. On the right side, three velocity profiles are displayed. They are taken at time t1 , t2 and t3 spotted on the left side figure by dotted lines. Enlargements displayed within the main axes are available for two figures on the right. They correspond to depths where velocity changes sign.
51
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Fig. 4 – On the left is plotted the evolution of the water content along a vertical versus time during a loading sequence. On the right side, three velocity profiles are displayed. They are taken at time t1 , t2 and t3 spotted on the left side figure by dotted lines. We also had on the top right axes the air pressure head profile in grey line.
volumes entering and escaping the filter from the top leads to a null balance. In case 2, the loading rate is 3 m3 m−2 h−1 and flush duration is 100 s. The loading rate is greater that the saturated conductivity (1.17 m3 m−2 h−1 ) and the first centimetres of the soil are quickly saturated. Then water starts to pond on the surface, disenabling air to escape through the surface. Numerically, this is done by switching top boundary conditions for air and water from constant pressure (ha,top = 0 m), also called Dirichlet condition, to constant flux (Qa,top = 0 m s−1 ), called Neumann condition, and from Neumann (Qw,top = Qloading ) to Dirichlet (hw,top n+1 = hw,top n + t Qloading ), respectively. Table 3 synthesizes the evolutions of boundary for both air and water phases. Air pressure builds up quickly in the first centimetres until the capillary pressure at the top reaches the bubbling pressure (also called air-entry pressure), we estimated using:
hae =
w w gd10
(6)
where w is air–water superficial tension (N m−1 ), g is the acceleration due to gravity (m s−2 ) and d10 is the grain size that is 10% finer by weight. At that time, air pressure inside the porous media is sufficient to desaturate the first centimetres and allows air to escape through the surface. Once more we switch the top air boundary condition and set it up to Dirichlet (ha,top = hae + hw ). The right top part of Fig. 4 gives the air pressure head and air Darcian velocity profiles. On the top of the sand, positive air velocities illustrate this movement upward but it is also noticeable that it is limited in the first centimetres. As we go deeper the airflow direction changes and we observe a downward flow. The infiltration front creates a high pressure point in the air phase. Above this point air flows upward, while under it air flows downward. The air pressure head profile illustrates this phenomenon. If we calculate the ratio between the maximum velocity of the downward flow and the maximum velocity of the upward flow (≈1.4) and compare it with the one obtained in the first case (≈0.03) we conclude that ponding increases the velocity of air pushed downward (Wanko, 2005). An immediate consequence of this
Table 3 – Evolution of the boundary conditions for air and water phases during a feeding sequence. Events
Flush starts Ponding starts Bubbling pressure reaches Air escapes stops Ponding ends
Air
Water
Top
Bottom
Top
Bottom
ha = 0 cm Qa = 0 Ha = hae + hw Qa = 0 m s−1 ha = 0 cm
ha = 0.13 cm ha = 0.13 cm ha = 0.13 cm ha = 0.13 cm ha = 0.13 cm
Qw = Qflush hw n+1 = hw n + t Qin hw n+1 = hw n + t Qin hw n+1 = hw n + t Qin Qw = 0 m s−1
hw = ha − hc hw = ha − hc hw = ha − hc hw = ha − hc hw = ha − hc
52
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Fig. 5 – Evolution of the water Darcian velocity on the top of the filter versus time. On the left is case 1 and on the right is case 2.
phenomenon is observed on the left side of Fig. 4. As the ponding starts, we observe on the deep levels of the sand (between depth 0.35 and 0.1) a decrease in the water content. Air forced downward by the incoming water presses on the water which accumulates at the interface between sand and gravel and makes it percolates. Air pressure decreases when most of the air has escaped the porous media and no more air flows through the surface. Then, air top boundary condition is moved back to a no-flux Neumann boundary condition. The new velocity profile is given by the middle graph on the left side of Fig. 4. At that time air is exclusively pushed downward. At the end of the ponding, air starts to be drawn up and this behaviour remains until the end of the resting period. The top boundary conditions are now atmospheric pressure for air phase and no flux for water. Over a loading sequence, the quantity of air entering by the top is 3.57 times greater than the quantity escaping by the top. The difference escapes through the drain. For both cases, we also calculate the ratio between the volume of air entering and the volume of supplied water. In case 1, ratio is equal to 0.99. Air pressure only slightly increases (limited compression) and no air is trapped within the system. The same volume of air is needed to replace the infiltrated water. In case 2, ratio is equal to 0.64. Air is trapped and stored by compression during ponding. When the latter stops, decompression leads to a redistribution of the stored mass of air which occupies empty spaces led by drained water. Thus, ponding reduces the volume of air drawn
up. Two-phase flow modelling not only allows us to have a better understanding of the air movement during the flush but it also allow us to better describe the water flow. Celia and Binning (1992) point out that air trapped under the surface during ponding reduces significantly the infiltration capacities. Fig. 5 compares the water Darcian velocities: when ponding occurs velocity is reduced by 55%. In their paper, Schwager and Boller (1997) estimate the dioxygen income by convection to 0.013 kg m−2 per day for and daily hydraulic load of 0.12 m. We obtain a daily income of 0.029 kg m−2 with a similar load of 0.1 m.
4. Dependency of dioxygen renewal by convection on hydraulic loading parameters We now investigate the dependency of dioxygen renewal by convection on the hydraulic loading parameters: daily hydraulic load, hydraulic loading rate and the number of flush. Of course dioxygen renewal not only depends on these parameters but we choose to focus on them since they may be easily changed by the operator even after the system has been built in order to optimize the quantity of dioxygen taken in the filter by convection. French guidelines (Molle et al., 2005) advise a hydraulic loading rate superior or equal to 0.5 m3 m−2 h−1 , a hydraulic load of 0.3 m per day distributed in about 6 flushes. We compute the mass of dioxygen entering from the top
53
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Table 4 – Benchmark scenarios for validation and models comparison. Benchmark 1 2 3 4 5 6 7 8 9 10 11 12 13
Daily hydraulic load (m) 0.1 0.3 0.3 0.3 0.3 0.3 0.5 0.5 0.5 0.5 0.7 0.7 0.9
Hydraulic loading rate (m s−1 )
Number of load
Ponding
1.39 × 10−4 1.39 × 10−4 1.39 × 10−4 2.78 × 10−4 1.39 × 10−4 2.78 × 10−4 2.78 × 10−4 8.33 × 10−4 2.78 × 10−4 8.33 × 10−4 8.33 × 10−4 8.33 × 10−4 8.33 × 10−4
4 12 4 4 6 6 4 6 8 8 10 12 12
no no no yes no yes yes yes yes yes yes yes yes
during a loading sequence by:
O2 vair (t, z = 0) dt
(7)
6. where fO2 is the bulk fraction of dioxygen in the atmosphere (20.95%) and O2 is the density of dioxygen. We ran a set of simulations where two parameters remain unchanged while the third varies over its entire domain. Daily hydraulic load and the number of flush per day seem to have strong influences on dioxygen income, while hydraulic loading a reduced rate one.
5.
0.029 0.085 0.062 0.048 0.082 0.074 0.069 0.099 0.117 0.120 0.155 0.168 0.189
mation of preferential flow but it does not affect the quantity of dioxygen entering by air convection.
mO2 = fO2
Dioxygen intake per day (kg m−2 )
Hydraulic loading rate
Fig. 6 depicts the evolution of incoming mass of dioxygen per flush and per day versus hydraulic loading rate for three chosen scenarios. The incoming mass of dioxygen per day remains constant in the three scenarios tested. A sufficient loading rate is necessary to allow a good distribution of the influent over the entire surface of the filter and to avoid for-
Daily hydraulic load
Influence of the daily hydraulic load is more predictable (Fig. 7): there is more dioxygen entering the filter as the daily hydraulic load rises. However, curves tend towards asymptotic values. Increasing the daily hydraulic load while keeping constant the hydraulic loading rate and the number of flush necessitates to increase the loading time. It leads to a situation where drainage occurs before loading stops. Then the flow rate at the bottom increases until it equals the loading rate. A steady state is reached and no more air is drawn up. This explains the asymptotic trend. For an hydraulic loading rate of 0.5 m3 m−2 h−1 and four flushes per day, this value is reached quickly and the total increase of the incoming mass of dioxygen per day when passing from 0.1 to 0.9 m is relatively limited: 0.042 kg m−2 . For an hydraulic loading rate of 3 m3 m−2 h−1 and 12 flushes per day, curve inflexion occurs latter (about 0.7 m) and the total increase of incoming mass of dioxygen per day
Fig. 6 – Incoming mass of dioxygen per day (right) and per flush (left) versus hydraulic loading rate. HLR = hydraulic loading rate; DL = daily load and NL = number of load.
54
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Fig. 7 – Incoming mass of dioxygen per day (right) and per flush (left) versus daily hydraulic load. HLR = hydraulic loading rate; DL = daily load and NL = number of load.
is much larger: 0.17 kg m−2 . Thus, increasing the loading time reduces the dioxygen intake if it exceeds the water residence time.
7.
more acute for high daily loading than for small ones where the two phenomena compensate. For a daily loading of 0.9 m and an hydraulic loading rate of 3 m3 m−2 h−1 the incoming mass of dioxygen increases from 0.072 kg m−2 to 0.19 kg m−2 when the number of flushes increases from 4 to 12.
Flush frequency
Incoming mass of dioxygen per flush and incoming mass of dioxygen per day versus the number of flush have opposite trends (see Fig. 8). The more frequent the flush, the less time sand has to drain out the water. A wet sand allows less dioxygen to enter by convection than does a dry sand. This is why curves representing the incoming mass of dioxygen per flush are decreasing as the number of flush augments. However, we notice that this decrease is less than the benefit brings by more frequent flushes. Increasing the number of flushes is a solution to the problem described in the previous section. This is
8. A validation benchmark dedicated to vertical flow filter The idea of development a benchmark grew up as we started to exchange results with other teams. Basically, benchmarks are required at two steps of a code development process: verification and validation. “Verification is the process of determining that a model implementation accurately represents the developer’s conceptual description of the model and the solution to the model. While validation is the process of determining the degree to which a model is an accurate representation of
Fig. 8 – Incoming mass of dioxygen per day (right) and per flush (left) versus daily hydraulic load. HLR = hydraulic loading rate; DL = daily load and NL = number of load.
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
the real world from the perspective of the intended uses of the model.” (AIAA, 1998) A verification test consists of comparing the numerical results with a benchmark of a correct answer provided by highly accurate solutions (Oberkampf and Trucano, 2002) like analytical solutions. In the field of twophase flow modelling an important work has been carried out ´ et al. (2007) to this by McWhorter and Sunada (1990) and Fucik end. A validation test relies on the usage of a benchmark of experimental data or simulation done with other numerical models. Binning and Celia (1992) compare their model with results obtained by Touma and Vauclin (1986) who carried out experiments on a bounded column. This may allow us to verify the ability of the model to take into account ponding but may not be verify the modelling of seepage front. We chose 13 different combinations of the loading parameters in order to obtain a test benchmark. We provide all the information necessary for modelling together with the results obtained with Binning’s model (see Table 4).
9.
Summary, conclusions and future work
Using numerical modelling, we estimate the dioxygen income by convection in vertical flow filter. In comparison with previous work, we introduce more realistic boundary conditions and we run simulation with various hydraulic loading parameters. Guidelines often advise to operate filters with a sufficient loading rate in order to pond over the entire surface. We show that ponding not only improves influent repartition but also allows air to penetrate deeper within the porous media. With the configuration chosen and the associated boundary condition, we found out that air may escape through the drain but hardly enters because of the seepage front between coarse sand and gravel. After a descriptive analysis we focus on the dependency of dioxygen renewal on hydraulic parameters. Hydraulic loading rate does not significantly affect the quantity of dioxygen entering the filter. Increasing the loading time reduces the dioxygen intake if the loading time; exceeds the residence time this problem may be overcome by increasing the number of flushes. Acknowledging the lack of data available for comparison, we introduce a benchmark of simulations for two-phase flow modelling of vertical flow filter. Numerical modelling is becoming more and more popular in the field of constructed wetland and especially for vertical flow filters. This is due to the recent development of a series of code which associate robust monophasic modelling with advanced pollution degradation models (Tomenko et al., 2007; Wanko et al., 2006; Langergraber and Simunek, 2005). While two phase modelling is complementary to this, it is less robust and requires more computing capacities (inadequate time or space step choice might lead the iterative procedure to not converge). Nonetheless it allows to describe a phenomenon that monophasic cannot reproduce and enable us to estimate the error produced by the latter. The two-phase flow model we use in this paper presents some serious limitations we would like to work out. Twodimensional modelling is required to investigate some phenomena like preferential flows or to test some engineered design choice. In order to make the model usable by more scientists we are improving the numerical convergence and
55
computing speed by the combined usage of discontinuous and mixed finite elements (Nayagum et al., 2004; Mosé et al., 1994). Finally, coupling transport-degradation equations for dioxygen will also allow taking into account gaseous diffusion and state, noting its relative importance to convection in the different simulations introduced here.
references
Abdou, H., Flury, M., 2004. Simulation of water flow and solute transport in free-drainage lysimeters and fields soils with heterogeneous structures. Eur. J. Soil Sci. 55, 229–241. AIAA, 1998. Guide for the Verification and Validation of Computational Fluid Dynamics Simulations. American Institute of Aeronautics and Astronautics, Reston, VA, AIAA-G-077-1998. Binning, P., 1994. Modeling unsaturated zone flow and contaminant transport in the air and water phases. Ph.D. Thesis at Princeton University. Brix, H., Arias, C., 2005. The use of vertical flow wetlands for on-site treatment of domestic wastewater: new Danish guidelines. Ecol. Eng. 25, 491–500. Celia, M., Binning, P., 1992. A mass conservative numerical solution for two-phase flow in porous media with application to unsaturated flow. Water Resour. Res. 28, 2819–2828. Cooper, P., Green, B., 1995. Reed bed treatment systems for sewage treatment in the United Kingdom—the first 10 year’s experience. Water Sci. Technol. 32, 317–327. ´ Fucik, R., Mikyska, J., Benes, M., Illangasekare, T., 2007. An improved semi-analytical solution for verification of numerical models of two-phase flow in porous media. Vadose Zone J. 6, 93–104. IFEN, 2006. L’assainissement en France en 1998 et 2001. Dossier de l’Ifen no D3, February 2006 [in French]. Kayser, K., Kunst, S., 2005. Processes in vertical-flow reed beds: nitrification oxygen transfer and soil clogging. Water Sci. Technol. 51, 177–184. Langergraber, G., Simunek, J., 2005. Modeling variably saturated flow and multicomponent reactive transport in constructed Wetlands. Vadose Zone J. 4, 924–938. Liénard, A., Guellaf, H., Boutin, C., 2001. Choice of the sand material for sand filters used for secondary treatment of wastewater. Water Sci. Technol. 44, 189–196. McWhorter, D., Sunada, D., 1990. Exact integral solutions for two-phase flow. Water Resour. Res. 26, 399–413. Meyer, P.D., Serne, R.J., 1999. Near-Field Hydrology Data Package for Immobilized Low-Activity Waste 2001 Performance Assessment. Pacific Northwest national Laboratory, Richland, Washington. Molle, P., Liénard, A., Boutin, C., Merlin, G., Iwema, A., 2005. How to treat raw sewage with constructed wetlands: an overview of the French Systems. Water Sci. Technol. 51, 11–21. Molle, P., Liénard, A., Grasmick, A., Iwema, A., 2006. Effect of reeds and feeding operations on hydraulic behaviour of vertical flow constructed wetlands under hydraulic overloads. Water Res. 10, 606–612. Mosé, R., Siegel, P., Ackerer, P., Chavent, G., 1994. Application of the mixed hybrid finite element approximation in a groundwater flow model: luxury or necessity? Water Resour. Res. 30, 3001–3012. Mualem, Y., 1973. A new model of predicting the hydraulic conductivity of porous media. Water Resour. Res. 12, 513–522. Nayagum, D., Schäfer, G., Mosé, R., 2004. Modeling two-phase incompressible flow in porous media using mixed hybrid and discontinuous finite elements. Comput. Geosci. 8, 49–73.
56
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Oberkampf, W.L., Trucano, T.G., 2002. Verification and validation in computational fluid dynamics. Prog. Aerospace Sci. 38, 209–272. Phuc, L., Morel-Seytroux, H., 1972. Effect of soil on air movement and compressibility on infiltration events. Soil Sci. Am. Proc. 36, 237–241. Platzer, C., Mauch, K., 1997. Soil clogging in vertical flow reed beds mechanisms, parameters, consequences and solutions? Water Sci. Technol. 35, 175–181. Prochaska, C.A., Zouboulis, A.I., Eskridge, K.M., 2007. Performance of pilot-scale vertical-flow constructed wetlands, as affected by season, substrate, hydraulic load and frequency of application of simulated urban sewage. Ecol. Eng. 31, 57–66. Schwager, A., Boller, M., 1997. Transport phenomena in intermittent filters. Water Sci. Technol. 35, 13–20. Tomenko, V., Ahmed, S., Popov, V., 2007. Modelling constructed wetland treatment system performance. Ecol. Model. 205, 355–364. Touma, J., Vauclin, M., 1986. Experimental and numerical analysis of two-phase infiltration in a partially saturated soil. Transp. Porous Media 1, 22–55.
Vachaud, G., Gaudet, J., Kuraz, V., 1974. Air and water flow during ponded infiltration in a vertical bounded column of soil. J. Hydrol. 22, 89–108. Van Genuchten, M.Th., 1980. A closed form equation for predicting the hydraulic conductivity in soils. Soil Sci. Soc. Am. J. 44, 892–898. Wanko, A., 2005. Etude des mécanismes de transfert et évaluation des capacités d’oxygénation et de traitement des dispositifs de traitement par cultures fixées sur supports fins alimentés en discontinu. Développement d’un modèle de biodégradation. Ph.D. Thesis. University Louis Pasteur, Strasbourg [in French]. Wanko, A., Mosé, R., Beck, C., 2005. Biological processing capacities and biomass growth in waste water treatment by infiltration on two kinds of sand. Water Air Soil Pollut. 165, 279–299. Wanko, A., Mosé, R., Carrayrou, J., Sadowski, A.G., 2006. Simulation of biodegradation in infiltration seepage, model development and hydrodynamic calibration. Water, Air Soil Pollut. 177, 19–43.
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/ecoleng
Diphasic modelling of vertical flow filter Nicolas Forquet ∗ , Adrien Wanko, Robert Mosé, Antoine-Georges Sadowski Laboratoire Systèmes Hydrauliques Ubrains, Ecole Nationale du Génie de l’Eau et de l’Environnement de Strasbourg, 1 quai Koch, B.P. 61039, 67070 Strasbourg Cedex, France
a r t i c l e
i n f o
a b s t r a c t
Article history:
Wastewater treatment quality in vertical flow filter fed intermittently depends on the oxygen
Received 25 April 2008
renewal capacities. Convective transport by the air phase is one source. Diphasic modelling
Received in revised form
has been shown to be an efficient tool to describe complicated movement of air and water
14 August 2008
during a feeding sequence and thus to estimate the convective dioxygen income. However,
Accepted 11 September 2008
earlier numerical works do not describe the full complexity of the system, like the seepage front existing at the interface with the drain, and are limited to one scenario. In this paper, we introduce an adaptation of an existing diphasic numerical model to the specificity of ver-
Keywords:
tical flow filter with a particular emphasis on the choice of boundary and initial conditions.
Diphasic modelling
Completely different behaviour was observed as ponding does or does not occur. Taking into
Convection
account the influence of air also enables us to understand a phenomenon like early seepage
Dioxygen
at the bottom due to air pushed downward and reduction of infiltration speed linked to air
Hydraulic parameters
compression in the top layers. We investigate relationships between the feeding parameters and the quantity of dioxygen entering the porous media by convection. Hydraulic loading rate does not seem to affect the quantity of dioxygen entering by convection. Increase of the daily hydraulic augments the dioxygen income. However, this phenomenon seems limited as it tends towards an asymptotic value. The number of flushes per day has a contradictory impact: it reduces the income of dioxygen per flush but have a positive impact over a day. Finally we initiate a benchmark of simulation scenarios, which might be a basis for model comparison and experimental validation. © 2008 Elsevier B.V. All rights reserved.
1.
Introduction
Efficient wastewater treatment for small communities and on-site systems is an on-going challenge in Europe, especially in France. Eighty-seven percent of French communities have less than 2000 inhabitants and together represent 15 million inhabitants (IFEN, 2006). Some of them are located on the top of catchment areas where the river ecosystem is particularly sensitive to pollution by ammonium. About 11 million inhabitants make use of on-site wastewater treatment (IFEN, 2006). Badly designed and poorly operating, many of them exhibit poor treatment performance and are a threat for groundwater. Thus, two current challenges are: cutting
∗
Corresponding author. Tel.: +33 388248271; fax: +33 388248283. E-mail address: [email protected] (N. Forquet). 0925-8574/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.ecoleng.2008.09.003
ammonium discharge from small wastewater treatment plant and reducing risks of diffuse pollution from on-site treatment systems. Vertical flow sand filter (VFSF) and vertical flow constructed wetland (VFCW) fed intermittently are among the most promising technologies to reach these goals. The survey studies done over the 20 last years (Brix and Arias, 2005; Molle et al., 2005; Cooper and Green, 1995) and pilot-scale experiments (Prochaska et al., 2007; Wanko et al., 2005) exhibit good efficiency to degrade ammonium and constituents that create biochemical oxygen demand (BOD) mainly because aerobic conditions maintained within the porous media. These systems are robust and require little maintenance. Molle et al. (2006) found that VFCW can handle
48
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
hydraulic overloads during rain events. However, phenomena occurring within these systems are not well understood. Thus, design criteria remain mainly empirical. As highlighted before, good performance depends on aerobic conditions: dioxygen is needed for removal of organic matter (e.g. BOD). To achieve a high ammonium removal, the dioxygen condition in the filter needs to be even better. Then dioxygen renewal is a keystone of any attempt to model transport and degradation of pollutants. We distinguish between dioxygen carried by convection and diffusion. Based on experimental data, Platzer and Mauch (1997) estimate the diffusion transfer rate at 1 gO2 m−2 h−1 except during 1.5 h per day when the soil is too saturated (typically after a flush). They assume that an equivalent volume of air replaces the water volume passing through the filter in order to calculate the amount of dioxygen entering by convection. Kayser and Kunst (2005) compare experimental mass balances with the theoretical approach of Platzer and conclude that the latter overestimates the influence of convection on the total dioxygen transfer. Platzer’s hypothesis considers plug flow for both air and water phases. This is inappropriate during flush especially when ponding occurs on the top. Vachaud et al. (1974), as well as Touma and Vauclin (1986), experimentally describe the situation where air cannot escape at the top and the subsequent compression of the air phase in the porous media. It not only changes the direction of air flow but also influences water infiltration speed. Their work was confirmed by the numerical studies of Phuc and Morel-Seytroux (1972) and Celia and Binning (1992) who used two-phase flow models. In the field of VFSF and VFCW, Schwager and Boller (1997) were the first to introduce two-phase flow modelling in order to represent the air phase velocity during a loading sequence (the flush duration plus the time before the next one) and to estimate dioxygen transfer. They conclude that diffusion is the dominant phenomenon (70% of dioxygen income). The potentials of their approach are great but seem to have been largely unexploited. The aim of this article is to continue the descriptive work initiated by Schwager and Boller in a more complex and realistic configuration and to complete it by a study of the evolution of the quantity of air renewed at each loading sequence as a function of different operating parameters (daily load, number of flush per day and loading rate). Finally, we introduce benchmark scenarios which might be used in future work to test the capability of the model to represent reality by comparisons with experimental data and with simulations by other models.
Table 1 – List of model constants and their values. Name
Value
Specific storativity, Ssw (m−1 ) Water dynamic viscosity, w (kg m−1 s−1 ) Air dynamic viscosity, a (kg m−1 s−1 ) Air density at the reference pressure, 0a (kg m−3 ) Air reference pressure head, h0a (m) Water density, 0w (kg m−3 ) Air–water superficial tension, w (N m−1 ) Acceleration of gravity, g (m s−2 )
1 × 10−4 1 × 10−3 1.82.10−5 1.2 10.1325 1000 7.3 × 10−2 9.81
where w is the water content (m3 m−3 ), Ssw is the specific storativity (m−1 ), hw and ha are the pressure heads for water and air, respectively (mH2 O), is the porosity, krw and kra are the saturation-dependent relative conductivities for water and air, respectively, k is the intrinsic permeability(m2 ), w and a are the dynamic viscosities of water and air, respectively (kg m−1 s−1 ), 0a is the air density at the reference pressure (kg m−3 ), h0a is the reference pressure (mH2 O), a and 0w are the air and water densities, respectively (kg m−3 ). The specific storativity term represents the volume of water released from storage per unit decline in hydraulic head. This usually neglected term has been included here to ensure the validity of Eq. (1) in the saturated zone. The storage effect is the result of a slight compression of the water. To account for air compressibility, air density is expressed as a function of air pressure (Touma and Vauclin, 1986):
a = 0a 1 +
Model description
We used the one-dimensional two-phase flow numerical model developed by Binning (1994) for the simulations. Combining mass conservation and Darcy velocity equations for each phase leads to a set of two equations which have to be solved simultaneously:
⎧ ⎨
−
ws wr + wr 1 + (˛hc )n 1−1/n
(1)
(3)
where ws is the water content at saturation (m3 m−3 ), wr is the residual water content (m3 m−3 ) and ˛ (m−1 ) and n are fitting parameters. In association with the Mualem (1973) theory on unsaturated conductivity, the latter relationship gives:
1/m m
krw (w ) = e 1/2 1 − (1 − e kra (w ) = (1 − e )
krw (w ) k ∂hw w ∂hw ∂ ∂w =0 + Ssw − −1 ∂t ∂t ∂z w ∂z k k ( ) ∂h ∂ ∂h ∂ ra w 0a a w a a ⎩ ( − ) a =0 − a − − ∂t ∂z a ∂z 0w h0a ∂t
(2)
Table 1 summarises values of model constants. The model has been implemented in high-level language to take advantage of vectorization. It allowed us to use smaller time step than in the original model by Binning, which is of great interest while modelling quick variations at the beginning and end of flushes. Fig. 1 shows the modelled filter. It consists of 50 cm of coarse sand over 10 cm of gravel. Water is injected on the top and drains freely at the bottom. Set of Eq. (1) requires specifying three soil-dependent relationships: the pressuresaturation and the relative permeability-saturation for air and for water. The capillary pressure head hc = ha − hw is related to the water content by (Van Genuchten, 1980): w =
2.
ha h0a
1/2
2
)
1/m 2m
1 − e
(4)
where m = 1−1/n and e = ( w − wr )/( ws − wr ). For coarse sand, we used experimental data collected on an alluvial sand, which fits the French guidelines for
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
49
Particular attention was given to the bottom boundary condition. The drain exposes it to atmospheric pressure and thus disrupted the pressure gradient (Abdou and Flury, 2004). As a result, a water-saturated zone must form at the bottom before water can escape. This boundary condition is known as seepage face. Downward seepage occurs only when the gravity force compensates the capillary retention potential: hc,N+1 = 0
(5)
where hc,N+1 is the capillary pressure at the last node. Since gravel and coarse sand have show a large difference in their constitutive relations, the capillary pressure will be discontinuous at the interface. A phenomenon similar to the seepage front occurs and leads to the development of a saturated layer at the bottom of sand. Numerical modelling is also highly sensitive to the choice of initial conditions. We started from dry material with a uniform water pressure head of −15 cm and applied a succession of loading sequence to reach a so-called pseudo-steady-state where variation of water content during a loading sequence is similar to the one at the previous loading sequence (see Fig. 2).
3.
Fig. 1 – Layout of the modelled vertical flow filter.
VFSF and VFCW (Liénard et al., 2001). For this purpose two experiments were done: a sand box experiment to obtain capillary-saturation fitting parameters and a measurement of the saturated conductivity on a small laboratory column to obtain the intrinsic permeability. Results are presented in Table 2 along with the values taken in literature for the gravel.
Table 2 – Sand and gravel parameter values. Data for sand come from experiments done on an alluvial sand collected near Strasbourg. Gravel data have been taken from Meyer and Serne (1999). Sand ws (water content at saturation) (m3 m−3 ) wr (residual water content) (m3 m−3 ) ˛ (fitting parameter) (m−1 ) N (fitting parameter) K (intrinsic permeability) (m2 ) ˚ (porosity) d10 (grain size that is 10% finer by weight) (mm) d60 (grain size that is 60% finer by weight) (mm) Density
Gravel
0.3276
0.42
0.02
0.005
0.03705 3.371 3.3078 × 10−11 0.35 0.21
493 2.29 1 × 10−7 0.45 –
0.57
–
1.7
–
Effects of air entrapment
This section is devoted to the description of air and water phase movements during a loading sequence after the pseudo-steady-state has been reached. We compare a case (1) where ponding does not occur with a case (2) where it does. In case 1, the loading rate is 0.5 m3 m−2 h−1 and the duration of the flush is 180 s. The loading rate does not exceed the infiltration capacities of the sand and the surface remains free of water. The left part of Fig. 3 shows the evolution of water content along a vertical versus time. On the top of the domain, water content increases until feeding stops, which is symbolised on the figure by the thick line labelled “End of flush”. Isolines of saturation give information on the saturation front progression. The variation on water content is smoother as we go deeper in the porous media. At the bottom of the figure, there are two narrow dark areas at depth 0.1 and 0 m which corresponds to high water content. It is consistent with the description made in the previous section of the seepage front at the interface between coarse sand and gravel and at the bottom boundary condition. On the right side of Fig. 3 are plotted three different air Darcian velocity profiles taken at three remarkable times. t1 is located after the flush starts. Air flows in two opposite directions. Above 0.52 m high positive velocities illustrates a fast upward movement of the air phase, which is escaping the porous media at the top. Under 0.52 m, the velocities are slightly negative, which means that only a small proportion of the air volume is forced downward by the incoming water front. At time t2 the flush is over. Air velocity on the top is falling to zero while an upward convective peak is still in action deeper in the porous media. Within a few seconds, the velocities on the top become negative, which means that air is drawn up at the same time that saturation decreases. This is the active step of air renewal as illustrated at time t3 . Velocities are much smaller than during the flush but keeps on until the end of the resting period. An integration of
50
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Fig. 2 – Evolution of water content along a vertical versus simulation time. Hydraulic loading rate = 0.5 m3 m−2 h−1 ; number of load per day = 12; daily hydraulic load = 0.3 m m−2 .
Fig. 3 – On the left is plotted the evolution of the water content along a vertical versus time during a loading sequence. On the right side, three velocity profiles are displayed. They are taken at time t1 , t2 and t3 spotted on the left side figure by dotted lines. Enlargements displayed within the main axes are available for two figures on the right. They correspond to depths where velocity changes sign.
51
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Fig. 4 – On the left is plotted the evolution of the water content along a vertical versus time during a loading sequence. On the right side, three velocity profiles are displayed. They are taken at time t1 , t2 and t3 spotted on the left side figure by dotted lines. We also had on the top right axes the air pressure head profile in grey line.
volumes entering and escaping the filter from the top leads to a null balance. In case 2, the loading rate is 3 m3 m−2 h−1 and flush duration is 100 s. The loading rate is greater that the saturated conductivity (1.17 m3 m−2 h−1 ) and the first centimetres of the soil are quickly saturated. Then water starts to pond on the surface, disenabling air to escape through the surface. Numerically, this is done by switching top boundary conditions for air and water from constant pressure (ha,top = 0 m), also called Dirichlet condition, to constant flux (Qa,top = 0 m s−1 ), called Neumann condition, and from Neumann (Qw,top = Qloading ) to Dirichlet (hw,top n+1 = hw,top n + t Qloading ), respectively. Table 3 synthesizes the evolutions of boundary for both air and water phases. Air pressure builds up quickly in the first centimetres until the capillary pressure at the top reaches the bubbling pressure (also called air-entry pressure), we estimated using:
hae =
w w gd10
(6)
where w is air–water superficial tension (N m−1 ), g is the acceleration due to gravity (m s−2 ) and d10 is the grain size that is 10% finer by weight. At that time, air pressure inside the porous media is sufficient to desaturate the first centimetres and allows air to escape through the surface. Once more we switch the top air boundary condition and set it up to Dirichlet (ha,top = hae + hw ). The right top part of Fig. 4 gives the air pressure head and air Darcian velocity profiles. On the top of the sand, positive air velocities illustrate this movement upward but it is also noticeable that it is limited in the first centimetres. As we go deeper the airflow direction changes and we observe a downward flow. The infiltration front creates a high pressure point in the air phase. Above this point air flows upward, while under it air flows downward. The air pressure head profile illustrates this phenomenon. If we calculate the ratio between the maximum velocity of the downward flow and the maximum velocity of the upward flow (≈1.4) and compare it with the one obtained in the first case (≈0.03) we conclude that ponding increases the velocity of air pushed downward (Wanko, 2005). An immediate consequence of this
Table 3 – Evolution of the boundary conditions for air and water phases during a feeding sequence. Events
Flush starts Ponding starts Bubbling pressure reaches Air escapes stops Ponding ends
Air
Water
Top
Bottom
Top
Bottom
ha = 0 cm Qa = 0 Ha = hae + hw Qa = 0 m s−1 ha = 0 cm
ha = 0.13 cm ha = 0.13 cm ha = 0.13 cm ha = 0.13 cm ha = 0.13 cm
Qw = Qflush hw n+1 = hw n + t Qin hw n+1 = hw n + t Qin hw n+1 = hw n + t Qin Qw = 0 m s−1
hw = ha − hc hw = ha − hc hw = ha − hc hw = ha − hc hw = ha − hc
52
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Fig. 5 – Evolution of the water Darcian velocity on the top of the filter versus time. On the left is case 1 and on the right is case 2.
phenomenon is observed on the left side of Fig. 4. As the ponding starts, we observe on the deep levels of the sand (between depth 0.35 and 0.1) a decrease in the water content. Air forced downward by the incoming water presses on the water which accumulates at the interface between sand and gravel and makes it percolates. Air pressure decreases when most of the air has escaped the porous media and no more air flows through the surface. Then, air top boundary condition is moved back to a no-flux Neumann boundary condition. The new velocity profile is given by the middle graph on the left side of Fig. 4. At that time air is exclusively pushed downward. At the end of the ponding, air starts to be drawn up and this behaviour remains until the end of the resting period. The top boundary conditions are now atmospheric pressure for air phase and no flux for water. Over a loading sequence, the quantity of air entering by the top is 3.57 times greater than the quantity escaping by the top. The difference escapes through the drain. For both cases, we also calculate the ratio between the volume of air entering and the volume of supplied water. In case 1, ratio is equal to 0.99. Air pressure only slightly increases (limited compression) and no air is trapped within the system. The same volume of air is needed to replace the infiltrated water. In case 2, ratio is equal to 0.64. Air is trapped and stored by compression during ponding. When the latter stops, decompression leads to a redistribution of the stored mass of air which occupies empty spaces led by drained water. Thus, ponding reduces the volume of air drawn
up. Two-phase flow modelling not only allows us to have a better understanding of the air movement during the flush but it also allow us to better describe the water flow. Celia and Binning (1992) point out that air trapped under the surface during ponding reduces significantly the infiltration capacities. Fig. 5 compares the water Darcian velocities: when ponding occurs velocity is reduced by 55%. In their paper, Schwager and Boller (1997) estimate the dioxygen income by convection to 0.013 kg m−2 per day for and daily hydraulic load of 0.12 m. We obtain a daily income of 0.029 kg m−2 with a similar load of 0.1 m.
4. Dependency of dioxygen renewal by convection on hydraulic loading parameters We now investigate the dependency of dioxygen renewal by convection on the hydraulic loading parameters: daily hydraulic load, hydraulic loading rate and the number of flush. Of course dioxygen renewal not only depends on these parameters but we choose to focus on them since they may be easily changed by the operator even after the system has been built in order to optimize the quantity of dioxygen taken in the filter by convection. French guidelines (Molle et al., 2005) advise a hydraulic loading rate superior or equal to 0.5 m3 m−2 h−1 , a hydraulic load of 0.3 m per day distributed in about 6 flushes. We compute the mass of dioxygen entering from the top
53
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Table 4 – Benchmark scenarios for validation and models comparison. Benchmark 1 2 3 4 5 6 7 8 9 10 11 12 13
Daily hydraulic load (m) 0.1 0.3 0.3 0.3 0.3 0.3 0.5 0.5 0.5 0.5 0.7 0.7 0.9
Hydraulic loading rate (m s−1 )
Number of load
Ponding
1.39 × 10−4 1.39 × 10−4 1.39 × 10−4 2.78 × 10−4 1.39 × 10−4 2.78 × 10−4 2.78 × 10−4 8.33 × 10−4 2.78 × 10−4 8.33 × 10−4 8.33 × 10−4 8.33 × 10−4 8.33 × 10−4
4 12 4 4 6 6 4 6 8 8 10 12 12
no no no yes no yes yes yes yes yes yes yes yes
during a loading sequence by:
O2 vair (t, z = 0) dt
(7)
6. where fO2 is the bulk fraction of dioxygen in the atmosphere (20.95%) and O2 is the density of dioxygen. We ran a set of simulations where two parameters remain unchanged while the third varies over its entire domain. Daily hydraulic load and the number of flush per day seem to have strong influences on dioxygen income, while hydraulic loading a reduced rate one.
5.
0.029 0.085 0.062 0.048 0.082 0.074 0.069 0.099 0.117 0.120 0.155 0.168 0.189
mation of preferential flow but it does not affect the quantity of dioxygen entering by air convection.
mO2 = fO2
Dioxygen intake per day (kg m−2 )
Hydraulic loading rate
Fig. 6 depicts the evolution of incoming mass of dioxygen per flush and per day versus hydraulic loading rate for three chosen scenarios. The incoming mass of dioxygen per day remains constant in the three scenarios tested. A sufficient loading rate is necessary to allow a good distribution of the influent over the entire surface of the filter and to avoid for-
Daily hydraulic load
Influence of the daily hydraulic load is more predictable (Fig. 7): there is more dioxygen entering the filter as the daily hydraulic load rises. However, curves tend towards asymptotic values. Increasing the daily hydraulic load while keeping constant the hydraulic loading rate and the number of flush necessitates to increase the loading time. It leads to a situation where drainage occurs before loading stops. Then the flow rate at the bottom increases until it equals the loading rate. A steady state is reached and no more air is drawn up. This explains the asymptotic trend. For an hydraulic loading rate of 0.5 m3 m−2 h−1 and four flushes per day, this value is reached quickly and the total increase of the incoming mass of dioxygen per day when passing from 0.1 to 0.9 m is relatively limited: 0.042 kg m−2 . For an hydraulic loading rate of 3 m3 m−2 h−1 and 12 flushes per day, curve inflexion occurs latter (about 0.7 m) and the total increase of incoming mass of dioxygen per day
Fig. 6 – Incoming mass of dioxygen per day (right) and per flush (left) versus hydraulic loading rate. HLR = hydraulic loading rate; DL = daily load and NL = number of load.
54
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Fig. 7 – Incoming mass of dioxygen per day (right) and per flush (left) versus daily hydraulic load. HLR = hydraulic loading rate; DL = daily load and NL = number of load.
is much larger: 0.17 kg m−2 . Thus, increasing the loading time reduces the dioxygen intake if it exceeds the water residence time.
7.
more acute for high daily loading than for small ones where the two phenomena compensate. For a daily loading of 0.9 m and an hydraulic loading rate of 3 m3 m−2 h−1 the incoming mass of dioxygen increases from 0.072 kg m−2 to 0.19 kg m−2 when the number of flushes increases from 4 to 12.
Flush frequency
Incoming mass of dioxygen per flush and incoming mass of dioxygen per day versus the number of flush have opposite trends (see Fig. 8). The more frequent the flush, the less time sand has to drain out the water. A wet sand allows less dioxygen to enter by convection than does a dry sand. This is why curves representing the incoming mass of dioxygen per flush are decreasing as the number of flush augments. However, we notice that this decrease is less than the benefit brings by more frequent flushes. Increasing the number of flushes is a solution to the problem described in the previous section. This is
8. A validation benchmark dedicated to vertical flow filter The idea of development a benchmark grew up as we started to exchange results with other teams. Basically, benchmarks are required at two steps of a code development process: verification and validation. “Verification is the process of determining that a model implementation accurately represents the developer’s conceptual description of the model and the solution to the model. While validation is the process of determining the degree to which a model is an accurate representation of
Fig. 8 – Incoming mass of dioxygen per day (right) and per flush (left) versus daily hydraulic load. HLR = hydraulic loading rate; DL = daily load and NL = number of load.
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
the real world from the perspective of the intended uses of the model.” (AIAA, 1998) A verification test consists of comparing the numerical results with a benchmark of a correct answer provided by highly accurate solutions (Oberkampf and Trucano, 2002) like analytical solutions. In the field of twophase flow modelling an important work has been carried out ´ et al. (2007) to this by McWhorter and Sunada (1990) and Fucik end. A validation test relies on the usage of a benchmark of experimental data or simulation done with other numerical models. Binning and Celia (1992) compare their model with results obtained by Touma and Vauclin (1986) who carried out experiments on a bounded column. This may allow us to verify the ability of the model to take into account ponding but may not be verify the modelling of seepage front. We chose 13 different combinations of the loading parameters in order to obtain a test benchmark. We provide all the information necessary for modelling together with the results obtained with Binning’s model (see Table 4).
9.
Summary, conclusions and future work
Using numerical modelling, we estimate the dioxygen income by convection in vertical flow filter. In comparison with previous work, we introduce more realistic boundary conditions and we run simulation with various hydraulic loading parameters. Guidelines often advise to operate filters with a sufficient loading rate in order to pond over the entire surface. We show that ponding not only improves influent repartition but also allows air to penetrate deeper within the porous media. With the configuration chosen and the associated boundary condition, we found out that air may escape through the drain but hardly enters because of the seepage front between coarse sand and gravel. After a descriptive analysis we focus on the dependency of dioxygen renewal on hydraulic parameters. Hydraulic loading rate does not significantly affect the quantity of dioxygen entering the filter. Increasing the loading time reduces the dioxygen intake if the loading time; exceeds the residence time this problem may be overcome by increasing the number of flushes. Acknowledging the lack of data available for comparison, we introduce a benchmark of simulations for two-phase flow modelling of vertical flow filter. Numerical modelling is becoming more and more popular in the field of constructed wetland and especially for vertical flow filters. This is due to the recent development of a series of code which associate robust monophasic modelling with advanced pollution degradation models (Tomenko et al., 2007; Wanko et al., 2006; Langergraber and Simunek, 2005). While two phase modelling is complementary to this, it is less robust and requires more computing capacities (inadequate time or space step choice might lead the iterative procedure to not converge). Nonetheless it allows to describe a phenomenon that monophasic cannot reproduce and enable us to estimate the error produced by the latter. The two-phase flow model we use in this paper presents some serious limitations we would like to work out. Twodimensional modelling is required to investigate some phenomena like preferential flows or to test some engineered design choice. In order to make the model usable by more scientists we are improving the numerical convergence and
55
computing speed by the combined usage of discontinuous and mixed finite elements (Nayagum et al., 2004; Mosé et al., 1994). Finally, coupling transport-degradation equations for dioxygen will also allow taking into account gaseous diffusion and state, noting its relative importance to convection in the different simulations introduced here.
references
Abdou, H., Flury, M., 2004. Simulation of water flow and solute transport in free-drainage lysimeters and fields soils with heterogeneous structures. Eur. J. Soil Sci. 55, 229–241. AIAA, 1998. Guide for the Verification and Validation of Computational Fluid Dynamics Simulations. American Institute of Aeronautics and Astronautics, Reston, VA, AIAA-G-077-1998. Binning, P., 1994. Modeling unsaturated zone flow and contaminant transport in the air and water phases. Ph.D. Thesis at Princeton University. Brix, H., Arias, C., 2005. The use of vertical flow wetlands for on-site treatment of domestic wastewater: new Danish guidelines. Ecol. Eng. 25, 491–500. Celia, M., Binning, P., 1992. A mass conservative numerical solution for two-phase flow in porous media with application to unsaturated flow. Water Resour. Res. 28, 2819–2828. Cooper, P., Green, B., 1995. Reed bed treatment systems for sewage treatment in the United Kingdom—the first 10 year’s experience. Water Sci. Technol. 32, 317–327. ´ Fucik, R., Mikyska, J., Benes, M., Illangasekare, T., 2007. An improved semi-analytical solution for verification of numerical models of two-phase flow in porous media. Vadose Zone J. 6, 93–104. IFEN, 2006. L’assainissement en France en 1998 et 2001. Dossier de l’Ifen no D3, February 2006 [in French]. Kayser, K., Kunst, S., 2005. Processes in vertical-flow reed beds: nitrification oxygen transfer and soil clogging. Water Sci. Technol. 51, 177–184. Langergraber, G., Simunek, J., 2005. Modeling variably saturated flow and multicomponent reactive transport in constructed Wetlands. Vadose Zone J. 4, 924–938. Liénard, A., Guellaf, H., Boutin, C., 2001. Choice of the sand material for sand filters used for secondary treatment of wastewater. Water Sci. Technol. 44, 189–196. McWhorter, D., Sunada, D., 1990. Exact integral solutions for two-phase flow. Water Resour. Res. 26, 399–413. Meyer, P.D., Serne, R.J., 1999. Near-Field Hydrology Data Package for Immobilized Low-Activity Waste 2001 Performance Assessment. Pacific Northwest national Laboratory, Richland, Washington. Molle, P., Liénard, A., Boutin, C., Merlin, G., Iwema, A., 2005. How to treat raw sewage with constructed wetlands: an overview of the French Systems. Water Sci. Technol. 51, 11–21. Molle, P., Liénard, A., Grasmick, A., Iwema, A., 2006. Effect of reeds and feeding operations on hydraulic behaviour of vertical flow constructed wetlands under hydraulic overloads. Water Res. 10, 606–612. Mosé, R., Siegel, P., Ackerer, P., Chavent, G., 1994. Application of the mixed hybrid finite element approximation in a groundwater flow model: luxury or necessity? Water Resour. Res. 30, 3001–3012. Mualem, Y., 1973. A new model of predicting the hydraulic conductivity of porous media. Water Resour. Res. 12, 513–522. Nayagum, D., Schäfer, G., Mosé, R., 2004. Modeling two-phase incompressible flow in porous media using mixed hybrid and discontinuous finite elements. Comput. Geosci. 8, 49–73.
56
e c o l o g i c a l e n g i n e e r i n g 3 5 ( 2 0 0 9 ) 47–56
Oberkampf, W.L., Trucano, T.G., 2002. Verification and validation in computational fluid dynamics. Prog. Aerospace Sci. 38, 209–272. Phuc, L., Morel-Seytroux, H., 1972. Effect of soil on air movement and compressibility on infiltration events. Soil Sci. Am. Proc. 36, 237–241. Platzer, C., Mauch, K., 1997. Soil clogging in vertical flow reed beds mechanisms, parameters, consequences and solutions? Water Sci. Technol. 35, 175–181. Prochaska, C.A., Zouboulis, A.I., Eskridge, K.M., 2007. Performance of pilot-scale vertical-flow constructed wetlands, as affected by season, substrate, hydraulic load and frequency of application of simulated urban sewage. Ecol. Eng. 31, 57–66. Schwager, A., Boller, M., 1997. Transport phenomena in intermittent filters. Water Sci. Technol. 35, 13–20. Tomenko, V., Ahmed, S., Popov, V., 2007. Modelling constructed wetland treatment system performance. Ecol. Model. 205, 355–364. Touma, J., Vauclin, M., 1986. Experimental and numerical analysis of two-phase infiltration in a partially saturated soil. Transp. Porous Media 1, 22–55.
Vachaud, G., Gaudet, J., Kuraz, V., 1974. Air and water flow during ponded infiltration in a vertical bounded column of soil. J. Hydrol. 22, 89–108. Van Genuchten, M.Th., 1980. A closed form equation for predicting the hydraulic conductivity in soils. Soil Sci. Soc. Am. J. 44, 892–898. Wanko, A., 2005. Etude des mécanismes de transfert et évaluation des capacités d’oxygénation et de traitement des dispositifs de traitement par cultures fixées sur supports fins alimentés en discontinu. Développement d’un modèle de biodégradation. Ph.D. Thesis. University Louis Pasteur, Strasbourg [in French]. Wanko, A., Mosé, R., Beck, C., 2005. Biological processing capacities and biomass growth in waste water treatment by infiltration on two kinds of sand. Water Air Soil Pollut. 165, 279–299. Wanko, A., Mosé, R., Carrayrou, J., Sadowski, A.G., 2006. Simulation of biodegradation in infiltration seepage, model development and hydrodynamic calibration. Water, Air Soil Pollut. 177, 19–43.