Journal of Contaminant Hydrology 173 (2015) 69–82
Contents lists available at ScienceDirect
Journal of Contaminant Hydrology journal homepage: www.elsevier.com/locate/jconhyd
Generation of dense plume fingers in saturated–unsaturated homogeneous porous media Clemens J.M. Cremer ⁎, Thomas Graf Institute of Fluid Mechanics and Environmental Physics in Civil Engineering, Leibniz Universität Hannover, Appelstr. 9A, 30167 Hannover, Germany
a r t i c l e
i n f o
Article history: Received 7 July 2014 Received in revised form 10 November 2014 Accepted 26 November 2014 Available online 9 December 2014 Keywords: Model Groundwater Saturated Unsaturated Density Numerical trigger Fingering
a b s t r a c t Flow under variable-density conditions is widespread, occurring in geothermal reservoirs, at waste disposal sites or due to saltwater intrusion. The migration of dense plumes typically results in the formation of vertical plume fingers which are known to be triggered by material heterogeneity or by variations in source concentration that causes the density variation. Using a numerical groundwater model, six perturbation methods are tested under saturated and unsaturated flow conditions to mimic heterogeneity and concentration variations on the pore scale in order to realistically generate dense fingers. A laboratory-scale sand tank experiment is numerically simulated, and the perturbation methods are evaluated by comparing plume fingers obtained from the laboratory experiment with numerically simulated fingers. Dense plume fingering for saturated flow can best be reproduced with a spatially random, time-constant perturbation of the solute source. For unsaturated flow, a spatially and temporally random noise of solute concentration or a random conductivity field adequately simulate plume fingering. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Leachate that invades soil, e.g. from waste disposal sites, frequently possesses higher solute concentration and therefore higher fluid density than ambient freshwater. Spatial variation of water density can lead to variable-density flow, which is a key process in the transport of contaminants, salt transport in coastal aquifers or hot geothermal plumes (Diersch and Kolditz, 2002; Illangasekare et al., 2006; Simmons, 2005; Simmons et al., 2001; Stöckl and Houben, 2012; Van Dam et al., 2009). The migration of dense plumes through a porous medium typically results in the formation of lobe-shaped instabilities or vertical plume fingers, which are known to have a nonnegligible effect on the propagation of the plume (Xie et al., 2011). Flow direction in solute plume fingers is downwards, counterbalanced by an upward flow of less dense fluid between ⁎ Corresponding author. E-mail address: [email protected] (C. Cremer). URL: http://www.hydromech.uni-hannover.de (C. Cremer).
http://dx.doi.org/10.1016/j.jconhyd.2014.11.008 0169-7722/© 2014 Elsevier B.V. All rights reserved.
the fingers. The generation of plume fingers can be attributed to material heterogeneity, which is known to trigger the formation of dense fingers (Goswami et al., 2011; Schincariol, 1998; Schincariol and Schwartz, 1990; Schincariol et al., 1997). In homogeneous media, however, fingers are also created even though the material is uniform with identical sizes of material grains (e.g. sand grains) (Goswami et al., 2011; Simmons et al., 2002; Van Dam et al., 2009; Wooding et al., 1997b). The reason can be assumed to be that pore-scale heterogeneity leading to different flow velocities also exists in homogeneous media due to two effects: (i) random variations of solute concentration leading to varying buoyancy effects, which results in different average flow velocities, and (ii) grains of identical size which may randomly arrange differently, thus forming tetrahedrons, hexahedrons or octahedrons, where each arrangement creates pores of varying diameters, thus also resulting in different pore-scale velocities. Usually, variable-density flow models do not consider effects (i) and (ii), such that the numerical model does not apply any external perturbation to trigger plume fingering. In that case, numerical errors of the computing process are expected to
70
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
perturb the flow field and to create unstable plume fingers. The problem with this approach, however, is that “these instabilities are not physically realistic and are essentially uncontrollable because their character depends on the extent to which numerical errors develop” (Schincariol et al., 1994). This has later been confirmed by Simmons et al. (1999) who noted that numerical errors “may occur by themselves and are uncontrollable as location and magnitude of the error may vary with time. Small changes to dispersion parameters, time step sizes, or spatial discretization throughout a series of simulations can cause different instabilities to form.” Xie et al. (2011) added that numerical errors are not easy to quantify. As an alternative to relying on internal numerical errors, deliberately triggering instabilities by external perturbations represents a controlled way to account for pore-scale heterogeneity. Perturbations can be introduced through a noise applied to the solute concentration or through a random conductivity field. List (1965) induced perturbations of the source concentration created by superposition of spatially oscillatory terms with sinusoidal form and performed a theoretical analysis on the relationship between the nature of the perturbation and the resulting instabilities. List (1965) identified the following parameters responsible for instability growth: the transverse Rayleigh number, the ratio of the transverse to the longitudinal Rayleigh numbers, and the nondimensional wave number. However, while List's stability theory provides a useful model for determining stability in a variable-density system, Schincariol (1998) criticized that the theory cannot account for dispersive dissipation of the dense plume as a function of travel distance and is therefore less useful in predicting the growth rate of instabilities in real systems. Schincariol and Schwartz (1990) conducted physical experiments on the laboratory-scale in homogeneous, layered, and lenticular porous media in order to further investigate the influence of hydraulic conductivity and density differences on the formation of instabilities. It was observed that small density differences around 0.8 mg cm−3 between a higher density water overlying a less dense water results in the formation of gravitational instabilities. Schincariol et al. (1994) numerically reproduced the physical experiment by Schincariol and Schwartz (1990). Numerically generated errors and their influence on the development of instabilities were examined using the Courant- and Peclet-stability criteria. Schincariol et al. (1994) concluded that the formation of unstable flow depends on the Rayleigh-number and on the wavelength of the perturbation. Schincariol et al. (1997) and Schincariol (1998) tested the applicability of the stability theory developed by List (1965) to variable-density flow in homogeneous and heterogeneous media. Numerical simulations were carried out, where perturbations were induced by lab-scale heterogeneity of the conductivity field. Schincariol et al. (1997) concluded that besides the Raleigh number and wavelength, statistical characteristics of the conductivity field (mean, variance, correlation length scales) play an important role in the development of instabilities. Simmons et al. (1999) numerically resimulated an experiment of an evaporating salt lake developed by Wooding et al. (1997a,b). A spatially and temporally random noise of the salt concentration in the lake was used to apply perturbation and to mimic field-scale heterogeneity. Using that perturbation technique, Simmons et al. (1999) could realistically reproduce
experimental results by Wooding et al. (1997a,b). Xie et al. (2011) adopted this perturbation mechanism to study the effects of dynamic solute boundary conditions on the migration of dense plumes using the modified solute Elder problem (Voss and Souza, 1987). That problem was later reused and slightly modified in multiple realizations to analyze the speed of free convective fingering (Xie et al., 2012). Weatherill et al. (2004) investigated the onset of convective roll instabilities due to thermal convection in a fluid-saturated porous layer that is uniformly being heated from below. The Horton–Rogers–Lapwood (HRL) problem (Horton and Rogers, 1945; Lapwood, 1948) was numerically resimulated with an initial perturbation at the domain center to initiate convection. A detailed investigation of that perturbation method is absent, and it therefore remains unknown whether the trigger used by Weatherill et al. (2004) really affects the onset of convection. Recently, Illangasekare et al. (2006) and Van Dam et al. (2009) obtained evidence of lobe-shaped instabilities from field measurements. This highlights the importance of unstable density-driven flow in natural field settings and furthermore the importance of an adequate representation of these phenomena in numerical models. The objective of the present study is to evaluate six different perturbation mechanisms in order to incorporate pore-scale heterogeneity into a numerical model such that dense fingers are realistically being generated under saturated and unsaturated flow conditions. While partly relying on the evaluation of well known perturbation mechanisms, this study also introduces several mechanisms that are new to variable-density flow simulations. To the authors' knowledge, (i) a spatially and temporally random noise of solute concentration, (ii) the application of a spatially random, time-constant perturbation of the solute source, and (iii) a random conductivity field have previously not been examined as perturbation mechanisms in the literature of variable-density flow. Success of all perturbation methods is validated by comparison of numerical results with a physical laboratory model conducted by (Simmons et al., 2002). 2. Physical model Simmons et al. (2002) conducted laboratory experiments of variable-density flow and salt transport in a 2D vertical sand tank. Ambient horizontal flow was not considered. The focus was put on variable-density free convective flow including the generation of dense salt fingers. Stained calcium chloride solutions of different densities were injected along the central part of the upper boundary. The experiments were carried out under fully saturated as well as partly saturated conditions. For the saturated case, injection of dense solutes created ponding along the source, such that the water table in the source was a few millimeters higher than those along the left and right top boundaries. The unsaturated experiment was carried out by draining the water level to 44 cm below the top surface of the sand. The present study focuses on the saturated high density (SHD) case and the unsaturated high density (UHD) case (Simmons et al., 2002) where the density of the source concentration was 1235 kg m−3. Fig. 1A and B presents the plume development of the physical SHD and UHD experiments after 40 min of dense
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
plume fingering. The density difference between freshwater initially filling the tank and the dense calcium chloride solution leads to downward density-driven migration within distinct plume fingers. Downward movement within these salt fingers is counterbalanced by upward movement of freshwater between salt fingers. In the UHD case, fingering did not take place in the unsaturated zone because of the small unsaturated hydraulic conductivity. Note that, in the UHD case, fingering started at the transition between the unsaturated and the saturated zones. Fig. 1C and D also show the number of fingers versus time for the SHD and UHD cases (Simmons et al., 2002). The number of plume fingers for the SHD consistently decreases with time. However, the UHD case shows an initial increasing and afterwards a decreasing number of fingers. According to Simmons et al. (2002), the time when plume fingers become more numerous coincides with the dense plume reaching the water table. While flow is dominated by capillarity in the unsaturated zone, flow in the saturated zone is dominated by density differences. The decrease in number of fingers can be attributed to the fact that multiple small fingers merge into smaller fingers over time (Fig. 1D). The evolution of the number of fingers with time will be used here to test different perturbation mechanisms. More
71
details on and an interpretation of the physical model can be found in Simmons et al. (2002). 3. Numerical model The numerical variably saturated variable-density flow and solute transport model HydroGeoSphere (Therrien et al., 2012) is used here to re-simulate the SHD and UHD laboratory experiments conducted by Simmons et al. (2002). Transient flow through a variably saturated porous medium is described by the Richards' equation (Therrien et al., 2012) ∂ðθS Sw Þ ¼ −∇q ∂t
ð1Þ
where θs [−] is saturated water content and Sw [−] is water saturation. The specific discharge q [L T−1] is given by the Darcy equation ρ−ρ0 ∇z q ¼ −K 0 kr ∇ h0 ρ0
ð2Þ
where K0 [L T−1] is the equivalent conductivity tensor, kr [−] is relative permeability with respect to the degree of water
(A)
(B)
(C)
(D)
Fig. 1. Physical experiments of (A) saturated high density flow (SHD), and (B) unsaturated high density flow (UHD) after 40 min of dense solute infiltration (Simmons et al., 2002). Number of dense plume fingers versus time for the physical (C) SHD and (D) UHD experiments.
72
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
saturation, h0 [L] is equivalent freshwater head, and ρ [M L−3] and ρ0 [M L−3] are the saltwater and freshwater density, respectively. Solute transport is described by the common advection– dispersion–diffusion equation (Bear and Bachmat, 1990) ∂ðϕcÞ ¼ −∇ ðqc−ϕD∇cÞ ∂t
ð3Þ
where c [−] is a relative solute concentration that varies from 0 to 1, ϕ [−] is porosity, and D [L2 T−1] the hydrodynamic dispersion tensor (Scheidegger, 1954). Water saturation (Sw) and relative permeability (Kr w) in the above equations are characterized by the functions proposed by Van Genuchten (1980) with parameters α = 14.5 m−1, β = 2.68, and Swr = 0.045 that are representative of medium sand (Carsel and Parrish, 1988) as used in the sand tank experiments by Simmons et al. (2002). Geometry of the numerical simulation domain is identical to that of the laboratory experiment, and extensions are lx = 1178 mm, ly = 53 mm and lz = 1060 mm. Left, right and bottom boundaries are assumed to be impermeable for flow. For the SHD case, the model represents ponding due to injection of the solute source by a constant equivalent freshwater head (Dirichlet) boundary condition that is h0 = 4.2 mm for the solute source and h0 = 0 mm on both sides of the source. For the UHD case, a specified flux (Neumann) boundary condition which varies in time is prescribed at the solute source and drain nodes that allow outflow if the water table reaches the top of the domain are implemented at both upper corners of the domain. For the SHD and UHD cases, all boundaries are zero-dispersive flux boundaries except the solute source where the constant relative concentration is c = 1 (corresponding to Cmax = 313,060 mg L−1, and ρmax = 1235 kg m−3). Initial conditions are h0 = 0 mm for saturated flow (SHD), h0 = 620 mm for unsaturated flow (UHD), and c = 0 for both cases. Total simulation time is 1 h, and time step size is adaptive such that relative concentration does not vary by more than 0.01 during a single time step. Dimensions, initial and boundary conditions are shown in Fig. 2, and physical parameter values are listed in Table 1. A grid convergence study was carried out, where the SHD case was simulated on increasingly fine spatial grids using 100 × 100, 200 × 200, 300 × 300, 400 × 400, 500 × 500 and 600 × 600 block elements. The minimum and maximum of the simulated concentration were observed over time for all simulations and found to fall within the realistic range of 0 ≤ c ≤ 1 (no negative concentrations or concentrations larger than 1) for the 500 × 500 grid. The domain of the SHD and UHD problems was therefore discretized by 500 × 500 block elements of size Δx = 2.356 mm, Δ y = 53 mm and Δ z = 2.12 mm. Using three-dimensional elements with only one element in the y-direction is essentially a 2D representation of a sand tank as previously done by others (Schincariol et al., 1994; Stöckl and Houben, 2012). 4. Perturbation methods The laboratory SHD and UHD experiments by Simmons et al. (2002) are numerically re-simulated here using six different numerical perturbation methods discussed below.
Fig. 2. Conceptual model, boundary and initial conditions. Following Simmons et al. (2002), the initial water level for the unsaturated high-density (UHD) case is at z = 620 mm. The domain is fully saturated for the saturated high-density (SHD) case.
The purpose of these methods is to numerically generate random variations in solute concentration or in the conductivity field (K-field) in order to trigger dense plume fingering. Goodness of the fit between physical and numerical results is measured using the number of simulated fingers at different output times (SHD: 10, 20, 30, 40, 50, 60 min; UHD: 6, 12, 15, 20, 26, 32, 40, 46 min). 4.1. Method A: homogeneous sand (reference case) This is the reference case where no numerical perturbation is applied (Fig. 3A). The K-field is homogeneous and the concentration is not perturbed. Concentration in the source is Table 1 Material parameters used in numerical simulations. Parameter
Value
Freshwater density Maximum fluid density Maximum CaCl2 concentration Porosity Hydraulic conductivity Fluid reference viscosity (μ 0)
1000 kg m−3 1235 kg m−3 313,060 mg L−1 36% 0.00163 m s−1 0.001 kg m−1 s−1
6
μ ðcÞ ¼ μ0 ∑ ai ci i¼0
a0 a1 a2 a3 a4 a5 a6 Longitudinal dispersivity Transverse dispersivity Free-solution diffusion coefficient
1.0015947 0.6269767 1.560481 −4.9759998 9.53726 −7.753216 2.4239454 0.001 m 0.0002 m 1.2 × 10−9 m2 s−1
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
(A)
(B)
(C)
(D)
(E)
(F)
73
Fig. 3. Tested perturbation methods: (A) homogeneous sand, (B) initial perturbation at domain center, (C) spatially random, time-constant perturbation of solute source, (D) spatially and temporally random perturbation of solute source, (E) spatially and temporally random noise of solute concentration and (F) random K-field.
74
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
spatially and temporally constant (c = 1). All following methods are modifications of method A. 4.2. Method B: initial perturbation at the domain center The initial concentration at the central nodes (at z = 0 and z = Δ z) of the simulation domain is increased from c = 0 to the small concentration increment XB (Fig. 3B). For the SHD and UHD cases, that increment is modified seven times in separate simulations (XB = 10−7, 10−6, 10−5, 10−4, 10−3, 10−2, 10−1), in order to understand its effect on plume fingering. Method B corresponds to the trigger used in Weatherill et al. (2004) to generate convective rolls in a closed box. Physically, this approach represents a small trace of salt initially present at the central point in the aquifer. The initial increase in salt concentration pre-defines the initial buoyancy-induced downward flow velocity at that point. This local perturbation is thought to trigger instability (finger) development in the entire sand box. 4.3. Method C: spatially random, time-constant perturbation of solute source For each node representing the solute source, a time-constant spatially random concentration is applied using c = 1 + Xc (rand(0) − 0.5) where rand(0) returns a real number between 0 and 1 (Fig. 3C). For the SHD and UHD simulations, the concentration increment Xc is modified seven times in separate simulations (Xc = 10−7, 10−6, 10−5, 10−4, 10−3, 10−2, 10−1). Because of the random character of method C, it is necessary to conduct more than one realization in order to exclude the effect of randomness on simulated fingering. Therefore, five realizations were conducted for each of the chosen trigger values. The idea of method C is to mimic random variations of the source concentration, which may occur in the laboratory experiment. Method C has not yet been applied in prior studies of variabledensity flow. 4.4. Method D: spatially and temporally random perturbation of solute source While method C assumes a time-constant spatially random concentration, method D assumes a time-variable spatially random perturbation to the source concentration in the form c(t) = 1 + XD (rand(0) − 0.5) (Fig. 3D) where rand(0) returns a time-dependent random number between 0 and 1. Five realizations for each of the seven XD values (XD = 10−7, 10−6, 10−5, 10−4, 10−3, 10−2, 10−1) are conducted to eliminate the effect of randomness. Method D has previously been applied by Simmons et al. (1999) and Wooding et al. (1997a) to generate dense plume fingering, but the significance of XD has not previously been studied. Method D is similar to the idea of method C to numerically reproduce physically occurring random concentration variations at the solute source. 4.5. Method E: spatially and temporally random noise of solute concentration
Dirichlet node I is perturbed by the expression cI = cI + XE(rand(0) − 0.5) where concentration increment XE is modified seven times (XE = 10−7, 10−6, 10−5, 10−4, 10−3, 10−2, 10−1) for the SHD and UHD cases in separate simulations. Physically, this approach accounts for constantly occurring, random fluctuations of fluid concentration within a real domain. Method E has not yet been applied in the literature of variable-density flow. 4.6. Method F: random K-field Random hydraulic conductivity fields (random K-fields) are being generated, and it is assumed that the random K variation represents pore-scale heterogeneity of the sand in the laboratory tank (Fig. 3F). The K-fields are being generated with a random K-field generator kindly provided by Olaf Cirpka (Universität Tübingen, Germany). By applying a variance to ln(K), a set of spatially variable K-values is generated where each value is associated with an element in the numerical grid, and where the K-values vary around the mean K-value of 0.00163 m s−1. Variance is chosen such that it induces physically insignificant but numerically significant heterogeneity to trigger plume fingering. Although the field generator offers the possibility to create correlated anisotropic random fields, solely random uncorrelated isotropic fields (i.e. white noise) are used in the present study because the laboratory soil is isotropic and homogeneous on the laboratory (meter) scale. Uncorrelated fields are created by choosing correlation lengths smaller than the cell size of the grid. Therefore, a random white noise of the K-field is applied to perturb variable-density flow. To determine the appropriate variance of ln(K) that reproduces the fingering of the SHD and UHD laboratory experiments, seven K-fields with different ln(K) variances (σln(K) = 10−7, 10−6, 10−5, 10−4, 10−3, 10−2, 10−1) are used to numerically simulate the SHD and UHD experiments. Variance 10−3 gives the best fit between physically and numerically simulated fingers. That variance is used to generate ten realizations in order to eliminate the effect of randomness of the K-field. Simmons et al. (2001) examined the effect of lenticular and layered K-fields on the generation and propagation of dense plume fingers on a large spatial scale (N100 m). It was found that heterogeneity can dissipate fingers once they are created. However, the effect of pore-scale heterogeneity (method F) on the generation of dense plume fingers has not yet been tested in the literature of variable-density flow. 5. Results and discussion The physical SHD and UHD cases were simulated using the six perturbation methods presented in Section 4. To quantify the quality of each perturbation method, errors between physically and numerically simulated number of plume fingers (f phys and f num, respectively) are evaluated in the L2-Norm: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u n 2 uX phys err ¼ t f i −f num i
ð4Þ
i¼1
Instability of flow is triggered by random variations of simulated concentrations on all non-Dirichlet transport nodes in the grid (Fig. 3E). The simulated concentration on the non-
where i is the output time, at which f phys (and therefore also f num) has been measured.
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
5.1. Saturated flow Results for saturated flow (SHD) are shown and interpreted in Figs. 4, 5 and 6. Fig. 4 presents the errors between physically and numerically simulated number of plume fingers evaluated in the L2-Norm for various magnitudes of each perturbation method. For the minimum error of each method, the simulated salt distribution after 40 min is shown in Fig. 5. For each perturbation method, Fig. 6 shows the number of dense plume fingers versus time for perturbation magnitudes giving the largest and smallest errors in the L2-Norm. For comparison, physically simulated number of fingers are also shown in Fig. 6.
75
Clearly, the unperturbed base case (method A) does not adequately reproduce fingering because of the large error (Fig. 4A) and because simulated salt distribution (Fig. 5A) is dissimilar to the real salt distribution (Fig. 1A). Without perturbation (method A), only few fingers are formed throughout the simulation (Fig. 6A), which are likely being created due to large concentration gradients at the edges of the salt source. The result of method A shows no resemblance to the laboratory experiment, where numerous small fingers are initially being generated which later merge into fewer larger fingers. This result suggests that one of the perturbation methods B to F is required to adequately simulate plume fingering.
(A)
(B)
(C)
(D)
(E)
(F)
Fig. 4. SHD case: Error between physically and numerically simulated number of plume fingers in the L2-Norm using perturbation methods (A) to (F) as shown in Fig. 3.
76
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
(A)
(B)
(C)
(D)
(E)
(F)
Fig. 5. SHD case: Results of numerical simulations at t = 40 min for all perturbation mechanisms with best fit using (A) no perturbation, (B) XB = 10−1, (C) XC = 10−2, (D) XD = 10−1, (E) XE = 10−1, (F) σln(K) = 10−3.
Errors for method B are high and do not show a relation to the magnitude of the applied trigger (Fig. 4B). The best fit is achieved for the trigger XB = 10−1. The corresponding simulated salt distribution after 40 min is shown in Fig. 5B. Only few fingers are generated whose amount seems to randomly rise and fall during the numerical simulation as shown in Fig. 6B. Neither a reproduction of the exact amount nor an approximation of the trend of fingers created in the laboratory experiment is achieved. The trigger magnitude for method C has a noticeable effect on the errors in the L2-Norm (Fig. 4C). The trend is that stronger
perturbation decreases the errors (Fig. 4C). The minimum error is found for XC = 10−2, and the corresponding salt distribution is shown in Fig. 5C. The trigger XC = 10−2 of method C is capable of realistically reproducing the decrease of number of salt fingers during the experiment, as evidenced in Fig. 6C. Application of a spatially and temporally perturbation of the solute source (method D) leads to smaller errors when a larger perturbation is applied, consequently giving the smallest error for 10−1 (Fig. 4D). The corresponding simulated salt distribution after 40 min is shown in Fig. 5D. As shown in Fig. 6D, the numerical simulation utilizing perturbation method D is able to
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
adequately reproduce the decrease of number of salt fingers throughout the experiment. Fig. 6D also demonstrates the impact of the perturbation strength, where the laboratory experiments are better reproduced for larger trigger values.
77
Utilizing a spatially and temporally random noise of solute concentration (method E) shows the same trend as method D, where the largest perturbation of 10−1 gives the smallest error (Fig. 4E). The corresponding salt distribution for method E is
(A)
(B)
(C)
(D)
(E)
(F)
Fig. 6. SHD case: Number of dense plume fingers versus time for all perturbation methods applied, comprising the perturbation magnitudes that give the worst and the best fit respectively.
78
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
shown in Fig. 5E. For perturbation method E, large trigger values reproduce the trend observed in the laboratory experiment. However, if a small noise is applied, the number of fingers seems to randomly vary in time during the experiment (Fig. 6E). The use of a random K-field (method F) produces the smallest error for a ln(K)-variance of 10−3 (Fig. 4F). The corresponding salt distribution after 40 min (Fig. 5F) is able to capture the convective features of the physical experiment (Fig. 1A). Fig. 6F shows the number of fingers for a small (10−6) and a larger (10−3) variance of ln(K) versus time. For the ln(K) variance of 10−3, the numbers of fingers formed in the numerical experiment and the physical experiment are in very
good agreement. This observation indicates that using random K-fields in the simulation of dense-plume transport realistically reproduces plume fingering in saturated media. 5.2. Unsaturated flow The results for unsaturated flow (UHD) are shown and interpreted in Figs. 7, 8 and 9. Fig. 7 shows the error between physically and numerically simulated numbers of plume fingers evaluated in the L2-Norm for various magnitudes of each perturbation method. For the minimum error of each method, the simulated salt distribution after 40 min is shown in
(A)
(B)
(C)
(D)
(E)
(F)
Fig. 7. UHD case: Error between physically and numerically simulated number of plume fingers in the L2-Norm using perturbation methods (A) to (F) as shown in Fig. 3.
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
Fig. 8. For each perturbation method, Fig. 9 shows the number of dense plume fingers versus time for perturbation methods giving the largest and smallest errors in the L2-Norm. For comparison, Fig. 9 also contains the number of physically simulated number of fingers. As evidenced by the large error in Fig. 7A and the simulated salt distribution in Fig. 8A, the unperturbed case (method A) does not adequately reproduce fingering observed in the laboratory experiment (Fig. 1B). The numerically simulated salt distribution after 40 min (Fig. 8A) only produces two large
79
fingers probably due to large concentration gradients at the edges of the solute source. It is therefore necessary to introduce one of the perturbation methods B to F to the numerical simulation. When an initial perturbation at the domain center is induced (method B), the lowest error is obtained with a perturbation magnitude of 10−1 (Fig. 7B). However, errors between the physical and the numerical experiment remain large for this trigger. Compared to the unperturbed sand (method A), the application of method B leads to the formation
(A)
(B)
(C)
(D)
(E)
(F)
Fig. 8. UHD case: Results of numerical simulations at t = 40 min for all perturbation mechanisms with best fit using (A) no perturbation, (B) XB = 10−1, (C) XC = 10−1, (D) XD = 10−1, (E) XE = 10−1, (F) σln(K) = 10−1.
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
20 18 16 14 12 10 8 6 4 2 0
(A)
20
30 40 time [min]
10
20
30 40 time [min]
20 18 16 14 12 10 8 6 4 2 0
50
50
20
20
50
30 40 time [min]
50
laboratory σln(K) =1E-1 σln(K) =1E-5
(F)
10
30 40 time [min] laboratory XD =1E-1 XD =1E-5
(D)
10
20 18 16 14 12 10 8 6 4 2 0
laboratory XB =1E-1 XB =1E-7
(B)
10
laboratory XE =1E-1 XE =1E-5
(E)
20 18 16 14 12 10 8 6 4 2 0
50
number of fingers
number of fingers
30 40 time [min] laboratory XC =1E-1 XC =1E-4
(C)
10
20 18 16 14 12 10 8 6 4 2 0
20
number of fingers
number of fingers
10
20 18 16 14 12 10 8 6 4 2 0
laboratory unperturbed number of fingers
number of fingers
80
20
30 40 time [min]
50
Fig. 9. UHD case: Number of dense plume fingers versus time for all perturbation methods applied, comprising the perturbation magnitudes that give the worst and the best fit respectively.
of more fingers after 40 min (Fig. 8B). However, the number of fingers produced by method B does not adequately reproduce results observed in the physical experiment (Fig. 1B). Perturbation methods C and D apply perturbations directly at the solute source on top of the unsaturated zone. For both, methods C and D, the smallest error is found for a trigger value
of 10−1 (Fig. 7C, D). However, it should be noted that the trigger value does not have a large influence on the error. Methods C and D result in the formation of more plume fingers in the central region of the plume front (Fig. 8C, D) in comparison to method A, but are not able to reproduce the finger pattern observed in the physical experiment. The amount of fingers
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
formed throughout the simulation for methods C and D is shown in Fig. 9C and D. Fig. 7E shows that a spatially and temporally random noise of solute concentration (method E) has a substantial impact on simulated plume fingering if the trigger magnitude is large. However, with increasing trigger magnitude errors substantially decrease (Fig. 7E). The best agreement between the physical and the numerical experiment is obtained with the largest trigger magnitude of 10−1. For this trigger, the numerically simulated salt distribution after 40 min is shown in Fig. 8E. The number of fingers produced by method E is larger than those of perturbation methods A to D, and it is almost identical to the number of fingers generated in the physical experiment (Figs. 1B and 8E). Perturbation method F (random K-field) gives the smallest error at a trigger magnitude of 10−1 (Fig. 7F). The associated plume fingering is shown in Fig. 8F. Similar to method E, use of a random K-field in method F realistically reproduces dense plume fingering under unsaturated conditions (Figs. 1B and 8F). Fig. 9F shows the number of dense plume fingers created with a random K-field. The trend in the amount of fingers produced in the physical experiment is in agreement with observations from the physical experiment. Both physical and numerical observations show an initially increasing and afterwards decreasing number of fingers. It should be noted that perturbation with a random K-field even for small variances (10−5, Fig. 9F) significantly improves the simulated
81
transport behavior and plume finger development compared to the unperturbed homogeneous reference case (Fig. 9A). 6. Summary and conclusion Perturbation mechanisms have been frequently used in recent literature to induce instabilities in the numerical representation of variable-density systems (Schincariol, 1998; Schincariol et al., 1997; Simmons et al., 1999; Weatherill et al., 2004; Xie et al., 2011, 2012). However, studies typically focused on a particular perturbation mechanism (e.g. Simmons et al. (1999), Weatherill et al. (2004), Xie et al. (2012)). Therefore, this study included various perturbation mechanisms documented in the literature of variable-density flow. Furthermore it has been extended by perturbation mechanisms that have not been used in the context of variable-density flow before. A quantitative investigation of the sensitivity of trigger magnitudes and a systematic evaluation of perturbation methods are still absent in the literature. Therefore, this study is aimed at systematically evaluating and comparing mechanisms to realistically create plume fingering. A laboratory experiment has been numerically re-simulated and different perturbation methods were applied. The quality of numerical results and hence quality of individual perturbation methods were evaluated by comparing simulated numerical results with laboratory results obtained by Simmons et al. (2002). Comparison was conducted by quantifying the number
(A)
(B)
(C)
(D)
Fig. 10. Results at t = 40 min for SHD (left column) and UHD (right column) physical (upper row) and numerical (bottom row) experiments (Simmons et al., 2002) for perturbation methods and trigger magnitudes giving the smallest error in the L2-Norm (Figs. 4, 7).
82
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
of dense plume fingers formed during the experiment and by visual comparison of salt distribution in the domain. To the authors' knowledge, this is the first study that compares and evaluates multiple perturbation methods for the simulation of saturated and unsaturated density-driven flow and solute transport in porous media. The main new key findings of this study are: 1. Perturbation is necessary to realistically reproduce dense plume migration under saturated and unsaturated conditions. 2. While a spatially random, time-constant perturbation of the solute source (method C) is appropriate to simulate plume fingers under saturated conditions (Fig. 10), that method fails under unsaturated conditions. 3. Plume fingering under saturated–unsaturated flow conditions can be reproduced best with a spatially and temporally random noise of solute concentration (method E, Fig. 10) and with a random K-field (method F). It can be concluded that perturbation is needed in numerical models to realistically reproduce dense plume migration. While all of the evaluated perturbation methods have an impact on the migration of a dense plume, only some methods were able to give a good approximation of real plume migration simulated in the laboratory experiment. It can also be concluded that pore-scale heterogeneity (represented by methods E and F) is the main driving mechanism for dense plume fingering in homogeneous porous media. The findings of this study can serve as a basis for later benchmarking of density-driven flow in porous media under saturated and unsaturated conditions. Additionally, it constitutes a progress for the detailed realistic representation of dense plume fingering. Future work could focus on the applicability of perturbation mechanisms to different field-scale problems and the use of field data to assess the quality of various perturbation mechanisms. Acknowledgment We are grateful to C.T. Simmons (Flinders University, Australia) for providing photographic material and for providing input on the laboratory experiment, and to O. Cirpka (Universität Tübingen, Germany) for providing a random field generator and advice concerning its use. We thank the editorial board of the Journal of Contaminant Hydrology (B. Sleep) for handling our manuscript. Furthermore, we are grateful for the constructive comments by the two anonymous reviewers which helped to improve the manuscript. References Bear, J., Bachmat, Y., 1990. Introduction to Modeling of Transport Phenomena in Porous Media vol. 4. Springer. Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24 (5), 755–769.
Diersch, H.J.G., Kolditz, O., 2002. High-density flow and transport in porous media: approaches and challenges. Adv. Water Resour. 25, 899–944. Goswami, R.R., Clement, T.P., Hayworth, J.H., 2011. Comparison of numerical techniques used for simulating variable-density flow and transport experiments. J. Hydrol. Eng. 17 (2), 272–282. Horton, C.W., Rogers, F.T., 1945. Convection currents in a porous medium. J. Appl. Phys. 16, 367–370. Illangasekare, T., Tyler, S.W., Clement, T.P., Villholth, K.G., Perera, A.P.G.R.L., Obeysekera, J., Jensen, K., 2006. Impacts of the 2004 tsunami on groundwater resources in Sri Lanka. Water Resour. Res. 42 (5). Lapwood, E.R., 1948. Convection of a fluid in a porous medium. Proc. Camb. Philol. Soc. 44, 508–521. List, E.J., 1965. The Stability and Mixing of a Density-stratified Horizontal Flow in a Saturated Porous Medium. (Dissertation (Ph.D.)). California Institute of Technology. Scheidegger, A.E., 1954. Statistical hydrodynamics in porous media. J. Appl. Phys. 25, 994. Schincariol, R.A., 1998. Dispersive mixing dynamics of dense miscible plumes: natural perturbation initiation by local-scale heterogeneities. J. Contam. Hydrol. 34, 247–271. Schincariol, R.A., Schwartz, F.W., 1990. An experimental investigation of variable density flow and mixing in homogeneous and heterogeneous media. Water Resour. Res. 26 (10), 2317–2329. Schincariol, R.A., Schwartz, F.W., Mendoza, C.A., 1994. On the generation of instabilities in variable density flow. Water Resour. Res. 30 (4), 913–927. Schincariol, R.A., Schwartz, F.W., Mendoza, C.A., 1997. Instabilities in variable density flows: stability and sensitivity analyses for homogeneous and heterogeneous media. Water Resour. Res. 33 (1), 31–41. Simmons, C.T., 2005. Variable density groundwater flow: from current challenges to future possibilities. Hydrogeol. J. 13, 116–119. Simmons, C.T., Narayan, K.A., Wooding, R.A., 1999. On a test case for densitydependent groundwater flow and solute transport models: the salt lake problem. Water Resour. Res. 35 (12), 3607–3620. Simmons, C.T., Fenstemaker, T.R., Sharp, J.M., 2001. Variable-density groundwater flow and solute transport in heterogeneous porous media: approaches, resolutions and future challenges. J. Contam. Hydrol. 52, 245–275. Simmons, C.T., Pierini, M.L., Hutson, J.L., 2002. Laboratory investigation of variable-density flow and solute transport in unsaturated–saturated porous media. Transp. Porous Media 47, 215–244. Stöckl, L., Houben, G., 2012. Flow dynamics and age stratification of freshwater lenses: experiments and modeling. J. Hydrol. 458–459, 9–15. Therrien, R., McLaren, R.G., Sudicky, E.A., Panday, S.M., 2012. HydroGeoSphere — A Three-dimensional Numerical Model Describing Fully-integrated Subsurface and Surface Flow and Solute Transport. Groundwater Simulations Group, University of Waterloo, Canada. Van Dam, R.L., Simmons, C.T., Hyndman, D.W., Wood, W.W., 2009. Natural free convection in porous media: first field documentation in groundwater. Geophys. Res. Lett. 36 (11). Van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. Voss, C.I., Souza, W.R., 1987. Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater–saltwater transition zone. Water Resour. Res. 23, 1851–1866. Weatherill, D., Simmons, C.T., Voss, C.I., Robinson, N.I., 2004. Testing densitydependent groundwater models: two-dimensional steady state unstable convection in infinite, finite and inclined porous layers. Adv. Water Resour. 27, 547–562. Wooding, R.A., Tyler, S.W., White, I., Anderson, P.A., 1997a. Convection in groundwater below an evaporating salt lake: 1. Onset of instability. Water Resour. Res. 33 (6), 1199–1217. Wooding, R.A., Tyler, S.W., White, I., Anderson, P.A., 1997b. Convection in groundwater below an evaporating salt lake: 2. Evolution of fingers or plumes. Water Resour. Res. 33 (6), 1219–1228. Xie, Y., Simmons, C.T., Werner, A.D., 2011. Speed of free convective fingering in porous media. Water Resour. Res. 47 (11). Xie, Y., Simmons, C.T., Werner, A.D., Diersch, H.J.G., 2012. Prediction and uncertainty of free convection phenomena in porous media. Water Resour. Res. 48 (2).
Contents lists available at ScienceDirect
Journal of Contaminant Hydrology journal homepage: www.elsevier.com/locate/jconhyd
Generation of dense plume fingers in saturated–unsaturated homogeneous porous media Clemens J.M. Cremer ⁎, Thomas Graf Institute of Fluid Mechanics and Environmental Physics in Civil Engineering, Leibniz Universität Hannover, Appelstr. 9A, 30167 Hannover, Germany
a r t i c l e
i n f o
Article history: Received 7 July 2014 Received in revised form 10 November 2014 Accepted 26 November 2014 Available online 9 December 2014 Keywords: Model Groundwater Saturated Unsaturated Density Numerical trigger Fingering
a b s t r a c t Flow under variable-density conditions is widespread, occurring in geothermal reservoirs, at waste disposal sites or due to saltwater intrusion. The migration of dense plumes typically results in the formation of vertical plume fingers which are known to be triggered by material heterogeneity or by variations in source concentration that causes the density variation. Using a numerical groundwater model, six perturbation methods are tested under saturated and unsaturated flow conditions to mimic heterogeneity and concentration variations on the pore scale in order to realistically generate dense fingers. A laboratory-scale sand tank experiment is numerically simulated, and the perturbation methods are evaluated by comparing plume fingers obtained from the laboratory experiment with numerically simulated fingers. Dense plume fingering for saturated flow can best be reproduced with a spatially random, time-constant perturbation of the solute source. For unsaturated flow, a spatially and temporally random noise of solute concentration or a random conductivity field adequately simulate plume fingering. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Leachate that invades soil, e.g. from waste disposal sites, frequently possesses higher solute concentration and therefore higher fluid density than ambient freshwater. Spatial variation of water density can lead to variable-density flow, which is a key process in the transport of contaminants, salt transport in coastal aquifers or hot geothermal plumes (Diersch and Kolditz, 2002; Illangasekare et al., 2006; Simmons, 2005; Simmons et al., 2001; Stöckl and Houben, 2012; Van Dam et al., 2009). The migration of dense plumes through a porous medium typically results in the formation of lobe-shaped instabilities or vertical plume fingers, which are known to have a nonnegligible effect on the propagation of the plume (Xie et al., 2011). Flow direction in solute plume fingers is downwards, counterbalanced by an upward flow of less dense fluid between ⁎ Corresponding author. E-mail address: [email protected] (C. Cremer). URL: http://www.hydromech.uni-hannover.de (C. Cremer).
http://dx.doi.org/10.1016/j.jconhyd.2014.11.008 0169-7722/© 2014 Elsevier B.V. All rights reserved.
the fingers. The generation of plume fingers can be attributed to material heterogeneity, which is known to trigger the formation of dense fingers (Goswami et al., 2011; Schincariol, 1998; Schincariol and Schwartz, 1990; Schincariol et al., 1997). In homogeneous media, however, fingers are also created even though the material is uniform with identical sizes of material grains (e.g. sand grains) (Goswami et al., 2011; Simmons et al., 2002; Van Dam et al., 2009; Wooding et al., 1997b). The reason can be assumed to be that pore-scale heterogeneity leading to different flow velocities also exists in homogeneous media due to two effects: (i) random variations of solute concentration leading to varying buoyancy effects, which results in different average flow velocities, and (ii) grains of identical size which may randomly arrange differently, thus forming tetrahedrons, hexahedrons or octahedrons, where each arrangement creates pores of varying diameters, thus also resulting in different pore-scale velocities. Usually, variable-density flow models do not consider effects (i) and (ii), such that the numerical model does not apply any external perturbation to trigger plume fingering. In that case, numerical errors of the computing process are expected to
70
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
perturb the flow field and to create unstable plume fingers. The problem with this approach, however, is that “these instabilities are not physically realistic and are essentially uncontrollable because their character depends on the extent to which numerical errors develop” (Schincariol et al., 1994). This has later been confirmed by Simmons et al. (1999) who noted that numerical errors “may occur by themselves and are uncontrollable as location and magnitude of the error may vary with time. Small changes to dispersion parameters, time step sizes, or spatial discretization throughout a series of simulations can cause different instabilities to form.” Xie et al. (2011) added that numerical errors are not easy to quantify. As an alternative to relying on internal numerical errors, deliberately triggering instabilities by external perturbations represents a controlled way to account for pore-scale heterogeneity. Perturbations can be introduced through a noise applied to the solute concentration or through a random conductivity field. List (1965) induced perturbations of the source concentration created by superposition of spatially oscillatory terms with sinusoidal form and performed a theoretical analysis on the relationship between the nature of the perturbation and the resulting instabilities. List (1965) identified the following parameters responsible for instability growth: the transverse Rayleigh number, the ratio of the transverse to the longitudinal Rayleigh numbers, and the nondimensional wave number. However, while List's stability theory provides a useful model for determining stability in a variable-density system, Schincariol (1998) criticized that the theory cannot account for dispersive dissipation of the dense plume as a function of travel distance and is therefore less useful in predicting the growth rate of instabilities in real systems. Schincariol and Schwartz (1990) conducted physical experiments on the laboratory-scale in homogeneous, layered, and lenticular porous media in order to further investigate the influence of hydraulic conductivity and density differences on the formation of instabilities. It was observed that small density differences around 0.8 mg cm−3 between a higher density water overlying a less dense water results in the formation of gravitational instabilities. Schincariol et al. (1994) numerically reproduced the physical experiment by Schincariol and Schwartz (1990). Numerically generated errors and their influence on the development of instabilities were examined using the Courant- and Peclet-stability criteria. Schincariol et al. (1994) concluded that the formation of unstable flow depends on the Rayleigh-number and on the wavelength of the perturbation. Schincariol et al. (1997) and Schincariol (1998) tested the applicability of the stability theory developed by List (1965) to variable-density flow in homogeneous and heterogeneous media. Numerical simulations were carried out, where perturbations were induced by lab-scale heterogeneity of the conductivity field. Schincariol et al. (1997) concluded that besides the Raleigh number and wavelength, statistical characteristics of the conductivity field (mean, variance, correlation length scales) play an important role in the development of instabilities. Simmons et al. (1999) numerically resimulated an experiment of an evaporating salt lake developed by Wooding et al. (1997a,b). A spatially and temporally random noise of the salt concentration in the lake was used to apply perturbation and to mimic field-scale heterogeneity. Using that perturbation technique, Simmons et al. (1999) could realistically reproduce
experimental results by Wooding et al. (1997a,b). Xie et al. (2011) adopted this perturbation mechanism to study the effects of dynamic solute boundary conditions on the migration of dense plumes using the modified solute Elder problem (Voss and Souza, 1987). That problem was later reused and slightly modified in multiple realizations to analyze the speed of free convective fingering (Xie et al., 2012). Weatherill et al. (2004) investigated the onset of convective roll instabilities due to thermal convection in a fluid-saturated porous layer that is uniformly being heated from below. The Horton–Rogers–Lapwood (HRL) problem (Horton and Rogers, 1945; Lapwood, 1948) was numerically resimulated with an initial perturbation at the domain center to initiate convection. A detailed investigation of that perturbation method is absent, and it therefore remains unknown whether the trigger used by Weatherill et al. (2004) really affects the onset of convection. Recently, Illangasekare et al. (2006) and Van Dam et al. (2009) obtained evidence of lobe-shaped instabilities from field measurements. This highlights the importance of unstable density-driven flow in natural field settings and furthermore the importance of an adequate representation of these phenomena in numerical models. The objective of the present study is to evaluate six different perturbation mechanisms in order to incorporate pore-scale heterogeneity into a numerical model such that dense fingers are realistically being generated under saturated and unsaturated flow conditions. While partly relying on the evaluation of well known perturbation mechanisms, this study also introduces several mechanisms that are new to variable-density flow simulations. To the authors' knowledge, (i) a spatially and temporally random noise of solute concentration, (ii) the application of a spatially random, time-constant perturbation of the solute source, and (iii) a random conductivity field have previously not been examined as perturbation mechanisms in the literature of variable-density flow. Success of all perturbation methods is validated by comparison of numerical results with a physical laboratory model conducted by (Simmons et al., 2002). 2. Physical model Simmons et al. (2002) conducted laboratory experiments of variable-density flow and salt transport in a 2D vertical sand tank. Ambient horizontal flow was not considered. The focus was put on variable-density free convective flow including the generation of dense salt fingers. Stained calcium chloride solutions of different densities were injected along the central part of the upper boundary. The experiments were carried out under fully saturated as well as partly saturated conditions. For the saturated case, injection of dense solutes created ponding along the source, such that the water table in the source was a few millimeters higher than those along the left and right top boundaries. The unsaturated experiment was carried out by draining the water level to 44 cm below the top surface of the sand. The present study focuses on the saturated high density (SHD) case and the unsaturated high density (UHD) case (Simmons et al., 2002) where the density of the source concentration was 1235 kg m−3. Fig. 1A and B presents the plume development of the physical SHD and UHD experiments after 40 min of dense
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
plume fingering. The density difference between freshwater initially filling the tank and the dense calcium chloride solution leads to downward density-driven migration within distinct plume fingers. Downward movement within these salt fingers is counterbalanced by upward movement of freshwater between salt fingers. In the UHD case, fingering did not take place in the unsaturated zone because of the small unsaturated hydraulic conductivity. Note that, in the UHD case, fingering started at the transition between the unsaturated and the saturated zones. Fig. 1C and D also show the number of fingers versus time for the SHD and UHD cases (Simmons et al., 2002). The number of plume fingers for the SHD consistently decreases with time. However, the UHD case shows an initial increasing and afterwards a decreasing number of fingers. According to Simmons et al. (2002), the time when plume fingers become more numerous coincides with the dense plume reaching the water table. While flow is dominated by capillarity in the unsaturated zone, flow in the saturated zone is dominated by density differences. The decrease in number of fingers can be attributed to the fact that multiple small fingers merge into smaller fingers over time (Fig. 1D). The evolution of the number of fingers with time will be used here to test different perturbation mechanisms. More
71
details on and an interpretation of the physical model can be found in Simmons et al. (2002). 3. Numerical model The numerical variably saturated variable-density flow and solute transport model HydroGeoSphere (Therrien et al., 2012) is used here to re-simulate the SHD and UHD laboratory experiments conducted by Simmons et al. (2002). Transient flow through a variably saturated porous medium is described by the Richards' equation (Therrien et al., 2012) ∂ðθS Sw Þ ¼ −∇q ∂t
ð1Þ
where θs [−] is saturated water content and Sw [−] is water saturation. The specific discharge q [L T−1] is given by the Darcy equation ρ−ρ0 ∇z q ¼ −K 0 kr ∇ h0 ρ0
ð2Þ
where K0 [L T−1] is the equivalent conductivity tensor, kr [−] is relative permeability with respect to the degree of water
(A)
(B)
(C)
(D)
Fig. 1. Physical experiments of (A) saturated high density flow (SHD), and (B) unsaturated high density flow (UHD) after 40 min of dense solute infiltration (Simmons et al., 2002). Number of dense plume fingers versus time for the physical (C) SHD and (D) UHD experiments.
72
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
saturation, h0 [L] is equivalent freshwater head, and ρ [M L−3] and ρ0 [M L−3] are the saltwater and freshwater density, respectively. Solute transport is described by the common advection– dispersion–diffusion equation (Bear and Bachmat, 1990) ∂ðϕcÞ ¼ −∇ ðqc−ϕD∇cÞ ∂t
ð3Þ
where c [−] is a relative solute concentration that varies from 0 to 1, ϕ [−] is porosity, and D [L2 T−1] the hydrodynamic dispersion tensor (Scheidegger, 1954). Water saturation (Sw) and relative permeability (Kr w) in the above equations are characterized by the functions proposed by Van Genuchten (1980) with parameters α = 14.5 m−1, β = 2.68, and Swr = 0.045 that are representative of medium sand (Carsel and Parrish, 1988) as used in the sand tank experiments by Simmons et al. (2002). Geometry of the numerical simulation domain is identical to that of the laboratory experiment, and extensions are lx = 1178 mm, ly = 53 mm and lz = 1060 mm. Left, right and bottom boundaries are assumed to be impermeable for flow. For the SHD case, the model represents ponding due to injection of the solute source by a constant equivalent freshwater head (Dirichlet) boundary condition that is h0 = 4.2 mm for the solute source and h0 = 0 mm on both sides of the source. For the UHD case, a specified flux (Neumann) boundary condition which varies in time is prescribed at the solute source and drain nodes that allow outflow if the water table reaches the top of the domain are implemented at both upper corners of the domain. For the SHD and UHD cases, all boundaries are zero-dispersive flux boundaries except the solute source where the constant relative concentration is c = 1 (corresponding to Cmax = 313,060 mg L−1, and ρmax = 1235 kg m−3). Initial conditions are h0 = 0 mm for saturated flow (SHD), h0 = 620 mm for unsaturated flow (UHD), and c = 0 for both cases. Total simulation time is 1 h, and time step size is adaptive such that relative concentration does not vary by more than 0.01 during a single time step. Dimensions, initial and boundary conditions are shown in Fig. 2, and physical parameter values are listed in Table 1. A grid convergence study was carried out, where the SHD case was simulated on increasingly fine spatial grids using 100 × 100, 200 × 200, 300 × 300, 400 × 400, 500 × 500 and 600 × 600 block elements. The minimum and maximum of the simulated concentration were observed over time for all simulations and found to fall within the realistic range of 0 ≤ c ≤ 1 (no negative concentrations or concentrations larger than 1) for the 500 × 500 grid. The domain of the SHD and UHD problems was therefore discretized by 500 × 500 block elements of size Δx = 2.356 mm, Δ y = 53 mm and Δ z = 2.12 mm. Using three-dimensional elements with only one element in the y-direction is essentially a 2D representation of a sand tank as previously done by others (Schincariol et al., 1994; Stöckl and Houben, 2012). 4. Perturbation methods The laboratory SHD and UHD experiments by Simmons et al. (2002) are numerically re-simulated here using six different numerical perturbation methods discussed below.
Fig. 2. Conceptual model, boundary and initial conditions. Following Simmons et al. (2002), the initial water level for the unsaturated high-density (UHD) case is at z = 620 mm. The domain is fully saturated for the saturated high-density (SHD) case.
The purpose of these methods is to numerically generate random variations in solute concentration or in the conductivity field (K-field) in order to trigger dense plume fingering. Goodness of the fit between physical and numerical results is measured using the number of simulated fingers at different output times (SHD: 10, 20, 30, 40, 50, 60 min; UHD: 6, 12, 15, 20, 26, 32, 40, 46 min). 4.1. Method A: homogeneous sand (reference case) This is the reference case where no numerical perturbation is applied (Fig. 3A). The K-field is homogeneous and the concentration is not perturbed. Concentration in the source is Table 1 Material parameters used in numerical simulations. Parameter
Value
Freshwater density Maximum fluid density Maximum CaCl2 concentration Porosity Hydraulic conductivity Fluid reference viscosity (μ 0)
1000 kg m−3 1235 kg m−3 313,060 mg L−1 36% 0.00163 m s−1 0.001 kg m−1 s−1
6
μ ðcÞ ¼ μ0 ∑ ai ci i¼0
a0 a1 a2 a3 a4 a5 a6 Longitudinal dispersivity Transverse dispersivity Free-solution diffusion coefficient
1.0015947 0.6269767 1.560481 −4.9759998 9.53726 −7.753216 2.4239454 0.001 m 0.0002 m 1.2 × 10−9 m2 s−1
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
(A)
(B)
(C)
(D)
(E)
(F)
73
Fig. 3. Tested perturbation methods: (A) homogeneous sand, (B) initial perturbation at domain center, (C) spatially random, time-constant perturbation of solute source, (D) spatially and temporally random perturbation of solute source, (E) spatially and temporally random noise of solute concentration and (F) random K-field.
74
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
spatially and temporally constant (c = 1). All following methods are modifications of method A. 4.2. Method B: initial perturbation at the domain center The initial concentration at the central nodes (at z = 0 and z = Δ z) of the simulation domain is increased from c = 0 to the small concentration increment XB (Fig. 3B). For the SHD and UHD cases, that increment is modified seven times in separate simulations (XB = 10−7, 10−6, 10−5, 10−4, 10−3, 10−2, 10−1), in order to understand its effect on plume fingering. Method B corresponds to the trigger used in Weatherill et al. (2004) to generate convective rolls in a closed box. Physically, this approach represents a small trace of salt initially present at the central point in the aquifer. The initial increase in salt concentration pre-defines the initial buoyancy-induced downward flow velocity at that point. This local perturbation is thought to trigger instability (finger) development in the entire sand box. 4.3. Method C: spatially random, time-constant perturbation of solute source For each node representing the solute source, a time-constant spatially random concentration is applied using c = 1 + Xc (rand(0) − 0.5) where rand(0) returns a real number between 0 and 1 (Fig. 3C). For the SHD and UHD simulations, the concentration increment Xc is modified seven times in separate simulations (Xc = 10−7, 10−6, 10−5, 10−4, 10−3, 10−2, 10−1). Because of the random character of method C, it is necessary to conduct more than one realization in order to exclude the effect of randomness on simulated fingering. Therefore, five realizations were conducted for each of the chosen trigger values. The idea of method C is to mimic random variations of the source concentration, which may occur in the laboratory experiment. Method C has not yet been applied in prior studies of variabledensity flow. 4.4. Method D: spatially and temporally random perturbation of solute source While method C assumes a time-constant spatially random concentration, method D assumes a time-variable spatially random perturbation to the source concentration in the form c(t) = 1 + XD (rand(0) − 0.5) (Fig. 3D) where rand(0) returns a time-dependent random number between 0 and 1. Five realizations for each of the seven XD values (XD = 10−7, 10−6, 10−5, 10−4, 10−3, 10−2, 10−1) are conducted to eliminate the effect of randomness. Method D has previously been applied by Simmons et al. (1999) and Wooding et al. (1997a) to generate dense plume fingering, but the significance of XD has not previously been studied. Method D is similar to the idea of method C to numerically reproduce physically occurring random concentration variations at the solute source. 4.5. Method E: spatially and temporally random noise of solute concentration
Dirichlet node I is perturbed by the expression cI = cI + XE(rand(0) − 0.5) where concentration increment XE is modified seven times (XE = 10−7, 10−6, 10−5, 10−4, 10−3, 10−2, 10−1) for the SHD and UHD cases in separate simulations. Physically, this approach accounts for constantly occurring, random fluctuations of fluid concentration within a real domain. Method E has not yet been applied in the literature of variable-density flow. 4.6. Method F: random K-field Random hydraulic conductivity fields (random K-fields) are being generated, and it is assumed that the random K variation represents pore-scale heterogeneity of the sand in the laboratory tank (Fig. 3F). The K-fields are being generated with a random K-field generator kindly provided by Olaf Cirpka (Universität Tübingen, Germany). By applying a variance to ln(K), a set of spatially variable K-values is generated where each value is associated with an element in the numerical grid, and where the K-values vary around the mean K-value of 0.00163 m s−1. Variance is chosen such that it induces physically insignificant but numerically significant heterogeneity to trigger plume fingering. Although the field generator offers the possibility to create correlated anisotropic random fields, solely random uncorrelated isotropic fields (i.e. white noise) are used in the present study because the laboratory soil is isotropic and homogeneous on the laboratory (meter) scale. Uncorrelated fields are created by choosing correlation lengths smaller than the cell size of the grid. Therefore, a random white noise of the K-field is applied to perturb variable-density flow. To determine the appropriate variance of ln(K) that reproduces the fingering of the SHD and UHD laboratory experiments, seven K-fields with different ln(K) variances (σln(K) = 10−7, 10−6, 10−5, 10−4, 10−3, 10−2, 10−1) are used to numerically simulate the SHD and UHD experiments. Variance 10−3 gives the best fit between physically and numerically simulated fingers. That variance is used to generate ten realizations in order to eliminate the effect of randomness of the K-field. Simmons et al. (2001) examined the effect of lenticular and layered K-fields on the generation and propagation of dense plume fingers on a large spatial scale (N100 m). It was found that heterogeneity can dissipate fingers once they are created. However, the effect of pore-scale heterogeneity (method F) on the generation of dense plume fingers has not yet been tested in the literature of variable-density flow. 5. Results and discussion The physical SHD and UHD cases were simulated using the six perturbation methods presented in Section 4. To quantify the quality of each perturbation method, errors between physically and numerically simulated number of plume fingers (f phys and f num, respectively) are evaluated in the L2-Norm: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u n 2 uX phys err ¼ t f i −f num i
ð4Þ
i¼1
Instability of flow is triggered by random variations of simulated concentrations on all non-Dirichlet transport nodes in the grid (Fig. 3E). The simulated concentration on the non-
where i is the output time, at which f phys (and therefore also f num) has been measured.
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
5.1. Saturated flow Results for saturated flow (SHD) are shown and interpreted in Figs. 4, 5 and 6. Fig. 4 presents the errors between physically and numerically simulated number of plume fingers evaluated in the L2-Norm for various magnitudes of each perturbation method. For the minimum error of each method, the simulated salt distribution after 40 min is shown in Fig. 5. For each perturbation method, Fig. 6 shows the number of dense plume fingers versus time for perturbation magnitudes giving the largest and smallest errors in the L2-Norm. For comparison, physically simulated number of fingers are also shown in Fig. 6.
75
Clearly, the unperturbed base case (method A) does not adequately reproduce fingering because of the large error (Fig. 4A) and because simulated salt distribution (Fig. 5A) is dissimilar to the real salt distribution (Fig. 1A). Without perturbation (method A), only few fingers are formed throughout the simulation (Fig. 6A), which are likely being created due to large concentration gradients at the edges of the salt source. The result of method A shows no resemblance to the laboratory experiment, where numerous small fingers are initially being generated which later merge into fewer larger fingers. This result suggests that one of the perturbation methods B to F is required to adequately simulate plume fingering.
(A)
(B)
(C)
(D)
(E)
(F)
Fig. 4. SHD case: Error between physically and numerically simulated number of plume fingers in the L2-Norm using perturbation methods (A) to (F) as shown in Fig. 3.
76
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
(A)
(B)
(C)
(D)
(E)
(F)
Fig. 5. SHD case: Results of numerical simulations at t = 40 min for all perturbation mechanisms with best fit using (A) no perturbation, (B) XB = 10−1, (C) XC = 10−2, (D) XD = 10−1, (E) XE = 10−1, (F) σln(K) = 10−3.
Errors for method B are high and do not show a relation to the magnitude of the applied trigger (Fig. 4B). The best fit is achieved for the trigger XB = 10−1. The corresponding simulated salt distribution after 40 min is shown in Fig. 5B. Only few fingers are generated whose amount seems to randomly rise and fall during the numerical simulation as shown in Fig. 6B. Neither a reproduction of the exact amount nor an approximation of the trend of fingers created in the laboratory experiment is achieved. The trigger magnitude for method C has a noticeable effect on the errors in the L2-Norm (Fig. 4C). The trend is that stronger
perturbation decreases the errors (Fig. 4C). The minimum error is found for XC = 10−2, and the corresponding salt distribution is shown in Fig. 5C. The trigger XC = 10−2 of method C is capable of realistically reproducing the decrease of number of salt fingers during the experiment, as evidenced in Fig. 6C. Application of a spatially and temporally perturbation of the solute source (method D) leads to smaller errors when a larger perturbation is applied, consequently giving the smallest error for 10−1 (Fig. 4D). The corresponding simulated salt distribution after 40 min is shown in Fig. 5D. As shown in Fig. 6D, the numerical simulation utilizing perturbation method D is able to
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
adequately reproduce the decrease of number of salt fingers throughout the experiment. Fig. 6D also demonstrates the impact of the perturbation strength, where the laboratory experiments are better reproduced for larger trigger values.
77
Utilizing a spatially and temporally random noise of solute concentration (method E) shows the same trend as method D, where the largest perturbation of 10−1 gives the smallest error (Fig. 4E). The corresponding salt distribution for method E is
(A)
(B)
(C)
(D)
(E)
(F)
Fig. 6. SHD case: Number of dense plume fingers versus time for all perturbation methods applied, comprising the perturbation magnitudes that give the worst and the best fit respectively.
78
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
shown in Fig. 5E. For perturbation method E, large trigger values reproduce the trend observed in the laboratory experiment. However, if a small noise is applied, the number of fingers seems to randomly vary in time during the experiment (Fig. 6E). The use of a random K-field (method F) produces the smallest error for a ln(K)-variance of 10−3 (Fig. 4F). The corresponding salt distribution after 40 min (Fig. 5F) is able to capture the convective features of the physical experiment (Fig. 1A). Fig. 6F shows the number of fingers for a small (10−6) and a larger (10−3) variance of ln(K) versus time. For the ln(K) variance of 10−3, the numbers of fingers formed in the numerical experiment and the physical experiment are in very
good agreement. This observation indicates that using random K-fields in the simulation of dense-plume transport realistically reproduces plume fingering in saturated media. 5.2. Unsaturated flow The results for unsaturated flow (UHD) are shown and interpreted in Figs. 7, 8 and 9. Fig. 7 shows the error between physically and numerically simulated numbers of plume fingers evaluated in the L2-Norm for various magnitudes of each perturbation method. For the minimum error of each method, the simulated salt distribution after 40 min is shown in
(A)
(B)
(C)
(D)
(E)
(F)
Fig. 7. UHD case: Error between physically and numerically simulated number of plume fingers in the L2-Norm using perturbation methods (A) to (F) as shown in Fig. 3.
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
Fig. 8. For each perturbation method, Fig. 9 shows the number of dense plume fingers versus time for perturbation methods giving the largest and smallest errors in the L2-Norm. For comparison, Fig. 9 also contains the number of physically simulated number of fingers. As evidenced by the large error in Fig. 7A and the simulated salt distribution in Fig. 8A, the unperturbed case (method A) does not adequately reproduce fingering observed in the laboratory experiment (Fig. 1B). The numerically simulated salt distribution after 40 min (Fig. 8A) only produces two large
79
fingers probably due to large concentration gradients at the edges of the solute source. It is therefore necessary to introduce one of the perturbation methods B to F to the numerical simulation. When an initial perturbation at the domain center is induced (method B), the lowest error is obtained with a perturbation magnitude of 10−1 (Fig. 7B). However, errors between the physical and the numerical experiment remain large for this trigger. Compared to the unperturbed sand (method A), the application of method B leads to the formation
(A)
(B)
(C)
(D)
(E)
(F)
Fig. 8. UHD case: Results of numerical simulations at t = 40 min for all perturbation mechanisms with best fit using (A) no perturbation, (B) XB = 10−1, (C) XC = 10−1, (D) XD = 10−1, (E) XE = 10−1, (F) σln(K) = 10−1.
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
20 18 16 14 12 10 8 6 4 2 0
(A)
20
30 40 time [min]
10
20
30 40 time [min]
20 18 16 14 12 10 8 6 4 2 0
50
50
20
20
50
30 40 time [min]
50
laboratory σln(K) =1E-1 σln(K) =1E-5
(F)
10
30 40 time [min] laboratory XD =1E-1 XD =1E-5
(D)
10
20 18 16 14 12 10 8 6 4 2 0
laboratory XB =1E-1 XB =1E-7
(B)
10
laboratory XE =1E-1 XE =1E-5
(E)
20 18 16 14 12 10 8 6 4 2 0
50
number of fingers
number of fingers
30 40 time [min] laboratory XC =1E-1 XC =1E-4
(C)
10
20 18 16 14 12 10 8 6 4 2 0
20
number of fingers
number of fingers
10
20 18 16 14 12 10 8 6 4 2 0
laboratory unperturbed number of fingers
number of fingers
80
20
30 40 time [min]
50
Fig. 9. UHD case: Number of dense plume fingers versus time for all perturbation methods applied, comprising the perturbation magnitudes that give the worst and the best fit respectively.
of more fingers after 40 min (Fig. 8B). However, the number of fingers produced by method B does not adequately reproduce results observed in the physical experiment (Fig. 1B). Perturbation methods C and D apply perturbations directly at the solute source on top of the unsaturated zone. For both, methods C and D, the smallest error is found for a trigger value
of 10−1 (Fig. 7C, D). However, it should be noted that the trigger value does not have a large influence on the error. Methods C and D result in the formation of more plume fingers in the central region of the plume front (Fig. 8C, D) in comparison to method A, but are not able to reproduce the finger pattern observed in the physical experiment. The amount of fingers
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
formed throughout the simulation for methods C and D is shown in Fig. 9C and D. Fig. 7E shows that a spatially and temporally random noise of solute concentration (method E) has a substantial impact on simulated plume fingering if the trigger magnitude is large. However, with increasing trigger magnitude errors substantially decrease (Fig. 7E). The best agreement between the physical and the numerical experiment is obtained with the largest trigger magnitude of 10−1. For this trigger, the numerically simulated salt distribution after 40 min is shown in Fig. 8E. The number of fingers produced by method E is larger than those of perturbation methods A to D, and it is almost identical to the number of fingers generated in the physical experiment (Figs. 1B and 8E). Perturbation method F (random K-field) gives the smallest error at a trigger magnitude of 10−1 (Fig. 7F). The associated plume fingering is shown in Fig. 8F. Similar to method E, use of a random K-field in method F realistically reproduces dense plume fingering under unsaturated conditions (Figs. 1B and 8F). Fig. 9F shows the number of dense plume fingers created with a random K-field. The trend in the amount of fingers produced in the physical experiment is in agreement with observations from the physical experiment. Both physical and numerical observations show an initially increasing and afterwards decreasing number of fingers. It should be noted that perturbation with a random K-field even for small variances (10−5, Fig. 9F) significantly improves the simulated
81
transport behavior and plume finger development compared to the unperturbed homogeneous reference case (Fig. 9A). 6. Summary and conclusion Perturbation mechanisms have been frequently used in recent literature to induce instabilities in the numerical representation of variable-density systems (Schincariol, 1998; Schincariol et al., 1997; Simmons et al., 1999; Weatherill et al., 2004; Xie et al., 2011, 2012). However, studies typically focused on a particular perturbation mechanism (e.g. Simmons et al. (1999), Weatherill et al. (2004), Xie et al. (2012)). Therefore, this study included various perturbation mechanisms documented in the literature of variable-density flow. Furthermore it has been extended by perturbation mechanisms that have not been used in the context of variable-density flow before. A quantitative investigation of the sensitivity of trigger magnitudes and a systematic evaluation of perturbation methods are still absent in the literature. Therefore, this study is aimed at systematically evaluating and comparing mechanisms to realistically create plume fingering. A laboratory experiment has been numerically re-simulated and different perturbation methods were applied. The quality of numerical results and hence quality of individual perturbation methods were evaluated by comparing simulated numerical results with laboratory results obtained by Simmons et al. (2002). Comparison was conducted by quantifying the number
(A)
(B)
(C)
(D)
Fig. 10. Results at t = 40 min for SHD (left column) and UHD (right column) physical (upper row) and numerical (bottom row) experiments (Simmons et al., 2002) for perturbation methods and trigger magnitudes giving the smallest error in the L2-Norm (Figs. 4, 7).
82
C.J.M. Cremer, T. Graf / Journal of Contaminant Hydrology 173 (2015) 69–82
of dense plume fingers formed during the experiment and by visual comparison of salt distribution in the domain. To the authors' knowledge, this is the first study that compares and evaluates multiple perturbation methods for the simulation of saturated and unsaturated density-driven flow and solute transport in porous media. The main new key findings of this study are: 1. Perturbation is necessary to realistically reproduce dense plume migration under saturated and unsaturated conditions. 2. While a spatially random, time-constant perturbation of the solute source (method C) is appropriate to simulate plume fingers under saturated conditions (Fig. 10), that method fails under unsaturated conditions. 3. Plume fingering under saturated–unsaturated flow conditions can be reproduced best with a spatially and temporally random noise of solute concentration (method E, Fig. 10) and with a random K-field (method F). It can be concluded that perturbation is needed in numerical models to realistically reproduce dense plume migration. While all of the evaluated perturbation methods have an impact on the migration of a dense plume, only some methods were able to give a good approximation of real plume migration simulated in the laboratory experiment. It can also be concluded that pore-scale heterogeneity (represented by methods E and F) is the main driving mechanism for dense plume fingering in homogeneous porous media. The findings of this study can serve as a basis for later benchmarking of density-driven flow in porous media under saturated and unsaturated conditions. Additionally, it constitutes a progress for the detailed realistic representation of dense plume fingering. Future work could focus on the applicability of perturbation mechanisms to different field-scale problems and the use of field data to assess the quality of various perturbation mechanisms. Acknowledgment We are grateful to C.T. Simmons (Flinders University, Australia) for providing photographic material and for providing input on the laboratory experiment, and to O. Cirpka (Universität Tübingen, Germany) for providing a random field generator and advice concerning its use. We thank the editorial board of the Journal of Contaminant Hydrology (B. Sleep) for handling our manuscript. Furthermore, we are grateful for the constructive comments by the two anonymous reviewers which helped to improve the manuscript. References Bear, J., Bachmat, Y., 1990. Introduction to Modeling of Transport Phenomena in Porous Media vol. 4. Springer. Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24 (5), 755–769.
Diersch, H.J.G., Kolditz, O., 2002. High-density flow and transport in porous media: approaches and challenges. Adv. Water Resour. 25, 899–944. Goswami, R.R., Clement, T.P., Hayworth, J.H., 2011. Comparison of numerical techniques used for simulating variable-density flow and transport experiments. J. Hydrol. Eng. 17 (2), 272–282. Horton, C.W., Rogers, F.T., 1945. Convection currents in a porous medium. J. Appl. Phys. 16, 367–370. Illangasekare, T., Tyler, S.W., Clement, T.P., Villholth, K.G., Perera, A.P.G.R.L., Obeysekera, J., Jensen, K., 2006. Impacts of the 2004 tsunami on groundwater resources in Sri Lanka. Water Resour. Res. 42 (5). Lapwood, E.R., 1948. Convection of a fluid in a porous medium. Proc. Camb. Philol. Soc. 44, 508–521. List, E.J., 1965. The Stability and Mixing of a Density-stratified Horizontal Flow in a Saturated Porous Medium. (Dissertation (Ph.D.)). California Institute of Technology. Scheidegger, A.E., 1954. Statistical hydrodynamics in porous media. J. Appl. Phys. 25, 994. Schincariol, R.A., 1998. Dispersive mixing dynamics of dense miscible plumes: natural perturbation initiation by local-scale heterogeneities. J. Contam. Hydrol. 34, 247–271. Schincariol, R.A., Schwartz, F.W., 1990. An experimental investigation of variable density flow and mixing in homogeneous and heterogeneous media. Water Resour. Res. 26 (10), 2317–2329. Schincariol, R.A., Schwartz, F.W., Mendoza, C.A., 1994. On the generation of instabilities in variable density flow. Water Resour. Res. 30 (4), 913–927. Schincariol, R.A., Schwartz, F.W., Mendoza, C.A., 1997. Instabilities in variable density flows: stability and sensitivity analyses for homogeneous and heterogeneous media. Water Resour. Res. 33 (1), 31–41. Simmons, C.T., 2005. Variable density groundwater flow: from current challenges to future possibilities. Hydrogeol. J. 13, 116–119. Simmons, C.T., Narayan, K.A., Wooding, R.A., 1999. On a test case for densitydependent groundwater flow and solute transport models: the salt lake problem. Water Resour. Res. 35 (12), 3607–3620. Simmons, C.T., Fenstemaker, T.R., Sharp, J.M., 2001. Variable-density groundwater flow and solute transport in heterogeneous porous media: approaches, resolutions and future challenges. J. Contam. Hydrol. 52, 245–275. Simmons, C.T., Pierini, M.L., Hutson, J.L., 2002. Laboratory investigation of variable-density flow and solute transport in unsaturated–saturated porous media. Transp. Porous Media 47, 215–244. Stöckl, L., Houben, G., 2012. Flow dynamics and age stratification of freshwater lenses: experiments and modeling. J. Hydrol. 458–459, 9–15. Therrien, R., McLaren, R.G., Sudicky, E.A., Panday, S.M., 2012. HydroGeoSphere — A Three-dimensional Numerical Model Describing Fully-integrated Subsurface and Surface Flow and Solute Transport. Groundwater Simulations Group, University of Waterloo, Canada. Van Dam, R.L., Simmons, C.T., Hyndman, D.W., Wood, W.W., 2009. Natural free convection in porous media: first field documentation in groundwater. Geophys. Res. Lett. 36 (11). Van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. Voss, C.I., Souza, W.R., 1987. Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater–saltwater transition zone. Water Resour. Res. 23, 1851–1866. Weatherill, D., Simmons, C.T., Voss, C.I., Robinson, N.I., 2004. Testing densitydependent groundwater models: two-dimensional steady state unstable convection in infinite, finite and inclined porous layers. Adv. Water Resour. 27, 547–562. Wooding, R.A., Tyler, S.W., White, I., Anderson, P.A., 1997a. Convection in groundwater below an evaporating salt lake: 1. Onset of instability. Water Resour. Res. 33 (6), 1199–1217. Wooding, R.A., Tyler, S.W., White, I., Anderson, P.A., 1997b. Convection in groundwater below an evaporating salt lake: 2. Evolution of fingers or plumes. Water Resour. Res. 33 (6), 1219–1228. Xie, Y., Simmons, C.T., Werner, A.D., 2011. Speed of free convective fingering in porous media. Water Resour. Res. 47 (11). Xie, Y., Simmons, C.T., Werner, A.D., Diersch, H.J.G., 2012. Prediction and uncertainty of free convection phenomena in porous media. Water Resour. Res. 48 (2).
