Comparison of steady-state and transient flow conditions on reactive transport of contaminants in the vadose soil zone

Journal of Hydrology 369 (2009) 225–233

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Comparison of steady-state and transient flow conditions on reactive transport of contaminants in the vadose soil zone David Kuntz *, Peter Grathwohl Center for Applied Geosciences, University of Tübingen, Sigwartstraße 10, 72074 Tübingen, Germany

a r t i c l e

i n f o

Keywords: Transient flow Kinetic sorption Deterministic modeling Damköhler Groundwater risk assessment Vadose zone

s u m m a r y This paper investigates transport and fate of reactive compounds in the unsaturated subsoil (vadose zone) using numerical simulations of steady-state and transient flow scenarios. The objective is to investigate whether steady-state flow simulations accurately describe reactive contaminant transport under transient conditions in the field. The focus lies on organic compounds; advection, kinetic sorption and biodegradation are considered as relevant processes. In model scenarios the impact of steady-state and transient seepage water flow on solute concentrations at the groundwater table was evaluated. Sorption and degradation kinetics relative to the advective transport velocity were characterized by dimensionless Damköhler numbers for degradation (Dad) and sorption (Das). Compound properties representing two frequent pollutants (Lindane and Phenanthrene) and two types of sources were investigated: a constant inflow concentration (infinite source, with Phenanthrene) and a decreasing concentration due to desorption of the compound (finite source, with Lindane). Total mass degraded vs. cumulative flux to the groundwater (the lower model boundary), as well as the maximum breakthrough concentration at the outflow were evaluated. For each source type, combinations of sorption and degradation kinetics were identified, at which the transient simulation showed deviations from the predictions based on the steady-state simulation. Results are plotted as the difference in the observed maximum breakthrough concentration and the difference in cumulative mass fluxes as functions of degradation and sorption kinetics (‘‘Damköhler plots”). Overall results indicate that steady-state flow conditions are appropriate for most field scenarios; only extreme infiltration events lead to higher pollutant concentrations in transient simulations because of the short residence time of seepage water in the unsaturated zone and thus reduced biodegradation and sorption. Ó 2009 Elsevier B.V. All rights reserved.

Introduction For groundwater risk assessment the long-term fate of contaminants in the unsaturated zone is important and at a certain point of compliance their average breakthrough concentration is the main variable of interest (Leeson et al., 2002; Reichard et al., 1990). Risk-assessment tools are required for decision-making on long-term management of e.g. spills and diffuse pollution. For the necessary transport calculations of contaminants in the vadose soil zone, commonly assumptions of steady-state flow and constant water content are applied in model calculations (e.g. Finkel, 1998; Höhener et al., 2006; Maier and Grathwohl, 2005). Transient seepage water flow, which is a time dependent function of rainfall, surface runoff and evapotranspiration, is usually unknown. Most models use a fraction of the annual rainfall that * Corresponding author. Present address: Sprollstr. 9/1, 72108 Rottenburg a.N., Germany. Tel.: +49 74 83 9 28 99 0. E-mail addresses: [email protected] (D. Kuntz), [email protected] (P. Grathwohl). 0022-1694/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2009.02.006

matches the groundwater recharge as a constant seepage water infiltration rate (=mean annual recharge). The recharge is often obtained from water balances of catchments or hydrograph separation (for methods and their uncertainty see e.g. Scanlon et al., 2002, and the references cited therein). It is generally assumed that this recharge simplification is valid for long-term predictions of contaminant transport. Jury and Tseng (1999) compared onedimensional solute transport of conservative and degrading compounds between transient numerical models and steady-state analytical approaches and found the predicted breakthrough of a single solute pulse sufficiently reproduced by simple representations of the water flow regime, provided that accurate estimate of average net water flux and mean water content are available. Regarding slow sorption/desorption kinetics and varying source scenarios, no evaluation of transient modeling approaches has been reported yet. Under transient flow conditions residence times of seepage water in the soil matrix change (depending on soil properties and rainfall intensities), and sorption as well as biodegradation of a compound can be significantly influenced.

226

D. Kuntz, P. Grathwohl / Journal of Hydrology 369 (2009) 225–233

Although deterministic simulations of compound transport in the vadose zone under transient conditions are no exception, the impact of transient flow under possibly varying reaction conditions has been implemented mainly in stochastic models so far. Mantoglou and Gelhar (1987) developed a framework for modeling largescale transient unsaturated flow in spatially variable soil, which Mantoglou (1992) extended by considering boundary conditions and non-stationarity effects. Small and Mular (1987) examined flow and transport in one-dimensional homogeneous soils under stochastic infiltration. Russo et al. (1994a,b) studied the effects of varying mean soil water content and periodic water application rates on the fate of non-reactive solutes in unsaturated soils. Especially the leaching of pesticides and herbicides in agricultural soils has received considerable attention regarding transient flow, both by stochastic approaches (e.g. Jury and Roth, 1990) and laboratory column experiments (Pot et al., 2005). Foussereau et al. (2000) developed an analytical model that predicts solute transport in heterogeneous soils under transient flow applying a Lagrangian representation of solute transport. For regional or local groundwater risk assessment, deterministic or process-based models are applied more often to predict e.g. the maximum concentration at the groundwater table. The impact of transient flow on concentration and fluxes of contaminants at the groundwater table (compared to steady-state flow simulations) is usually neglected. However, especially when sorption equilibrium is not reached and biodegradation occurs in well permeable soils, concentration and fluxes at the groundwater table can be incorrectly predicted by models assuming steady-state flow conditions. The objective of this study is to examine the validity of the steady-state assumption of seepage water flow in deterministic modeling, and to delimit the conditions under which this simplification would lead to unreliable predictions. The transport behavior of two selected model compounds is compared between steady-state and transient flow conditions in model scenarios. These numerical experiments are employed in order to investigate extreme flow conditions in combination with a range of sorption/desorption and biodegradation rate constants. Such studies are essential for the upscaling of local models to the catchment and river basin scale (Jury and Tseng, 1999).

For the numerical solution of flow and transport equations the software code HYDRUS-1D version 3.0 was used (Šimùnek et al., 2005). Unsaturated water flow is based on the relationship developed by Richards (1931):

  @h @wm @ @w kðwm Þ m  kðwm Þ ¼ @wm @t @z @z

ð1Þ

with h being the soil’s volumetric water content (m3 m3), wm the matric potential head (m), k the unsaturated hydraulic conductivity depending on it (m s1), z the depth (m) and t the time (s). Modeling water flow in the vadose zone with Richards’ equation requires knowledge of the medium’s hydraulic properties, namely the specific water capacity, which is derived from the retention function h(wm), and the hydraulic conductivity function k(wm), which is calculated by applying Mualem’s (1976) integral Eq. (3) solved with the widely accepted water retention function of van Genuchten (1980) Eq. (2):

Se ðwm Þ ¼

1 þ jawm jn

1

kðSe Þ ¼ kf Ske  ½1  ð1  S1=m Þm 2 e

ð3Þ

with kf being the saturated hydraulic conductivity (m s1), while k (–) is an empirical parameter, set to 0.5, accounting for the water content dependency of pore connectivity and tortuosity (Mualem, 1976). Transport of a conservative tracer is calculated by the combined advection–dispersion equation for one-dimensional flow:

  @hcw @ @cw @qcw hDh ¼  @z @t @z @z

ð4Þ

with cw being the aqueous solute concentration (kg m3), Dh a hydrodynamic dispersion coefficient (m2 s1) that combines mechanical dispersion with molecular diffusion (see e.g. Crank, 1975), and q the vertical Darcy velocity (m s1), which is calculated by solving the Richards equation (Eq. 1). Biodegradation is implemented following first-order kinetics (or ‘‘half-life” kinetics), which is an accepted approximation when the substrate concentration is low compared to the available degrading biomass (Alexander, 1999):

rbio ¼ c  h  cw

ð5Þ 3

m

wm < 0 wm  0

ð2Þ

with a (m1), n (–) and m (–) being shape parameters (m is set to 1  1/n) and Se (–) the relative water saturation ranging from one

1

with rbio being the degradation rate in (kg m s ) and c a first-order degradation rate constant (s1). Mass transfer of organic compounds in sorption processes is considered to be limited by intraparticle diffusion. Diffusion equations, however, increase computing time significantly, but they can be approximated by a first-order approach (Grathwohl, 1998; van Genuchten and Wagenet, 1989):

@cs ¼ x  ðcs;eq  cs Þ @t

ð6Þ

Hereby cs is the solid phase concentration (kg kg1), x is a first-order rate constant (s1) and cs,eq is the solid phase concentration in equilibrium with the local aqueous phase concentration given by a linear sorption isotherm:

cs;eq ¼ K d  cw

ð7Þ 3

Transport modeling

(

to zero between the water content at full saturation (hS) and a residual water content (hR). The unsaturated conductivity is then given by:

1

Kd is the distribution coefficient (m kg ). A widely adopted procedure of predicting Kd for organic compounds presumes sorption of organic compounds is due to partitioning into organic matter and KOM usually is normalized to organic carbon (KOC, Karickhoff et al., 1979; Schwarzenbach and Westall, 1981; Karickhoff, 1984):

K d ¼ K OM  fOM ¼ K OC  fOC

ð8Þ

fOC is the fraction of organic carbon of the soil (–). KOC is commonly derived from tabulated water–octanol distribution coefficients (KOW) utilizing a linear free-energy relationship (Karickhoff, 1981; Chiou et al., 1983):

log K OC ¼ a  logðK OW þ bÞ

ð9Þ

a and b are empirical parameters, here set to 1 and 0.21, respectively (Karickhoff et al., 1979). Preferential flow could bypass degradation and sorption in the soil matrix and concentrations as well as fluxes to the subsoil may be higher than predicted without preferential flow. This is especially relevant for scenarios with a strong variance in infiltration. Under such conditions macropores may be activated at high soil moisture contents due to heavy precipitation events, which is significantly lower at steady-state flow calculations. Regarding transient calculations together with preferential flow the reader is referred to Vogel et al. (2000), Ray et al. (2004), Dusek et al. (2006) and Šimunek et al. (2003).

227

D. Kuntz, P. Grathwohl / Journal of Hydrology 369 (2009) 225–233

Setup of the model To evaluate the relevance of transient flow for contaminant transport in the vadose zone, we first defined a one-dimensional soil profile with a contaminant source. Two infiltration functions with the same annual water volume were defined, one as steadystate flow, the other generating regular infiltration pulses (Fig. 1). For the transient infiltration, we chose two infiltration events per year, each of 3 days duration (infiltration rate: 100 mm day1). These extreme infiltration events can be considered a ‘‘worst case” scenario for comparison of steady-state and transient flow conditions. Surface processes including precipitation, runoff, evapotranspiration, interception or rootwater uptake were not considered and a total net infiltration of 600 mm y1 was assumed in each scenario. Thus, the net infiltration equals the groundwater recharge. The model domain was set to a total length of 10 cm to maintain reasonable computing times. The top boundary condition was a Neumann boundary condition with either constant or transient inflow according to Fig. 1. The bottom boundary condition was set to ‘‘free drainage” (also called ‘‘seepage face”), representing an observation point above the groundwater table in a continuous soil profile. Initial condition for the steady-state simulation was the water content distribution developing under the same steady-state flow scenario. For the transient simulations the initial condition was set to a drained profile approximately 90 days after an infiltration event. Ranges of sorption/desorption and degradation kinetics were evaluated by changing the respective rate constants. This led to several hundred model realizations for each source type that covered all relevant combinations of different reaction kinetics. For each realization, we compared the steady-state and the extreme transient flow regime described above (Fig. 1). Damköhler numbers (see ‘‘Kinetic state of a model scenario – the Damköhler concept”) were used to evaluate and to quantify the differences between steady-state and transient model realizations as discussed in ‘‘Evaluating the relevance of transient flow”. Soil profile and source types Two source types were evaluated:  The finite source represents sorbed contaminants in the topsoil. The source is finite and contaminants are released by slow desorption, leading to a decreasing concentration with time. The mass available to either degradation or groundwater contamination Min (kg) at any time is the total mass desorbed in the source. Properties of the soil material in the source and

below were assumed identical (only one soil over the total model domain). Sorption and degradation kinetics were set the same in the source and in the soil below.  The infinite source assumes a constant contaminant concentration leaving the source. In the model this was represented by a uniform model domain (as in type one), together with a constant concentration inflow boundary condition (c0). The mass available to either degradation or groundwater contamination Min (kg) at any time is the total mass infiltrated. In both cases the soil below the source (finite source) or the total model domain (infinite source) was initially uncontaminated. The initial aqueous phase concentration was zero. Both source types and the concentration and mass nomenclature are summarized in Fig. 2. For the finite source, the model was run until the source contamination was completely removed. In the infinite source, starting at an initially clean soil profile, the contamination front will move downwards and is retarded by sorption and decreased in concentration by degradation. In steady-state flow eventually a stationary contamination profile is reached. When the steady-state distribution in the solid phase was not reached during the model time for slow sorption kinetics, the reduction of breakthrough concentration due to sorption is also very small. Thus the evaluation of breakthrough behavior or mass balances was not affected by the state of the solid phase contaminant distribution. The soil properties used in the simulations are listed in Table 1. Model compounds Biodegradation and sorption/desorption are defined by the degradation rate constants, the sorption distribution coefficient and the respective rate constant Eqs. (5)–(7). All three depend not only on compound properties but also on soil and flow characteristics. While KOC is well known for many compounds, their kinetic behavior regarding sorption and biodegradation is more difficult to predict in field scenarios. For the finite source a compound that exhibits relatively weak sorption was considered (e.g. Lindane). For the infinite source a stronger sorbing model compound was chosen (Phenanthrene). The physicochemical properties of these compounds and the sources thereof are listed in Table 1. Kinetic state of a model scenario – the Damköhler concept In the natural environment, soil characteristics are quite heterogeneous and differ from site to site. Physico-chemical parameters of organic pollutants cover a wide range, which influence sorption and biodegradation.

600 TR: 100 mm/d

80

400

Cumulative (Transient) Cumulative (Steady State) Infiltration (Transient) Infiltration (Steady State)

40

200

SS: 1,64 mm/d

0

0 0

60

120

180

240

300

Time [days] Fig. 1. Infiltration functions for steady-state and transient flow.

360

Cumulative infiltration [mm]

Infiltration rate [mm/day]

120

228

D. Kuntz, P. Grathwohl / Journal of Hydrology 369 (2009) 225–233 Table 1 Primary model parameters used in this study.

600mm/a (SS/TR)

Primary soil parameters

Infinite source c0 (const) / Min Finite source

28%

Mdeg / Msorb 72%

Model Domain

c0 / Min

cout / Mout

Hydraulic properties at 600 mm annual infiltration Steady-state Water content (h) Darcy velocity (q) Seepage velocity (v)

cout / Mout

The efficiency of reactions organic compounds undergo in a natural soil system depends on the reaction rate constants and the time available for the reaction. Under steady-state conditions the hydraulic soil properties and water infiltration rates define the residence time of seepage water in the soil. Sorption/desorption and degradation time scales are inverse proportional to the respective rate constants. The ratio of residence time to reaction time is expressed as Damköhler numbers (Damköhler, 1936; Jennings and Kirkner, 1984):

reaction rate transport timescale ¼ transport rate reaction timescale

ð10Þ

Calculation of the Damköhler number requires the definition of the timescales for advective transport, for biodegradation and for sorption. The residence time (tres) of seepage water in any soil profile is defined by the soil thickness (z) and the seepage velocity (v), which can be calculated by dividing the Darcy velocity (here: 600 mm per year in the steady-state scenario) by the soil’s water content:

t res ¼

z

v

z ¼ h q

ð11Þ

0.134 0.6 m y1 4.5 m y1

Model time (tmax) Domain depth (z)

34 years 10 cm

Compound parameters Model compound example Initial source concentration (cs,0) Inflow concentration (cw,0) Compound KOC used Sorption distribution coefficient (Kd)

Finite source Lindane 0.001 kg kg1 0 5000 l kg1a 10 l kg1

Infinite source Phenanthrene 0 0.46 mg l1 18000 l kg1b 36 l kg1

Reaction parameters

Residence time at steadystate flow (tres) Sorption direction Sorption equilibrium criterion (r) Sorption Damköhler range (Das) Sorption rate constant range (x) Degradation Damköhler range (Dad) Degradation rate constant range (c)

a

Finite source Source 2.28 d Desorption 105% of Kd

Soil 5.85 d Sorption 95% of Kd 0.01–100

1.23e9 to 1.23e5 s1

Infinite source Soil 8.13 d Sorption 95% of Kd 2.88e10 to 2.88e6 s1

0.01–100 1.42e8 to 1.42e4 s1

Calculated after Eq. (9), KOW from Xiao et al. (2004). Calculated after Eq. (9), KOW from Verschueren (1996).

cs ðtÞ ¼ cs;eq þ

Dc0 xtR e R

ð13Þ

cs,eq is the solid phase concentration at equilibrium:

ð12Þ

i denotes the order of the reaction, eliminating the concentration for i = 1. The residence time tres in Eq. (12) must represent the time available for the reaction. Here we follow the common assumption, that compounds are bioavailable to degradation only in the aqueous phase. Thus the bioavailable residence time even of solutes that are retarded by sorption equals the residence time of the seepage water, although they effectively travel slower than the percolating water. For sorption the reaction time depends on the distribution coefficient Kd and the rate constant x. Following Eqs. (6)–(9), we define the reaction timescale as the time required to achieve sorption equilibrium or near-equilibrium conditions. In a finite soil volume with no-flow the development of the concentration in the solid phase cs(t) after separation of the variables and integration of Eq. (6) is given by:

Transient (2  3 days) 0.291 36.5 m y1 125 m y1

Model setup parameters

b

For steady-state water flow under unit gradient conditions (no matric potential gradient) the water content is constant and can be calculated by Eq. (3), assuming the Darcy velocity equals the unsaturated hydraulic conductivity k(Se). Considering biodegradation the Damköhler number for first-order reactions is then defined as (Damköhler, 1936): i1 Dad ¼ t res  c  cw

42.7% 1.66 g cm3 9106 m s1 0.2% 0.427 0.065 0.0244 cm1 3.457 2.4 mm

Mdeg / Msorb

Fig. 2. Finite source (left) and infinite source (right). c0: aqueous inflow concentration, cout: aqueous concentration at the lower model boundary, Mi: mass fluxes input (in), degradation (deg), sorption (sorb) and outflow (out).

Da ¼

Porosity (e) Bulk density (rb) Saturated hydraulic conductivity (kf) Organic carbon content (fOC) Saturated water content (hS) Residual water content (hR) Van Genuchten parameter a Van Genuchten parameter n Dispersion length (al)

cs;eq ¼

K d  qb cs;0 þ cw;0 R h

ð14Þ

Dc0 denotes the initial concentration difference driving the process, defined as:

Dc0 ¼ cs;0  K d cw;0

ð15Þ

R is the retardation factor (R = 1 + qbKd/h), cs,0 the initial solid phase concentration (kg kg1), cw,0 the initial aqueous phase concentration (kg m3) and qb the bulk density (kg m3). The other terms (Kd, x, h) were defined above Eqs. (7), (6), and (1). The complementary development of the concentration in the aqueous phase is given by:

cw ðtÞ ¼ cw;eq 

Dc0 qb xtR  e R h

ð16Þ

D. Kuntz, P. Grathwohl / Journal of Hydrology 369 (2009) 225–233

cw,eq is the aqueous phase concentration at equilibrium (cs,eq/Kd). By defining the ratio of solid phase concentration to aqueous phase concentration r = cs(t)/cw(t) under near-equilibrium conditions close to Kd (r = 0.95 Kd for sorption and r = 1.05 Kd for desorption) we may solve for the reaction time t:

tðrÞ ¼ 

" # 1 r  cw;eq  cs;eq  ln Dc   0 xR 1 þ r qhb R

ð17Þ

This leads to the Damköhler number for sorption after Eq. (10):

Das ¼ t res =tðrÞ

ð18Þ

With the Damköhler numbers from Eqs. (12) and (18) we are able to characterize and compare soil profiles and contamination scenarios regarding their sorption and degradation kinetics independently from the primary parameters. We can evaluate the relevance of transient flow for transport in dependence of reaction kinetics as a function of these two Damköhler numbers, in order to provide results that are transferable to any scenario in which reaction times and residence times can be approximated. The validity of this concept has been shown recently by a sensitivity analysis that varied all primary parameters (Kuntz, 2008). Sorption behavior in a soil column was found to be similar over a wide range of primary parameters of the scenario, provided the Damköhler numbers were equal. Previous research has shown that for Damköhler numbers above 100, the reaction (e.g. sorption/desorption) can be regarded as an equilibrium process, while for Damköhler numbers below 0.01 kinetics are so slow, that they are almost negligible (Bold, 2004). Hence, the Damköhler numbers evaluated in this study range from 0.01 to 100, and for a selected set of primary parameters the degradation and sorption rate constants are calculated accordingly. Table 1 summarizes all primary parameters used for modeling as well as the resulting residence times at steady-state and the reaction rate constant range used to cover Damköhler numbers from 0.01 to 100.

229

conditions, while under transient flow and changing water saturation, the residence time of seepage water and thus the Damköhler numbers vary in space and time. In the text and figures of this paper, Damköhler numbers refer to the equivalent steady-state conditions with transient Damköhler numbers being lower during infiltration events and higher in no-flow periods. Evaluating the relevance of transient flow For the evaluation whether transient flow conditions have a significant influence, ‘‘concentration” and ‘‘fate” of the respective compound are considered as follows: (1) ‘‘Concentration”: the maximum breakthrough concentration observed at the lower model boundary, i.e. a point of compliance. This is the common variable of interest for groundwater risk assessment and used for action/no action decisions. (2) ‘‘Fate”: the total mass lost to degradation vs. the total massflow over the lower model boundary. A substance can be degraded, leached out at the lower model boundary or remain in the soil profile during the model period. Mass balances are independent of the breakthrough behavior and a better indicator of behavior and fate of contaminants. The solute concentration at the outflow of the model domain is monitored for each numerical parameter realization. It is expressed as the concentration relative to the potential maximum concentration (cw/cw,0). The potential maximum concentration cw,0 is the aqueous phase concentration in sorption equilibrium within the source (upper soil) at t = t0 for the finite source (FS), and the inflow concentration for the infinite source (IS):

cw cw  K d ðFSÞ ¼ cw;0 cs;source

cw cw ðISÞ ¼ cw;0 cw;inflow

ð19Þ

The difference between the observed maximum concentration (dc/ c0) is calculated by:

cw;max cw;max cw;max ¼ ðTRÞ  ðSSÞ cw;0 cw;0 cw;0

Damköhler numbers in steady-state and transient flow

d

For an annual infiltration rate of 600 mm/y and soil hydraulic properties as introduced in Table 1, the residence time of the percolating water under steady-state flow calculates to 195 h for a 10 cm column at a constant water content of 13.4%. For the infinite source type, to evaluate sorption Damköhler ranges from 0.01 to 100, the sorption rate constant was varied according to Eqs. (17) and (18) from 2.88e10 s1 to 2.88e6 s1. For the finite source type the sorption rate constant calculates different due to the different equilibrium sorption distribution and the difference in sorption direction between source and soil below (rate constant varies from 1.23e9 s1 to 1.23e5 s1). In the corresponding transient scenarios the rate constants were kept at the same values, but the residence time of the percolating water changes greatly during infiltration events and the no-flow periods in between. The actual Damköhler numbers for the transient scenarios thus differ significantly from the steady-state scenarios. Assuming a constant water content in the column is reached during an infiltration event, the infiltration intensity of 100 mm/d (=36.5 m/y) leads to residence times in the column of approximately 7 h at a water content of 29.1% (in the no-flow periods between infiltration events, the residence time and the Damköhler number approach infinity). Considering this, an infinite source type soil system where sorption is close to equilibrium (Damköhler number of 100) at steady-state flow, is shifted towards non-equilibrium conditions during the infiltration events with sorption Damköhler numbers of approximately 3.6. We conclude that meaningful Damköhler numbers can be defined only for homogeneous soils under steady-state flow

Hereby ‘‘TR” and ‘‘SS” are not variables but indicate the transient and the steady-state scenario, respectively. Overall degradation vs. transport to the groundwater (‘‘Fate”) is quantified by the amount degraded (Mdeg, [kg]) and the amount leaving the lower model boundary (Mout, [kg]). To keep different reaction rates and both source types comparable, we only consider the mass available (Min, [kg]) to degradation or transport, which is the desorbed mass in the finite source and the infiltrated (and thus total) mass for the infinite source. For any model realization, Mdeg and Mout are competitive and can thus be compared by a simple relative difference (DF):

DF SS=TR ¼

M out  Mdeg M in

ð20Þ

ð21Þ

Assuming the dissolved mass and the mass sorbed below the source are very small compared to Min at the end of the model time, DF ranges from 1 (total degradation, no groundwater contamination) to +1 (no degradation, all mass passed the lower model boundary). A value of zero for DF indicates an equal amount of mass lost to degradation and to the groundwater. To compare steady-state vs. transient transport flow conditions, DF is calculated for each flow regime (DFSS and DFTR) and the difference between both is quantified, here termed as evaluation function for ‘‘fate” (Ef):

Ef ¼ DF SS  DF TR

ð22Þ

Negative values of Ef, denote lower overall degradation at transient flow compared to steady-state flow. For positive values of Ef, tran-

230

D. Kuntz, P. Grathwohl / Journal of Hydrology 369 (2009) 225–233

sient flow ‘‘enhances” degradation. For Ef close to zero, no significant difference exists between steady-state and transient flow transport regarding compound fate. The theoretical range of Ef is from 2 (total degradation in SS with no degradation in TR) to +2 (vice versa). From a groundwater risk assessment point of view, negative Ef values indicate that groundwater contamination is underestimated by steady-state flow models.

tions for which the latter must be viewed as potentially sensitive to transient flow. Finite source type The finite source type exhibits for fast sorption kinetics a typical breakthrough curve, with retardation depending on the compound’s Kd value. While for steady-state the breakthrough curve is smooth, it occurs in pulses during transient flow. For slow sorption/desorption, less mass desorbs in the source, but sorption in deeper parts of the profile is also slower. At very slow sorption/ desorption, the percolating water takes up only trace amounts at its passage through the source and consequently only very little is sorbed below. This leads to a slowly declining solute concentration until the source is depleted. Concerning the concentration at the point of compliance (cout, Fig. 3a) we find no difference between steady-state and transient flow (dc/c0 = 0) for degradation Damköhler numbers >1 due to total degradation of the solute. Towards slow or no degradation, differences are more pronounced (0.3 < dc/c0 < 0.2) depending on the sorption kinetics. Towards fast sorption kinetics (Das > 10) the maximum observed breakthrough concentration is greater in the steady-state scenario, due to the additional dispersion in the transient flow

Results – Damköhler plots For each source type 100 model realizations were calculated in steady-state and transient flow mode. The 100 scenarios are equally distributed in log space over a Damköhler range from 0.01 to 100 in both sorption (Das) and degradation kinetics (Dad) at steady-state flow (SS). Eqs. (20) and (22) are plotted in Fig. 3 as functions of both Damköhler numbers for each source type. The isoline-values refer to the difference between the steady-state flow scenario and the transient reference scenario of this work (cp. ‘‘Setup of the model”), and are therefore only a qualitative measure of the impact of transient flow and not directly transferable to individual field sites or model scenarios. These Damköhler plots allow fast identification of combinations of sorption and degradation kinetics for which steady-state models provide accurate predictions of compound transport, and combina-

Infinite Source

Finite Source Concentration c/c0

2

Concentration c/c0

2

10

10

a

b no difference

1

Dadegradation (SS)

10

1

10 0

0

10

1 0.

0.5

5 0.3 0.2

-

0.125

-1

> 0.15

-2

10

-0. 2 < -0.25

0.1

10

-2

-1

10

ferenc no dif

10

0

1

10

2

10

c 1

-2

10

-1

10

dif fer

10

0.75

1

2

10

10

d en

1

10

ce

0 0

0

10

Fate Ef

2

10

no

10

10

10

Fate Ef

2

10

0.2 5 0.5

< -1.2

0

10

-0.9

6

-0.

> 1.0

3 -0.

-1

-1

10

10

n -2

10

e

-2

-1

10

Dadegradation (SS)

> 0.65

0

10

iff od

ere

nce

-2

-2

10

-1

10

0

10

Dasorption (SS)

1

10

2

10

10

-2

10

-1

10

0

10

1

10

2

10

Dasorption (SS)

Fig. 3. Damköhler plots for the finite (left) and infinite (right) source, plotting differences between steady-state and transient flow regarding breakthrough concentration Eq. (20) and compound mass balances Eq. (22). Damköhler numbers specify the steady-state scenario only.

231

D. Kuntz, P. Grathwohl / Journal of Hydrology 369 (2009) 225–233

mode. Towards equilibrium sorption this difference vanishes, although zero difference was not covered by the evaluation (sorption kinetics is shifted towards non-equilibrium during transient infiltration events). Even at sorption Damköhler values of 100 the maximum breakthrough concentration in steady-state flow is approximately 18%-points higher than in the transient scenario (cp. Fig. 4a). The maximum difference was found around sorption Damköhler numbers of 11, this point is however also dependent on the transient infiltration function (Fig. 1) and may not be transferred directly to other transient scenarios. For very slow sorption/desorption (Das < 0.1), the source leaches slowly and the solute concentration in the water is low (c/c0 < 15%) in steady-state flow scenarios. In the transient scenario the average concentration is even lower due to compound dispersion by infiltration pulses and degradation in no-flow periods. Peak concentrations however reach a breakthrough concentration of up to 25% of the potential maximum, thus exceeding the concentration of the steadystate calculation at the early events by 10–15%-points (cp. Fig. 4b). Summarizing the evaluation of breakthrough concentration for the finite source we observe that transient flow affects the breakthrough curve similar to additional dispersion, thus stretching the curve in time and lowering the peak concentration. Strong infiltration events travel faster through the soil profile and transport the dissolved compound faster into deeper regions of the soil compared to steady-state flow. When considering peak concentrations, the highest concentration will be observed in a steady-state flow scenario for fast sorption kinetics, while for slow sorption kinetics pulsed breakthrough events in transient scenarios may exceed the equivalent steady-state flow concentration. When we observe only average concentrations and neglect single peak events in transient flow, we conclude that steady-state flow simulations may generally overestimate the breakthrough concentration. Concerning compound fate, for slow or no degradation eventually all mass reaches the groundwater (Mout = Min, Mdeg = 0) and for fast degradation all mass is degraded (Mout = 0, Mdeg = Min) independent of advection (Fig. 3c). At fast or equilibrium sorption kinetics (sorption Damköhler numbers >10), the total time for contaminant leaching is equal for both flow conditions, all mass is desorbed quickly and thus bioavailable. Therefore the final mass balance (cp. Fig. 4c) shows no difference between steady-state and transient flow. At Damköhler numbers <1 for sorption we observe that transient flow ‘‘enhances” degradation, with a maximum around Damköhler numbers of 0.20.4 for degradation. Specifically: at Das = 0.03 and Dad = 0.6 we observe 40% of the total desorbed mass being degraded under steadystate flow while this increases to 89% under transient flow (Fig. 4d). After all mass leached from the source (Mout + Mdeg = Min) Eq. (21) yields DFSS = 0.2 and DFTR = 0.78, and thus Ef = 0.98. Note that the enhancement of degradation under transient flow is due to the extreme transient infiltration conditions. At steadystate flow, most of the mass that desorbs in the source is continuously transported across the lower model boundary at slow sorption kinetics. During transient infiltration events, the seepage water flow rate is higher than in steady-state flow, but due to slow desorption only little mass leaches per event. With only two infiltration events in the transient flow scenario, the time until all mass is desorbed is longer compared to the steady-state flow mode. This in turn allows for more degradation in the source between infiltration events. We expect differences between transient and steadystate flow to be greatly reduced under natural conditions, where infiltration intensities are lower and intervals shorter. Infinite source type For the infinite source scenario a concentration front migrates into an initially clean soil profile. Under steady-state flow, eventu-

50%

(a) Finite source - Breakthrough (Dad =0.1 Das=100)

40% 30%

Concentration plots cout/c0 (SS)

20%

cout/c0 (TR)

10% 0% 30%

(b) Finite source - Breakthrough (Dad =0.03 Das=0.05) 20%

10%

0% 100% 80%

Mass balance plots Mdeg (SS)

60%

(c) Finite source Mass Balance (Dad=0.1 Das=50)

Mout (SS)

40%

Mdeg (TR) Mout (TR)

20% 0% 100%

(d) Finite source - Mass Balance (Dad =0.6 Das=0.03)

80% 60% 40% 20% 0% 100%

(e) Infinite source - Breakthrough (Dad =4.6 Das=0.2)

80% 60% 40% 20% 0% 100%

(f) Infinite source - Breakthrough (Dad =4.6 Das=25)

80% 60% 40% 20% 0% 100%

(g) Infinite source - Mass Balance (Dad =4.6 Das=0.2)

80% 60% 40% 20% 0% 100%

(h) Infinite source - Mass Balance (Dad =4.6 Das=25)

80% 60% 40% 20% 0%

0

5

10

15

20

25

30

35

Years Fig. 4. Breakthrough concentration and mass balance for selected scenarios (steady-state, SS: dashed lines; transient, TR: solid lines). Damköhler numbers specify the steady-state scenario only.

ally a stationary concentration profile over depth develops. Its shape only depends on degradation kinetics. For fast sorption the contamination front advances equally in the aqueous and solid phase, at slow sorption kinetics we observe early breakthrough while the solid phase takes up the contaminant more slowly until

232

D. Kuntz, P. Grathwohl / Journal of Hydrology 369 (2009) 225–233

sorption equilibrium is reached. In this source type transient flow only yields higher breakthrough concentrations compared to steady-state conditions. Under steady-state flow the breakthrough curve is smooth and reaches a constant concentration plateau. Strong infiltration events in transient flow can bypass degradation due to shorter residence times of the seepage water in the soil profile. For very fast (Dad > 50) and very slow (Dad < 0.1) degradation kinetics, no difference exists in the maximum breakthrough concentration between steady-state and transient flow due to either total or no degradation (Fig. 3b). The largest difference in the observed maximum concentration between both flow conditions can be expected when degradation is fast enough to degrade the total contaminant mass just before the seepage water reaches the lower model boundary. This is the case at degradation Damköhler numbers around 4–5 (Fig. 3b). For steady-state flow we observe no breakthrough, while transient flow generates breakthrough pulses. The concentration of these breakthrough pulses depends on the strength of the infiltration events and on sorption. For no or slow sorption and strong infiltration events a high breakthrough concentration is reached. At Das = 0.2 and Dad = 4.6 under steady-state flow the concentration at the lower model boundary (cout) reaches about 3% of the inflow concentration (c0), while at transient flow single events break through with about 80% of the inflow concentration (Fig. 4e). At faster sorption kinetics the breakthrough pulses are buffered by sorption in the soil profile. Towards fast sorption the breakthrough curve under transient flow approaches the steady-state flow breakthrough curve and the difference in the maximum breakthrough concentration diminishes (Fig. 4f) and vanishes at true equilibrium sorption. Contaminant fate for the infinite source scenario is directly coupled to breakthrough concentration (Fig. 4g and h). If under steadystate no breakthrough occurs (all mass is degraded: Mdeg = Min), transient conditions lead to breakthrough pulses shifting the mass balance towards more groundwater contamination (Mout). This is reflected in the Damköhler plot (Fig. 3d) which is similar to the concentration information.

pulses occurring under transient conditions are not covered. For compounds that are degraded between Damköhler numbers of 0.5 and 50, steady-state flow results in long residence times in the soil profile, allowing for degradation of most of the contaminant mass, while extreme infiltration events can bypass degradation. This effect is most pronounced for weak or slow sorbing compounds. At sorption Damköhler numbers >1.0 transient pulses are effectively buffered. While we observe differences in mass balance and the breakthrough concentration of up to 80%, the given values here represent one extreme transient scenario, under more natural infiltration conditions the difference will be less pronounced. In conclusion, the assumption of steady-state flow conditions for the simulation of reactive transport of pollutants in the unsaturated zone is appropriate for a wide range of sorption/desorption and degradation rate constants (‘‘no difference” regions in Fig. 3). When threshold concentrations at a specific low level must not be exceeded even for short time periods, then transient simulations may be needed, but only if degradation or sorption are very slow. As formerly mentioned, the simulations and the resulting Damköhler plots provided here are based on an artificial transient scenario. The relevance of transient flow under natural conditions is investigated in Kuntz (2008) for several precipitation time series at various European locations, considering also evapotranspiration and root water uptake for realistic infiltration rates. It could be shown, that even for extreme precipitation locations, considering both total annual precipitation as well as variation in precipitation, the differences between transient simulations and steady-state simulation remain low (<15% of the total mass/potential maximum concentration) when considering long-term fate and average breakthrough concentrations. Only at specific cases resembling the infinite source type, steady-state simulations may severely underestimate breakthrough concentrations at single infiltration events.

Conclusions

This work was supported by the European Integrated Project AquaTerra (GOCE 505428) under the Community’s Sixth Framework Programme.

In scenarios where a topsoil contamination releases contaminants by slow desorption (finite source) transient seepage water flow has two significant effects on the breakthrough history over the total removal: single extreme infiltration events lead to higher breakthrough concentration peaks than predicted by steady-state flow models, and the compound spreads faster over the soil profile. This effect is more pronounced at slow sorption kinetics, leads to longer time for total contaminant removal from the source and decreases the maximum average concentration. No or little differences in mass balances between steady-state and transient models can be expected under natural infiltration conditions. Considering the breakthrough concentration, a maximum difference of about 20–30% (of the potential maximum concentration cw,0) between transient and steady-state flow was observed for degradation kinetics with Damköhler numbers <0.1. For sorption kinetics with Damköhler numbers <0.1 a transient simulation may show higher peak concentrations compared to a steady-state simulation while for sorption Damköhler numbers >1 steady-state flow calculations can overestimate breakthrough concentrations even compared to transient infiltration peaks. The difference, however, will decrease further under natural infiltration conditions due to the extreme character of the input function generating transient flow in this work. For infinite source scenarios steady-state flow conditions may underestimate groundwater contamination because concentration

Acknowledgement

References Alexander, M., 1999. Biodegradation and Bioremediation, second ed. Academic Press, Ithaca, New York. ISBN: 0-12-049861-8. Bold, S., 2004. Process-based Prediction of the Long-term Risk of Groundwater Pollution by Organic Non-volatile Contaminants. Dissertation, University of Tübingen TGA C72, Tübingen. Chiou, C.T., Porter, P.E., Schmedding, D.W., 1983. Partition equilibria of nonionic organic compounds between soil organic matter and water. Environmental Science & Technology 17 (4), 227–231. Crank, J., 1975. The Mathematics of Diffusion, second ed. Oxford University Press, New York. Damköhler, G., 1936. Einflüsse der Strömung, Diffusion und des Wärmeüberganges auf die Leistung von Reaktionsöfen. Zeitschrift Elektrochemie 42 (12), 846–862. Dusek, J., Vogel, T., Lichner, L., Cipakova, A., Dohnal, M., 2006. Simulated cadmium transport in macroporous soil during heavy rainstorm using dual-permeability approach. Biologia 61, 251–254. Finkel, M., 1998. Quantitative Beschreibung des Transports von Polyzyklischen Aromatischen Kohlenwasserstoffen (PAK) und Tensiden in porösen Medien. Dissertation, University of Tübingen, TGA C47, Tübingen. Foussereau, X., Graham, W.D., Ashie Akpoji, G., Destouni, G., Rao, P.S.C., 2000. Stochastic analysis of transport in unsaturated heterogeneous soils under transient flow regimes. Water Resources Research 36 (4), 911–921. Grathwohl, P., 1998. Diffusion in Natural Porous Media. Kluwer Academic Publisher, Boston. ISBN: 0-7923-8102-5. Höhener, P., Dakhel, N., Christophersen, M., Broholm, M., Kjeldsen, P., 2006. Biodegradation of hydrocarbons vapors: comparison of laboratory studies and field investigations in the vadose zone at the emplaced fuel source experiment, Airbase Vrløse, Denmark. Journal of Contaminant Hydrology 88 (3–4), 337– 358.

D. Kuntz, P. Grathwohl / Journal of Hydrology 369 (2009) 225–233 Jennings, A.A., Kirkner, D.J., 1984. Instantaneous equilibrium approximation analysis. Journal of Hydraulic Engineering 110 (12), 1700–1717. Jury, W.A., Roth, R., 1990. Transfer Function and Solute Movement through Soil: Theory and Applications. Birkhäuser, Basel, Switzerland. Jury, W.A., Tseng, P.-H., 1999. Options for modeling ground water pollution potential by dissolved chemicals. Geophysical Monograph 108, 21–29. Karickhoff, S.W., 1984. Organic pollutant sorption in aquatic systems. Journal of Hydraulic Engineering 110 (6), 707–720. Karickhoff, S.W., 1981. Semi-empirical estimation of sorption of hydrophobic pollutants on natural sediments and soils. Chemosphere 10 (8), 833– 846. Karickhoff, S.W., Brown, D.S., Scott, T.A., 1979. Sorption of hydrophobic pollutants on natural sediments. Water Research 13 (3), 241–248. Kuntz, D., 2008. The Relevance of Transient Flow on Reactive Contaminant Transport in the Unsaturated Soil Zone. Dissertation, University of Tübingen, TGA Reihe C: Hydro-, Ingenieur- und Umwelt-geologie, Institute of Geosciences. ISSN: 0935-4948 (Print), ISSN: 1610-4706 (Internet). Leeson, J., Edwards, A., Smith, J.W.N., Potter, H.A.B., Bourn, M., 2002. Hydrological Risk Assessments for Landfills and the Derivation of Groundwater Control and Trigger Levels. Environmental Agency, Solihull, U.K. Maier, U., Grathwohl, P., 2005. Natural attenuation in the unsaturated zone and shallow groundwater: coupled modeling of vapor phase diffusion, biogeochemical processes and transport across the capillary fringe. In: Nützmann, G., Viotti, P., Aagard, P. (Eds.), Reactive Transport in Soil and Groundwater. Springer, Berlin, pp. 141–155. Mantoglou, A., 1992. A theoretical approach for modeling unsaturated flow in spatially variable soils: effective flow models in finite domains and nonstationarity. Water Resources Research 28 (1), 251–267. Mantoglou, A., Gelhar, L.W., 1987. Stochastic modeling of large scale transient unsaturated flow systems. Water Resources Research 23 (1), 37–46. Mualem, Y., 1976. New model for predicting the hydraulic conductivity of unsaturated porous media. Water Resources Research 12 (3), 513–522. Pot, V., Šimunek, J., Benoit, P., Cosquet, A., Yra, A., Martínez-Cordón, M.-J., 2005. Impact of rainfall intensity on the transport of two herbicides in undisturbed grassed filter strip soil cores. Journal of Contaminant Hydrology 81, 63–88. Ray, C., Vogel, T., Dusek, J., 2004. Modeling depth-variant and domain-specific sorption and biodegradation in dual-permeability media. Journal of Contaminant Hydrology 70 (1–2), 63–87. Reichard, E., Cranor, C., Raucher, R., Zaponi, G., 1990. Groundwater Contamination Risk Assessment – A Guide to Understanding and Managing Uncertainties.

233

International Association of Hydrological Sciences (IAHS) publication 196, Wallingford, ISBN: 0-947571-72-8. Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. Physics 1, 318–333. Russo, D., Zaidel, J., Laufer, A., 1994a. Stochastic analysis of solute transport in partially saturated heterogeneous soil, 1, numerical experiments. Water Resources Research 30, 769–779. Russo, D., Zaidel, J., Laufer, A., 1994b. Stochastic analysis of solute transport in partially saturated heterogeneous soil, 2, prediction of solute spreading and breakthrough. Water Resources Research 30, 781–790. Scanlon, B.R., Healy, R.W., Cook, P.G., 2002. Choosing appropriate techniques for quantifying groundwater recharge. Hydrogeology Journal 10 (1), 18–39. Schwarzenbach, R.P., Westall, J., 1981. Transport of nonpolar organic compounds from surface water of groundwater. Laboratory sorption studies. Environmental Science & Technology 15 (11), 1360–1367. Šimunek, J., Jarvis, N.J., van Genuchten, M.Th., Gärdenäs, A., 2003. Nonequilibrium and preferential flow and transport in the vadose zone: review and case study. Journal of Hydrology 272, 14–35. Šimùnek, J., van Genuchten, M.Th., Šejna, M., 2005. The HYDRUS-1D Software Package for Simulating the One-dimensional Movement of Water, Heat, and Multiple Solutes in Variably-saturated Media, Version 3.0. HYDRUS Software Series 1, Department of Environmental Science, University of California, Riverside, CA. Small, M.J., Mular, J.R., 1987. Long-term pollutant degradation in the unsaturated zone with stochastic rainfall infiltration. Water Resources Research 23, 2246– 2256. van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44 (5), 892–898. van Genuchten, M.Th., Wagenet, R.J., 1989. Two-site/two-region models for pesticide transport and degradation: theoretical development and analytical solutions. Soil Science Society of America Journal 53, 1303–1310. Verschueren, K., 1996. Handbook of Environmental Data on Organic Chemicals, third ed. John Wiley & Sons, Canada. ISBN: 0-471-28659-1. Vogel, T., Gerke, H.H., Zhang, R., Van Genuchten, M.T., 2000. Modeling flow and transport in a two-dimensional dual-permeability system with spatially variable hydraulic properties. Journal of Hydrology 238 (1–2), 78–89. Xiao, H., Li, N., Wania, F., 2004. Compilation, evaluation, and selection of physical– chemical property data for alpha-, beta-, and gamma-hexachlorocyclohexane. Journal of Chemical and Engineering Data 49 (2), 173–185.