Quantifying geological uncertainty for flow and transport modeling in multi-modal heterogeneous formations

In this work, we address the problem of characterizing the heterogeneity and uncertainty of hydraulic properties for complex geological settings. Hereby, we distinguish between two scales of heterogeneity, namely the hydrofacies structure and the intrafacies variability of the hydraulic properties. We employ multiple-point geostatistics to characterize the hydrofacies architecture. The multiple-point statistics are borrowed from a training image that is designed to reflect the prior geological conceptualization. The intrafacies variability of the hydraulic properties is represented using conventional two-point correlation methods, more precisely, spatial covariance models under a multi-Gaussian spatial law. We address the different levels and sources of uncertainty in characterizing the subsurface heterogeneity, and explore their effect on groundwater flow and transport predictions. Typically, uncertainty is assessed by way of many images, termed realizations, of a fixed statistical model. However, in many cases, sampling from a fixed stochastic model does not adequately represent the space of uncertainty. It neglects the uncertainty related to the selection of the stochastic model and the estimation of its input parameters. We acknowledge the uncertainty inherent in the definition of the prior conceptual model of aquifer architecture and in the estimation of global statistics, anisotropy, and correlation scales. Spatial bootstrap is used to assess the uncertainty of the unknown statistical parameters. As an illustrative example, we employ a synthetic field that represents a fluvial setting consisting of an interconnected network of channel sands embedded within finer-grained floodplain material. For this highly non-stationary setting we quantify the groundwater flow and transport model prediction uncertainty for various levels of hydrogeological uncertainty. Results indicate the importance of accurately describing the facies geometry, especially for transport predictions.     

Modeling groundwater velocity uncertainty in nonstationary composite porous media

Despite the intensive research over the past decades in the field of stochastic subsurface hydrology, our ability to analyze and model heterogeneous groundwater systems remains limited. Most existing theories are either too restrictive to handle practical complexity or too expensive to be applied to realistic problem sizes. In this paper we present approximate, closed-form equations that allow modeling 2D nonstationary flows in statistically inhomogeneous aquifers, including composite aquifers containing multiple zones characterized by different statistical models. The composite representation has the effect of decreasing the variance of deviations from the mean, relaxing the limitation of the small-perturbation assumption. The simple formulas are illustrated with a number of examples and compared with a corresponding first-order nonstationary numerical analysis and Monte Carlo simulation. The results show that, despite the gross simplifications, the closed-form equations are robust and able to capture complex variance dynamics, reproducing surprisingly well the first-order numerical solutions and the Monte Carlo simulation even in highly nonstationary, variable situations.     

Impact of log-transmissivity variogram structure on groundwater flow and transport predictions

We analyze the impact of the choice of the variogram model adopted to characterize the spatial variability of natural log-transmissivity on the evaluation of leading (statistical) moments of hydraulic heads and contaminant travel times and trajectories within mildly (randomly) heterogeneous two-dimensional porous systems. The study is motivated by the fact that in several practical situations the differences between various variogram types and a typical noisy sample variogram are small enough to suggest that one would often have a hard time deciding which of the tested models provides the best fit. Likewise, choosing amongst a set of seemingly likely variogram models estimated by means of geostatistical inverse models of flow equations can be difficult due to lack of sensitivity of available model discrimination criteria. We tackle the problem within the framework of numerical Monte Carlo simulations for mean uniform and radial flow scenarios. The effect of three commonly used isotropic variogram models, i.e., Gaussian, Exponential and Spherical, is analyzed. Our analysis clearly shows that (ensemble) mean values of the quantities of interest are not considerably influenced by the variogram shape for the range of parameters examined. Contrariwise, prediction variances of the quantities examined are significantly affected by the choice of the variogram model of the log-transmissivity field. The spatial distribution of the largest/lowest values of the relative differences observed amongst the tested models depends on a combination of variogram shape and parameters and relative distance from internal sources and the outer domain boundary. Our findings suggest the need of developing robust techniques to discriminate amongst a set of seemingly equally likely alternative variogram models in order to provide reliable uncertainty estimates of state variables.     

Statistical characteristics of flow as indicators of channeling in heterogeneous porous and fractured media

We introduce two new channeling indicators D_ic and D_cc based on the Lagrangian distribution of flow rates. On the basis of the participation ratio, these indicators characterize the extremes of both the flow-tube width distribution and the flow rate variation along flow lines. The participation ratio is an indicator biased toward the larger values of a distribution and is equal to the normalized ratio of the square of the first-order moment to the second-order moment. Compared with other existing indicators, they advantageously provide additional information on the flow channel geometry, are consistently applicable to both porous and fractured media, and are generally less variable for media generated using the same parameters than other indicators. Based on their computation for a broad range of porous and fracture permeability fields, we show that they consistently characterize two different geometric properties of channels. D_ic gives a characteristic scale of low-flow zones in porous media and a characteristic distance between effectively flowing structures in fractured cases. D_cc gives a characteristic scale of the extension of high-flow zones in porous media and a characteristic channel length in fractured media. D_ic is mostly determined by channel density and permeability variability. D_cc is, however, more affected by the nature of the correlation structure like the presence of permeability channels or fractures in porous media and the length distribution in fracture networks.     

Stochastic analysis of spatial variability in unconfined groundwater flow

In this paper, spatial variability in steady one-dimensional unconfined groundwater flow in heterogeneous formations is investigated. An approach to deriving the variance of the hydraulic head is developed using the nonlinear filter theory. The nonlinear governing equation describing the one-dimensional unconfined groundwater flow is decomposed into three linear partial differential equations using the perturbation method. The linear and quadratic frequency response functions are obtained from the first- and second-order perturbation equations using the spectral method. Furthermore, under the assumption of the exponential covariance function of log hydraulic conductivity, the analytical solutions of both the spectrum and the variance of the hydraulic head produced from the linear system are derived. The results show that the variance derived herein is less than that of Gelhar (1977). The reason is that the log transmissivity is linearized in Gelhars work. In addition, the analytical solutions of both the spectrum and the variance of the hydraulic head produced from the quadratic system are derived as well. It is found that the correlation scale and the trend in mean of log hydraulic conductivity are important to the dimensionless variance ratio.     

Coupling multiscale finite element method for consolidation analysis of heterogeneous saturated porous media

The multiscale finite element method is developed for solving the coupling problems of consolidation of heterogeneous saturated porous media under external loading conditions. Two sets of multiscale base functions are constructed, respectively, for the pressure field of fluid flow and the displacement field of solid skeleton. The coupling problems are then solved with a multiscale numerical procedure in space and time domain. The heterogeneities induced by permeabilities and mechanical parameters of the saturated porous media are both taken into account. Numerical experiments are carried out for different cases in comparison with the standard finite element method. The numerical results show that the coupling multiscale finite element method can be successfully used for solving the complicated coupling problems. It reduces greatly the computing effort in both memory and time for transient problems.     

Knowledge,transparency, and refutability in groundwater models,an example from the Death Valley regional groundwater flow system

This work demonstrates how available knowledge can be used to build more transparent and refutable computer models of groundwater systems. The Death Valley regional groundwater flow system, which surrounds a proposed site for a high level nuclear waste repository of the United States of America, and the Nevada National Security Site (NNSS), where nuclear weapons were tested, is used to explore model adequacy, identify parameters important to (and informed by) observations, and identify existing old and potential new observations important to predictions. Model development is pursued using a set of fundamental questions addressed with carefully designed metrics. Critical methods include using a hydrogeologic model, managing model nonlinearity by designing models that are robust while maintaining realism, using error-based weighting to combine disparate types of data, and identifying important and unimportant parameters and observations and optimizing parameter values with computationally frugal schemes. The frugal schemes employed in this study require relatively few (10–1000 s), parallelizable model runs. This is beneficial because models able to approximate the complex site geology defensibly tend to have high computational cost. The issue of model defensibility is particularly important given the contentious political issues involved.     

Density-dependent dispersion in heterogeneous porous media Part II: Comparison with nonlinear models

The results of a series of high-resolution numerical experiments are used to test and compare three nonlinear models for high-concentration-gradient dispersion. Gravity stable miscible displacement is considered. The first model, introduced by Hassanizadeh, is a modification of Fick’s law which involves a second-order term in the dispersive flux equation and an additional dispersion parameter β. The numerical experiments confirm the dependency of β on the flow rate. In addition, a dependency on travelled distance is observed. The model can successfully be applied to nearly homogeneous media (σ² = 0.1), but additional fitting is required for more heterogeneous media.The second and third models are based on homogenization of the local scale equations describing density-dependent transport. Egorov considers media that are heterogeneous on the Darcy scale, whereas Demidov starts at the pore-scale level. Both approaches result in a macroscopic balance equation in which the dispersion coefficient is a function of the dimensionless density gradient. In addition, an expression for the concentration variance is derived. For small σ², Egorov’s model predictions are in satisfactory agreement with the numerical experiments without the introduction of any new parameters. Demidov’s model involves an additional fitting parameter, but can be applied to more heterogeneous media as well.     

Flow and transport through two-dimensional isotropic media of binary conductivity distribution. Part 1: NUMERICAL methodology and semi-analytical solutions

Flow and transport take place in a heterogeneous medium made up from inclusions of conductivity K submerged in a matrix of conductivity K ₀. We consider two-dimensional isotropic media, with circular inclusions of uniform radii, that are placed at random and without overlap in the matrix. The system is completely characterized by the conductivity contrast =K/K ₀ and by the volume fraction n. The flow is uniform in the mean, of velocity U=const. The derivation of the velocity field is achieved by a numerical method of high accuracy, based on analytical elements. Approximate analytical solutions are derived by a few methods: composite elements, effective medium, dilute systems and first-order approximation in logconductivity variance. The latter was employed by Rubin (1995), while the dilute system approximation was used by Eames and Bush (1999) and Dagan and Lessoff (2001). Transport is solved in a Lagrangean framework, with trajectories determined numerically from the velocity field, by particle tracking. Results for the velocity variance and for the longitudinal macrodispersivity, for a few values of and n, are presented in Part 2.     

The stability of density-driven flows in saturated heterogeneous porous media

This work concludes the investigations into the stability of haline flows in saturated porous media. In the first part [33] a stability criterion for density-driven flow in a saturated homogeneous medium was derived excluding dispersion. In the second part [34], the effects of dispersion were included. The latter criterion made reasonable predictions of the stability regimes (indicated by the number of fingers present) as a function of density and dispersivity variations. We found out that destabilising variables caused an increase in the number of fingers and vice versa. The investigation is extended here for the effects of the medium heterogeneity. The cell problem derived via homogenization theory [20] is solved and its solution used to evaluate the elements of the macrodispersion tensor as functions of time for flow aligned parallel to gravity. The longitudinal coefficient exhibits asymptotic behaviour for favourable and moderately unfavourable density contrasts while it grows indefinitely for higher density contrasts. The range of densities stabilised by medium heterogeneities can thus be estimated from the behaviour of the coefficient. The d³f software program is used for the numerical simulations. The code uses the cell-centred finite volume and the implicit Euler techniques for the spatial and temporal discretisations respectively.     

Flow and transport through two-dimensional isotropic media of binary conductivity distribution. Part 2: NUMERICAL simulations and comparison with theoretical results

In the present part the results of numerical simulations of flow and transport in media made up from circular inclusions of conductivity K that are submerged in a matrix of conductivity K ₀, subjected to uniform mean velocity, are presented. This is achieved for a few values of =K/K ₀ (0.01, 0.1 and 10) and of the volume fraction n (0.05, 0.1 and 0.2). The numerical simulations (NS) are compared with the analytical approximate models presented in Part 1: the composite elements (CEA), the effective medium (EMA), the dilute system (DSA) and the first-order in the logconductivity variance (FOA). The comparison is made for the longitudinal velocity variance and for the longitudinal macrodispersivity. This is carried out for n<0.2, for which the theoretical and simulation models represent the same structure of random and independent inclusions distribution. The main result is that transport is quite accurately modeled by the EMA and CEA for low , for which _L is large, whereas in the case of =10, the EMA matches the NS for n<0.1. The first-order approximation is quite far apart from the NS for the values of examined . This material is based upon work supported by the National Science Foundation under Grant No. 0218914. Authors also wish to thank the Center of Computational Research, University at Buffalo for assistance in running numerical simulations.     

Numerical modeling of two-phase flow in heterogeneous permeable media with different capillarity pressures

Contrast in capillary pressure of heterogeneous permeable media can have a significant effect on the flow path in two-phase immiscible flow. Very little work has appeared on the subject of capillary heterogeneity despite the fact that in certain cases it may be as important as permeability heterogeneity. The discontinuity in saturation as a result of capillary continuity, and in some cases capillary discontinuity may arise from contrast in capillary pressure functions in heterogeneous permeable media leading to complications in numerical modeling. There are also other challenges for accurate numerical modeling due to distorted unstructured grids because of the grid orientation and numerical dispersion effects. Limited attempts have been made in the literature to assess the accuracy of fluid flow modeling in heterogeneous permeable media with capillarity heterogeneity. The basic mixed finite element (MFE) framework is a superior method for accurate flux calculation in heterogeneous media in comparison to the conventional finite difference and finite volume approaches. However, a deficiency in the MFE from the direct use of fractional flow formulation has been recognized lately in application to flow in permeable media with capillary heterogeneity. In this work, we propose a new consistent formulation in 3D in which the total velocity is expressed in terms of the wetting-phase potential gradient and the capillary potential gradient. In our formulation, the coefficient of the wetting potential gradient is in terms of the total mobility which is smoother than the wetting mobility. We combine the MFE and discontinuous Galerkin (DG) methods to solve the pressure equation and the saturation equation, respectively. Our numerical model is verified with 1D analytical solutions in homogeneous and heterogeneous media. We also present 2D examples to demonstrate the significance of capillary heterogeneity in flow, and a 3D example to demonstrate the negligible effect of distorted meshes on the numerical solution in our proposed algorithm.     

An efficient,high-order probabilistic collocation method on sparse grids for three-dimensional flow and solute transport in randomly heterogeneous porous media

In this study, a probabilistic collocation method (PCM) on sparse grids is used to solve stochastic equations describing flow and transport in three-dimensional, saturated, randomly heterogeneous porous media. The Karhunen–Loève decomposition is used to represent log hydraulic conductivity Y=ln_Ks$Y=\ln K_s$. The hydraulic head h   and average pore-velocity v$\mathbf{v}$ are obtained by solving the continuity equation coupled with Darcy’s law with random hydraulic conductivity field. The concentration is computed by solving a stochastic advection–dispersion equation with stochastic average pore-velocity v$\mathbf{v}$ computed from Darcy’s law. The PCM approach is an extension of the generalized polynomial chaos (gPC) that couples gPC with probabilistic collocation. By using sparse grid points in sample space rather than standard grids based on full tensor products, the PCM approach becomes much more efficient when applied to random processes with a large number of random dimensions. Monte Carlo (MC) simulations have also been conducted to verify accuracy of the PCM approach and to demonstrate that the PCM approach is computationally more efficient than MC simulations. The numerical examples demonstrate that the PCM approach on sparse grids can efficiently simulate solute transport in randomly heterogeneous porous media with large variances.     

A note on benchmarking of numerical models for density dependent flow in porous media

Verification of numerical models for density dependent flow in porous media (DDFPM) by the means of appropriate benchmark problems is a very important step in developing and using these models. Recently, Infinite Horizontal Box (IHB) problem was suggested as a possible benchmark problem for verification of DDFPM codes. IHB is based on Horton–Rogers–Lapwood (HRL) problem. Suitability of this problem for the benchmarking purpose has been investigated in this paper. It is shown that the wavelength of instabilities fails to be a proper criterion to be considered for this problem. However, the threshold of instability formation has been found to be appropriate for benchmarking purpose.     

Conceptual and numerical models for groundwater flow in an arid inland river basin

The Heihe River Basin (HRB) is an inland watershed in northwest China with a total area of approximately 130,000 km², stretching from the Qilian Mountains in the south to the oases and agricultural fields in the middle and further to the Gobi desert in the north bordering Mongolia. As part of a major ecohydrological research initiative to provide a stronger scientific underpinning for sustainable water management in arid ecosystems, a regional‐scale integrated ecological and hydrological model is being developed, incorporating the knowledge based on the results of environmental isotope tracer analysis and the multiscale observation datasets. The first step in the model development effort is to construct and calibrate a groundwater flow model for the middle and lower HRB where the oases and vegetation along the Heihe river corridor are highly dependent on groundwater. In this study, the software tool ‘Arc Hydro Groundwater’ is used to build and visualize a hydrogeological data model for the HRB that links all relevant spatiotemporal hydrogeological data in a unified geodatabase within the ArcGIS environment. From the conceptual model, a regional‐scale groundwater flow model has been developed using MODFLOW‐2005. Critical considerations in developing the flow model include the representation of mountainous terrains and fluvial valleys by individual model layers, treatment of aquifer heterogeneities across multiple scales and selection of proper observation data and boundary conditions for model calibration. This paper discusses these issues in the context of the Heihe River Basin, but the results and insights from this study will have important implications for other large, regional groundwater modelling studies, especially in arid and semiarid inland river basins. Copyright © 2014 John Wiley & Sons, Ltd.     

Numerical investigation of the anisotropic hydraulic conductivity behavior in heterogeneous porous media

The behavior of the mean equivalent hydraulic conductivity normal and parallel to stratification (K₁, and K₂, respectively) is studied here through Monte Carlo simulations of three-dimensional, steady-state flow in statistically anisotropic, bounded, and heterogeneous media. For water flow normal to stratification in strongly heterogeneous porous media (²_Y=3) the value of K₁ is not unique; it ranges from an arithmetic to a geometric, and finally, to a harmonic mean behavior depending on field dimensions, and medium anisotropy. For a fixed anisotropy ratio and variance of Y = ln K, the larger the distance, in the direction perpendicular to stratification, over which water flow takes place, the faster the rate at which, K_H, behavior is approached. However, even for large anisotropy ratios, harmonic mean behavior appears to be a good approximation only for aquifer thickness L₁ that is large enough to allow stratified flow to occur. For small aquifer thickness (L₁/₁<8, where ₁ is the integral scale normal to stratification) the limiting behavior, for large anisotropy ratios, appears to be, instead, that of two-dimensional flow, i.e., water flows primarily parallel to the planes of stratification. When the aquifer thickness is very small compared to the horizontal dimensions (and with relative similar integral scales in the three directions) a behavior resembling arithmetic mean conditions is exhibited, i.e., water flow takes place through heterogeneous, vertical, soil volumes. The geostatistical expressions of Desbarats (1992a) for upscaling hydraulic conductivity values were utilized and closed form empirical relations were developed for the main components of the upscaled hydraulic conductivity tensor.     

Comparison of theory and experiment for solute transport in weakly heterogeneous bimodal porous media

In this work, the influence of non-equilibrium effects on solute transport in a weakly heterogeneous medium is discussed. Three macro-scale models (upscaled via the volume averaging technique) are investigated: (i) the two-equation non-equilibrium model, (ii) the one-equation asymptotic model and (iii) the one-equation local equilibrium model. The relevance of each of these models to the experimental system conditions (duration of the pulse injection, dispersivity values…) is analyzed. The numerical results predicted by these macroscale models are compared directly with the experimental data (breakthrough curves). Our results suggest that the preasymptotic zone (for which a non-Fickian model is required) increases as the solute input pulse time decreases. Beyond this limit, the asymptotic regime is recovered. A comparison with the results issued from the stochastic theory for this regime is performed. Results predicted by both approaches (volume averaging method and stochastic analysis) are found to be consistent.     

Modelling of groundwater flow in heterogeneous porous media by finite element method

The HySuf‐FEM code (Hydrodynamic of Subsurface Flow by Finite Element Method) is proposed in this article in order to estimate the spatial variability of the transmissivity values of the Berrechid aquifer (Morocco). The calibration of the model is based on the hydraulic head, hydraulic conductivity and recharge. Three numerical tests are used to validate the model and verify its convergence. The first test case consists in using the steady analytical solution of the Poisson equation. In the second, the model has been compared with the hydrogeological system which is characterized by an unconfined monolayer (isotropic layer) and computed by using PMWIN‐MODFLOW software. The third test case is based on the comparison between the results of HySuf‐FEM and the multiple cell balance method in the aquifer system with natural boundaries case. Good agreement between the Hydrodynamic of Subsurface Flow, the numerical tests and the spatial distribution of the thickening of the hydrogeological system is deduced from the analysis and the interpretations of hydrogeological wells. Copyright © 2011 John Wiley & Sons, Ltd.