Lattice Boltzmann method with two relaxation times for advection–diffusion equation: Third order analysis and stability analysis

The objectives of this study are to investigate the third order accuracy and linear stability of the lattice Boltzmann method (LBM) with the two-relaxation-time collision operator (LTRT) for the advection–diffusion equation (ADE) and compare the LTRT model with the single-relaxation-time (LBGK) model. While the LBGK has been used extensively, the LTRT appears to be a more flexible model because it uses two relaxation times. The extra relaxation time can be used to improve solution accuracy and/or stability. This study conducts a third order Chapman–Enskog expansion on the LTRT to recover the macroscopic differential equations up to the third order. The dependency of third order terms on the relaxation times is obtained for different types of equilibrium distribution functions (EDFs) and lattices. By selecting proper relaxation times, the numerical dispersion can be significantly reduced. Furthermore, to improve solution accuracy, this study introduces pseudo-velocities to develop new EDFs to reduce the second order numerical diffusion. This study also derives stability domains based on the lattice Peclet number and Courant number for different types of lattices, EDFs and different values of relaxation times, while conducting linear stability analysis on the LTRT. Numerical examples demonstrate the improvement of the LTRT solution accuracy and stability by selecting proper relaxation times, lattice Peclet number and Courant number.     

Anisotropic dispersive Henry problem

The Henry problem has played a key role in our understanding of seawater intrusion into coastal aquifers and in benchmarking density dependent flow codes. This paper seeks to modify Henry’s problem to ensure sensitivity to density variations and vertical salinity profiles that resemble field observations. In the proposed problem, the “dispersive Henry problem”, mixing is represented by means of the traditional Scheidegger dispersion tensor (dispersivity times water flux). Anisotropy in the hydraulic conductivity is acknowledged and Henry’s seaside boundary condition of prescribed salt concentration is replaced by a flux dependent boundary condition, which represents more realistically salt transport across the seaside boundary. This problem turns out to be very sensitive to density variations and its solution gets closer to reality. However, an improvement in the traditional Henry problem (gain in sensitivity and realism) can be also achieved if the value of the Peclet number is significantly reduced.Although the dispersive problem lacks an analytical solution, it can shed light on flow in coastal aquifers. It provides significant information about the factors controlling seawater penetration, width of the mixing zone and influx of seawater. The width of the mixing zone depends basically on dispersion with longitudinal and transverse dispersion controlling different parts of the mixing zone but displaying similar overall effects. Toe penetration is mainly controlled by the horizontal permeability and by the geometric mean of the dispersivities. Finally, transverse dispersivity and the geometric mean of the hydraulic conductivity are the leading parameters controlling the amount of saltwater that enters the aquifer.     

Laplace-transform analytic element solution of transient flow in porous media

A Laplace-transform analytic element method (LT-AEM) is described for the solution of transient flow problems in porous media. Following Laplace transformation of the original flow problem, the analytic element method (AEM) is used to solve the resultant time-independent modified Helmholtz equation, and the solution is inverted numerically back into the time domain. The solution is entirely general, retaining the mathematical elegance and computational efficiency of the AEM while being amenable to parallel computation. It is especially well suited for problems in which a solution is required at a limited number of points in space–time, and for problems involving materials with sharply contrasting hydraulic properties. We illustrate the LT-AEM on transient flow through a uniform confined aquifer with a circular inclusion of contrasting hydraulic conductivity and specific storage. Our results compare well with published analytical solutions in the special case of radial flow.     

Eulerian spatial moments for solute transport in three-dimensional heterogeneous,dual-permeability media

A Eulerian analytical method is developed for nonreactive solute transport in heterogeneous, dual-permeability media where the hydraulic conductivities in fracture and matrix domains are both assumed to be stochastic processes. The analytical solution for the mean concentration is given explicitly in Fourier and Laplace transforms. Instead of using the fast fourier transform method to numerically invert the solution to real space (Hu et al., 2002), we apply the general relationship between spatial moments and concentration (Naff, 1990; Hu et al., 1997) to obtain the analytical solutions for the spatial moments up to the second for a pulse input of the solute. Owing to its accuracy and efficiency, the analytical method can be used to check the semi-analytical and Monte Carlo numerical methods before they are applied to more complicated studies. The analytical method can be also used during screening studies to identify the most significant transport parameters for further analysis. In this study, the analytical results have been compared with those obtained from the semi-analytical method (Hu et al., 2002) and the comparison shows that the semi-analytical method is robust. It is clearly shown from the analytical solution that the three factors, local dispersion, conductivity variation in each domain and velocity convection flow difference in the two domains, play different roles on the solute plume spreading in longitudinal and transverse directions. The calculation results also indicate that when the log-conductivity variance in matrix is 10 times less than its counterpart in fractures, it will hardly influence the solute transport, whether the conductivity field is matrix is treated as a homogeneous or random field.     

Two-region non-Darcian flow toward a well in a confined aquifer

We have derived an analytical solution for two-region flow toward a well in a confined aquifer based on a linearization method. The two-region flow includes Izbash non-Darcian flow near the well and Darcian flow in the rest of the aquifer. The wellbore storage is also considered. The type curves in the non-Darcian and Darcian flow domains are obtained by a numerical Laplace inversion method incorporated in MATLAB programs. We have compared our results with the one-region Darcian flow model (Theis). Our solutions agree with those of Sen [Sen Z. Type curves for two-region well flow. J Hydr Eng 1988;114(12):1461–84] which were obtained using the Boltzmann transform at late times for fully turbulent flow, while some difference has been found at early and moderate times. We have defined a dimensionless non-Darcian hydraulic conductivity term which is shown to be a key parameter for analyzing the two-region flow. A smaller dimensionless non-Darcian hydraulic conductivity results in a larger drawdown in the non-Darcian flow region at late times. However, the dimensionless non-Darcian hydraulic conductivity does not affect the slope of the dimensionless drawdown versus the logarithmic dimensionless time in the non-Darcian flow region at late times. The dimensionless non-Darcian hydraulic conductivity does not affect the late time drawdown in the Darcian flow region.     

A closed form slug test theory for high permeability aquifers

We incorporate a linear estimate of casing friction into the analytical slug test theory of Springer and Gelhar (1991) for high permeability aquifers. The modified theory elucidates the influence of inertia and casing friction on consistent, closed form equations for the free surface, pressure, and velocity fluctuations for overdamped and underdamped conditions. A consistent, but small, correction for kinetic energy is included as well. A characteristic velocity linearizes the turbulent casing shear stress so that an analytical solution for attenuated, phase shifted pressure fluctuations fits a single parameter (damping frequency) to transducer data from any depth in the casing. Underdamped slug tests of 0.3, 0.6, and 1 m amplitudes at five transducer depths in a 5.1 cm diameter PVC well 21 m deep in the Plymouth-Carver Aquifer yield a consistent hydraulic conductivity of 1.5 x 10(-3) m/s. The Springer and Gelhar (1991) model underestimates the hydraulic conductivity for these tests by as much as 25% by improperly ascribing smooth turbulent casing friction to the aquifer. The match point normalization of Butler (1998) agrees with our fitted hydraulic conductivity, however, when friction is included in the damping frequency. Zurbuchen et al. (2002) use a numerical model to establish a similar sensitivity of hydraulic conductivity to nonlinear casing friction.     

Saltwater Upconing Due to Cyclic Pumping by Horizontal Wells in Freshwater Lenses

This article deals with the quantification of saltwater upconing below horizontal wells in freshwater lenses using analytical solutions as a computationally fast alternative to numerical simulations. Comparisons between analytical calculations and numerical simulations are presented regarding three aspects: (1) cyclic pumping; (2) dispersion; and (3) finite horizontal wells in a finite domain (a freshwater lens). Various hydrogeological conditions and pumping regimes within a dry half year are considered. The results show that the influence of elastic and phreatic storage (which are not taken into account in the analytical solutions) on the upconing of the interface is minimal. Furthermore, the analytical calculations based on the interface approach compare well with numerical simulations as long as the dimensionless interface upconing is below 1/3, which is in line with previous studies on steady pumping. Superimposing an analytical solution for mixing by dispersion below the well over an analytical solution based on the interface approach is appropriate in case the vertical flow velocity around the interface is nearly constant but should not be used for estimating the salinity of the pumped groundwater. The analytical calculations of interface upconing below a finite horizontal well compare well with the numerical simulations in case the distance between the horizontal well and the initial interface does not vary significantly along the well and in case the natural fluctuation of the freshwater lens is small. In order to maintain a low level of salinity in the well during a dry half year, the dimensionless analytically calculated interface upconing should stay below 0.25.     

Flow in unsaturated random porous media,nonlinear numerical analysis and comparison to analytical stochastic models

This work presents a rigorous numerical validation of analytical stochastic models of steady state unsaturated flow in heterogeneous porous media. It also provides a crucial link between stochastic theory based on simplifying assumptions and empirical field and simulation evidence of variably saturated flow in actual or realistic hypothetical heterogeneous porous media. Statistical properties of unsaturated hydraulic conductivity, soil water tension, and soil water flux in heterogeneous soils are investigated through high resolution Monte Carlo simulations of a wide range of steady state flow problems in a quasi-unbounded domain. In agreement with assumptions in analytical stochastic models of unsaturated flow, hydraulic conductivity and soil water tension are found to be lognormally and normally distributed, respectively. In contrast, simulations indicate that in moderate to strong variable conductivity fields, longitudinal flux is highly skewed. Transverse flux distributions are leptokurtic. the moments of the probability distributions obtained from Monte Carlo simulations are compared to modified first-order analytical models. Under moderate to strong heterogeneous soil flux conditions (σ²_y≥1), analytical solutions overestimate variability in soil water tension by up to 40% as soil heterogeneity increases, and underestimate variability of both flux components by up to a factor 5. Theoretically predicted model (cross-)covariance agree well with the numerical sample (cross-)covarianaces. Statistical moments are shown to be consistent with observed physical characteristics of unsaturated flow in heterogeneous soils.©1998 Elsevier Science Limited. All rights reserved     

Composite analytical solutions for a soil vapour extraction system

Soil vapour extraction (SVE) is a common remediation technique for cleaning up unsaturated soils contaminated by volatile organic compounds (VOCs). Analytical solutions, which result from simple mathematical models, can allow the fast approximation of the time‐dependent effluent concentration and the gaining of insight into the processes that take place during soil remediation. Deriving the analytical solutions to advection–dispersion equations that simultaneously take into account the mechanical dispersion and molecular diffusion is very difficult because of the variable dependence of governing equations' coefficients. In this study, we first present two simplified analytical solutions that only consider mechanical dispersion or molecular diffusion. The two developed analytical solutions are compared with the numerical solution that simultaneously considers both mechanical dispersion and molecular diffusion to examine the applicability of the two simplified analytical solutions and distinguishes the individual contribution of the mechanical dispersion and molecular diffusion to total VOCs transport in an SVE system. Results show that dispersion plays an important role during SVE decontamination and neither the diffusion‐dominated solution nor the dispersion‐dominated solution can agree well with the numerical solution when both mechanical dispersion and molecular diffusion have significant contributions to the total VOCs transport flux. A composite analytical solution that linearly couples the diffusion‐ and dispersion‐dominated analytical solutions, which is proposed herein to eliminate the discrepancy between the analytical solutions and the numerical solution. Results indicate that the proposed composite analytical solution agrees well with the numerical solution and is an effective tool for quickly and accurately evaluating the time‐dependent effluent concentration for parameters of the different ranges of interest in an SVE remedial system. Copyright © 2006 John Wiley & Sons, Ltd.     

Coordinate mapping of analytical contaminant transport solutions to non-uniform flow fields

Existing analytical solutions to 2D and 3D contaminant transport problems are limited by the mathematically convenient assumption of uniform flow. An approximate method is developed herein for coordinate mapping of 2D (vertically-averaged) transport solutions to non-uniform steady-state irrotational and divergence-free flow fields in single-layer aquifers. The method enables existing analytical transport solutions to be applied to aquifer systems with wells, non-uniform saturated thickness, surface water features, and (to a limited degree) heterogeneous hydraulic conductivity and recharge. This mass-conservative coordinate mapping approach is inexact in its approximation of the dispersion process but is still sufficiently accurate for many simple flow systems. The degree of model error is directly proportional to the variation of velocity magnitude within the domain. These mapped analytical solutions are compared to numerical simulation results and the coordinate mapping errors are investigated. The methods described herein may be used in the traditional capacity of analytical transport models, i.e., screening and preliminary site assessment, without sacrificing accuracy by assuming locally uniform flow conditions or applying an ad-hoc coordinate transformation. The solutions benefit from the traditional advantages of analytical methods, particularly the removal of artifacts due to spatial and temporal discretization: no time-stepping or numerical discretization is required.     

Two-dimensional power series solution for non-axisymmetrical transport in a radially convergent tracer test with scale-dependent dispersion

It has been known for many years that dispersivity increases with solute travel distance in a subsurface environment. The increase of dispersivity with solute travel distance results from the significant variation of hydraulic properties of heterogeneous media and was identified in the literature as scale-dependent dispersion. This study presents an analytical solution for describing two-dimensional non-axisymmetrical solute transport in a radially convergent flow tracer test with scale-dependent dispersion. The power series technique coupling with the Laplace and finite Fourier cosine transform has been applied to yield the analytical solution to the two-dimensional, scale-dependent advection–dispersion equation in cylindrical coordinates with variable-dependent coefficients. Comparison between the breakthrough curves of the power series solution and the numerical solutions shows excellent agreement at different observation points and for various ranges of scale-related transport parameters of interest. The developed power series solution facilitates fast prediction of the breakthrough curves at any observation point.     

An effective scale-dependent dispersivity deduced from a purely convective flow field

Abstract

In the case of straight flow but with hydraulic conductivity varying in a transverse direction, the distribution of hydraulic conductivity has been determined for which the breakthrough curve due to convection only will have the same analytical form as the onedimensional convection/dispersion equation solution at the outlet end of a porous medium. That distribution is found exactly and it is very similar to the lognormal distribution. This result is significant since field evidence indicates that the logarithm of hydraulic conductivity is normally distributed. For the case considered, a simple relation between dispersivity and the coefficient of variation of hydraulic conductivity is found. One can thus determine very simply dispersivity in terms of the parameters of the distribution of hydraulic conductivity. This is particularly useful to estimate dispersivity in various cells of finite difference or finite element models when the distribution of hydraulic conductivity is not stationary, i.e. varies in space.     

Vulnerability of three spatially variable soil groups for solute leaching

Solute leaching in unsaturated soil is influenced by the variability in hydraulic functions (water retention and conductivity) that govern the flow process. Variability in measured soil hydraulic functions of a coarse-, medium- and fine-textured soil group was quantified with the scaling theory of similar media. Solute leaching in these soils was calculated with Monte Carlo simulation assuming, successively, hydraulic conductivity, K, volumetric water content, 0, and pressure head, h, to be constant. In addition to variability in hydraulic functions, variability in the solute retardation factor was also taken into account. To examine this effect five solutes were considered: a conservative solute (chloride), a non-retarded solute subject to decay (nitrate), a retarded solute that does not decay (cadmium) and two organic solutes which are retarded but have different sorption and decay parameters (the pesticide atrazine and a chlorinated hydrocarbon). The numerical results obtained with Monte Carlo simulation were in a number of instances verified with analytical solutions. The three soil groups distinguished showed considerable differences in vulnerability for leaching of the five solutes, emphasizing the importance of the effect of variability in soil hydraulic functions when studying solute leaching. Numerical and analytical results showed good agreement. Therefore, in relatively simple situations analytical solutions are attractive. However, in complicated situations, analytical solutions are cumbersome and numerical solutions are the only realistic alternative.     

Numerical analysis of one-dimensional unsaturated flow in layered soils

The paper deals with numerical solutions to the Richards equation to simulate one-dimensional flow processes in the unsaturated zone of layered soil profiles. The equation is expressed in the pressure-based form and a finite-difference algorithm is developed for accurately estimating the values of the hydraulic conductivity between two neighboring nodes positioned in different soil layers, often referred to as the interlayer hydraulic conductivity. The algorithm is based upon flux conservation and continuity of pressure potential at the interface between two consecutive layers, and does not add significantly to simulation run time. The validity of the model is established for a number of test problems by comparing numerical results with the analytical solutions developed by Srivastava and Yeh²⁹ which hold for vertical infiltration towards the water table in a two-layer soil profile. The results show a significant reduction in relative mass balance errors when using the proposed model. Some specific insights into its numerical performance are also gained by comparisons with a numerical model in which the more common geometric averaging operator acts on the interlayer conductivities.     

Solutions for Non‐Darcian Flow to an Extended Well in Fractured Rock

We have investigated non‐Darcian flow to a vertical fracture represented as an extended well using a linearization procedure and a finite difference method in this study. Approximate analytical solutions have been obtained with and without the consideration of fracture storage based on the linearization procedure. A numerical solution for such a non‐Darcian flow case has also been obtained with a finite difference method. We have compared the numerical solution with the approximate analytical solutions obtained by the linearization method and the Boltzmann transform. The results indicate that the linearized solution agrees generally well with the numerical solution at late times, and underestimates the dimensionless drawdown at early times, no matter if the fracture storage is considered or not. When the fracture storage is excluded, the Boltzmann transform solution overestimates the dimensionless drawdown during the entire pumping period. The dimensionless drawdowns in the fracture with fracture storage for different values of dimensionless non‐Darcian hydraulic conductivity β approach the same asymptotic value at early times. A larger β value results in a smaller dimensionless drawdown in both the fracture and the aquifer when the fracture storage is included. The dimensionless drawdown is approximately proportional to the square root of the dimensionless time at late times.     

Modelling of the groundwater flow and of tracer movement in the porous and fissured media: Chalk Aquifer (Northern part of Paris Basin,France)

A groundwater flow model has been developed in order to study the chalk aquifer of Paris Basin, based on most of the geological and hydrological available data. The numerical processes are intended to modelling the groundwater flow in the Senonian (Late Cretaceous) formations and to visualize the tracer movement in groundwater resources in the experimental site of LaSalle Beauvais (northern part Paris Basin). Both objectives were achieved as follows: (i) the comprehension of the spatial distribution of the hydraulic conductivity in the chalk aquifer taking into account the characteristics of the hydrogeological system and (ii) the use of the analytical solution for describing one‐dimensional to two‐dimensional solute transport in a unidirectional steady‐state flow tracer with scale‐dependent dispersion. Advection and diffusion mechanisms are taken into account. Comparison between the breakthrough curves of the analytical and the numerical solutions provided an excellent agreement for various ranges of scale‐related transport parameters of interest. The developed power series solution facilitates fast prediction of the breakthrough curves at each observation point. Thus, the derived new solutions are widely applicable and are very useful for the validation of numerical transport. The numerical approach is carried out by MT3DMS, a Modular 3‐D Multi‐Species Transport Model for Simulation of Advection, Dispersion, and Chemical Reactions of Contaminants in Groundwater Systems, and based on total variation‐diminishing method using the ULTIMATE algorithm. The estimation of the infected surface could constitute an approach in water management and allows to prevent the risks of pollution and to manage the groundwater resource from a durable development perspective. Copyright © 2016 John Wiley & Sons, Ltd.     

A comparison of solute-transport solution techniques and their effect on sensitivity analysis and inverse modeling results

Five common numerical techniques for solving the advection-dispersion equation (finite difference, predictor corrector, total variation diminishing, method of characteristics, and modified method of characteristics) were tested using simulations of a controlled conservative tracer-test experiment through a heterogeneous, two-dimensional sand tank. The experimental facility was constructed using discrete, randomly distributed, homogeneous blocks of five sand types. This experimental model provides an opportunity to compare the solution techniques: the heterogeneous hydraulic-conductivity distribution of known structure can be accurately represented by a numerical model, and detailed measurements can be compared with simulated concentrations and total flow through the tank. The present work uses this opportunity to investigate how three common types of results--simulated breakthrough curves, sensitivity analysis, and calibrated parameter values--change in this heterogeneous situation given the different methods of simulating solute transport. The breakthrough curves show that simulated peak concentrations, even at very fine grid spacings, varied between the techniques because of different amounts of numerical dispersion. Sensitivity-analysis results revealed: (1) a high correlation between hydraulic conductivity and porosity given the concentration and flow observations used, so that both could not be estimated; and (2) that the breakthrough curve data did not provide enough information to estimate individual values of dispersivity for the five sands. This study demonstrates that the choice of assigned dispersivity and the amount of numerical dispersion present in the solution technique influence estimated hydraulic conductivity values to a surprising degree.     

Sea‐level rise impact on fresh groundwater lenses in two‐layer small islands

The fresh groundwater lenses (FGLs) of small islands can be highly vulnerable to climate change impacts, including sea‐level rise (SLR). Many real cases of atoll or sandy islands involve two‐layer hydrogeological conceptualizations. In this paper, the influential factors that affect FGLs in two‐layer small islands subject to SLR are investigated. An analytical solution describing FGLs in circular islands, composed of two geological layers, is developed for the simplified case of steady‐state and sharp‐interface conditions. An application of the developed model is demonstrated to estimate the FGL thickness of some real‐world islands by comparison with existing FGL thickness data. Furthermore, numerical modelling is applied to extend the analysis to consider dispersion effects and to confirm comparable results for both cases. Sensitivity analyses are used to assess the importance of land‐surface inundation caused by SLR, relative to other parameters (i.e. thickness of aquifer layers, hydraulic conductivity, recharge rate and land‐surface slope) that influence the FGL. Dimensionless parameters are used to generalize the findings. The results demonstrate that land‐surface inundation has a considerable impact on a FGL influenced by SLR, as expected, although the FGL volume is more sensitive to recharge, aquifer thickness and hydraulic conductivity than SLR impacts, considering typical parameter ranges. The methodology presented in this study provides water resource managers with a rapid‐assessment tool for evaluating the likely impacts of SLR and accompanying LSI on FGLs.     

Aquifer Properties Determined from Two Analytical Solutions

In the analysis of pumping test data, the quality of the determined aquifer parameters can be greatly improved by using a proper model of the aquifer system. Moench (1995) provided an analytical solution for flow to a well partially penetrating an unconfined aquifer. His solution, in contrast to the Neuman solution (1974), accounts for the noninstantaneous decline of the water table (delayed yield). Consequently, the calculated drawdown in these two solutions is different under certain circumstances, and this difference may therefore affect the computation of aquifer properties from pumping test data. This paper uses an inverse computational method to calculate four aquifer parameters as well as a delayed yield parameter, α₁ from pumping test data using both the Neuman (1974) and Moench (1995) solutions. Time-drawdown data sets from a pumping test in an unconfined alluvial aquifer near Grand Island, Nebraska, were analyzed. In single-well analyses, horizontal hydraulic conductivity values derived from the Moench solution are lower, but vertical hydraulic conductivity values are higher than those calculated from the Neuman solution. However, the hydraulic conductivity values in composite-well analyses from both solutions become very close. Furthermore, the Neuman solution produces similar hydraulic conductivity values in the single-well and composite-well analyses, but the Moench solution does not. While variable α₁, seems to play a role in affecting the computation of aquifer parameters in the single-well analysis, a much smaller effect was observed in the composite-well analysis. In general, specific yield determined using the Moench solution could be slightly higher than the values from the Neuman solution; however, they are still lower than the realistic values for sand and gravel aquifers.     

Uncertainty Estimation in One‐Dimensional Heat Transport Model for Heterogeneous Porous Medium

In many practical applications, the rates for ground water recharge and discharge are determined based on the analytical solution developed by Bredehoeft and Papadopulos (1965) to the one‐dimensional steady‐state heat transport equation. Groundwater flow processes are affected by the heterogeneity of subsurface systems; yet, the details of which cannot be anticipated precisely. There exists a great deal of uncertainty (variability) associated with the application of Bredehoeft and Papadopulos' solution (1965) to the field‐scale heat transport problems. However, the quantification of uncertainty involved in such application has so far not been addressed, which is the objective of this wok. In addition, the influence of the statistical properties of log hydraulic conductivity field on the variability in temperature field in a heterogeneous aquifer is also investigated. The results of the analysis demonstrate that the variability (or uncertainty) in the temperature field increases with the correlation scale of the log hydraulic conductivity covariance function and the variability of temperature field also depends positively on the position.