A new upscaling method for fractured porous media

We present a method to determine equivalent permeability of fractured porous media. Inspired by the previous flow-based upscaling methods, we use a multi-boundary integration approach to compute flow rates within fractures. We apply a recently developed multi-point flux approximation Finite Volume method for discrete fracture model simulation. The method is verified by upscaling an arbitrarily oriented fracture which is crossing a Cartesian grid. We demonstrate the method by applying it to a long fracture, a fracture network and the fracture network with different matrix permeabilities. The equivalent permeability tensors of a long fracture crossing Cartesian grids are symmetric, and have identical values. The application to the fracture network case with increasing matrix permeabilities shows that the matrix permeability influences more the diagonal terms of the equivalent permeability tensor than the off-diagonal terms, but the off-diagonal terms remain important to correctly assess the flow field.     

Simulation of density-driven flow in fractured porous media

We study density-driven flow in a fractured porous medium in which the fractures are represented as manifolds of reduced dimensionality. Fractures are assumed to be thin regions of space filled with a porous material whose properties differ from those of the porous medium enclosing them. The interfaces separating the fractures from the embedding medium are assumed to be ideal. We consider two approaches: (i) the fractures have the same dimension, d, as the embedding medium and are said to be d-dimensional; (ii) the fractures are considered as (d − 1)-dimensional manifolds, and the equations of density-driven flow are found by averaging the d-dimensional laws over the fracture width. We show that the second approach is a valid alternative to the first one. For this purpose, we perform numerical experiments using finite-volume discretization for both approaches. The results obtained by the two methods are in good agreement with each other.     

An inverse problem methodology to identify flow channels in fractured media using synthetic steady-state head and geometrical data

We present a methodology for identifying highly-localized flow channels embedded in a significantly less permeable medium using steady-state head and geometrical data. This situation is typical of fractured media where flows are often strongly channeled at the scales of interest (10 m–1 km). The objective is to identify both geometrical and hydraulic characteristics of the conducting structures. Channels are identified in decreasing order of importance by successive optimizations of an objective function. The identification strategy takes advantage of the hierarchical flow organization to restrict the dimension of the solution space of each individual optimization step. The characteristics of the secondary channels are strongly determined by the main flow channels. The latter are slightly modified by the secondary channels through the addition of a regularization term to the main channel characteristics in the objective function. As the objective function is strongly non-convex with numerous local minima, inversion is performed using a stochastic algorithm (simulated annealing). We assess the possibilities of the hierarchical identification strategy on simple synthetic steady-state flow configurations where hydraulic data are made up of 25 regularly spaced heads and of the boundary conditions. Those flow structures that are dominated by at most two simple channels can be identified with these head data only. Configurations comprising up to three complex and interconnected channels can still be identified with additional geometrical information including the distances of piezometers to their closest channel. The capabilities of the hierarchical identification strategy are limited to flow structures dominated by at most three equivalent flow channels. We finally discuss the perspectives of application of the method to transient-state data obtained on a more restricted number of piezometers.     

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.     

Exact transverse macro dispersion coefficients for transport in heterogeneous porous media

We study transport through heterogeneous media. We derive the exact large scale transport equation. The macro dispersion coefficients are determined by additional partial differential equations. In the case of infinite Peclet numbers, we present explicit results for the transverse macro dispersion coefficients. In two spatial dimensions, we demonstrate that the transverse macro dispersion coefficient is zero. The result is not limited on lowest order perturbation theory approximations but is an exact result. However, the situation in three spatial dimensions is very different: The transverse macro dispersion coefficients are finite – a result which is confirmed by numerical simulations we performed.     

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.     

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.     

Non-Fickian mass transport in fractured porous media

The paper provides an introduction to fundamental concepts of mathematical modeling of mass transport in fractured porous heterogeneous rocks. Keeping aside many important factors that can affect mass transport in subsurface, our main concern is the multi-scale character of the rock formation, which is constituted by porous domains dissected by the network of fractures. Taking into account the well-documented fact that porous rocks can be considered as a fractal medium and assuming that sizes of pores vary significantly (i.e. have different characteristic scales), the fractional-order differential equations that model the anomalous diffusive mass transport in such type of domains are derived and justified analytically. Analytical solutions of some particular problems of anomalous diffusion in the fractal media of various geometries are obtained. Extending this approach to more complex situation when diffusion is accompanied by advection, solute transport in a fractured porous medium is modeled by the advection-dispersion equation with fractional time derivative. In the case of confined fractured porous aquifer, accounting for anomalous non-Fickian diffusion in the surrounding rock mass, the adopted approach leads to introduction of an additional fractional time derivative in the equation for solute transport. The closed-form solutions for concentrations in the aquifer and surrounding rocks are obtained for the arbitrary time-dependent source of contamination located in the inlet of the aquifer. Based on these solutions, different regimes of contamination of the aquifers with different physical properties can be readily modeled and analyzed.     

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.     

Density-dependent dispersion in heterogeneous porous media Part I: A numerical study

In this paper, we describe carefully conducted numerical experiments, in which a dense salt solution vertically displaces fresh water in a stable manner. The two-dimensional porous media are weakly heterogeneous at a small scale. The purpose of these simulations, conducted for a range of density differences, is to obtain accurate concentration profiles that can be used to validate nonlinear models for high-concentration-gradient dispersion. In this part we focus on convergence of the computations, in numerical and statistical sense, to ensure that the uncertainty in the results is small enough.Concentration variances are computed, which give estimates of the uncertainty in local concentration values. These local variations decrease with increasing density contrast. For tracer transport, obtained longitudinal dispersivities are in accordance with analytical findings. In the case of high-density contrasts, stabilizing gravity forces counteract the growth of dispersive fingers, decreasing the effective width of the transition zone. For small log-permeability variances, the decrease of the apparent dispersivity that is found is in agreement with laboratory results for homogeneous columns.     

Application of implicit sub-time stepping to simulate flow and transport in fractured porous media

In general, the accuracy of numerical simulations is determined by spatial and temporal discretization levels. In fractured porous media, the time step size is a key factor in controlling the solution accuracy for a given spatial discretization. If the time step size is restricted by the relatively rapid responses in the fracture domain to maintain an acceptable level of accuracy in the entire simulation domain, the matrix tends to be temporally over-discretized. Implicit sub-time stepping applies smaller sub-time steps only to the sub-domain where the accuracy requirements are less tolerant and is most suitable for problems where the response is high in only a small portion of the domain, such as within and near the fractures in fractured porous media. It is demonstrated with illustrative examples that implicit sub-time stepping can significantly improve the simulation efficiency with minimal loss in accuracy when simulating flow and transport in fractured porous media. The methodology is successfully applied to density-dependent flow and transport simulations in a Canadian Shield environment, where the flow and transport is dominated by discrete, highly conductive fracture zones.     

A numerical method for two-phase flow in fractured porous media with non-matching grids

We propose a novel computational method for the efficient simulation of two-phase flow in fractured porous media. Instead of refining the grid to capture the flow along the faults or fractures, we represent the latter as immersed interfaces, using a reduced model for the flow and suitable coupling conditions. We allow for non matching grids between the porous matrix and the fractures to increase the flexibility of the method in realistic cases. We employ the extended finite element method for the Darcy problem and a finite volume method that is able to handle cut cells and matrix-fracture interactions for the saturation equation. Moreover, we address through numerical experiments the problem of the choice of a suitable numerical flux in the case of a discontinuous flux function at the interface between the fracture and the porous matrix. A wrong approximate solution of the Riemann problem can yield unphysical solutions even in simple cases.     

An efficient numerical model for incompressible two-phase flow in fractured media

Various numerical methods have been used in the literature to simulate single and multiphase flow in fractured media. A promising approach is the use of the discrete-fracture model where the fracture entities in the permeable media are described explicitly in the computational grid. In this work, we present a critical review of the main conventional methods for multiphase flow in fractured media including the finite difference (FD), finite volume (FV), and finite element (FE) methods, that are coupled with the discrete-fracture model. All the conventional methods have inherent limitations in accuracy and applications. The FD method, for example, is restricted to horizontal and vertical fractures. The accuracy of the vertex-centered FV method depends on the size of the matrix gridcells next to the fractures; for an acceptable accuracy the matrix gridcells next to the fractures should be small. The FE method cannot describe properly the saturation discontinuity at the matrix–fracture interface. In this work, we introduce a new approach that is free from the limitations of the conventional methods. Our proposed approach is applicable in 2D and 3D unstructured griddings with low mesh orientation effect; it captures the saturation discontinuity from the contrast in capillary pressure between the rock matrix and fractures. The matrix–fracture and fracture–fracture fluxes are calculated based on powerful features of the mixed finite element (MFE) method which provides, in addition to the gridcell pressures, the pressures at the gridcell interfaces and can readily model the pressure discontinuities at impermeable faults in a simple way. To reduce the numerical dispersion, we use the discontinuous Galerkin (DG) method to approximate the saturation equation. We take advantage of a hybrid time scheme to alleviate the restrictions on the size of the time step in the fracture network. Several numerical examples in 2D and 3D demonstrate the robustness of the proposed model. Results show the significance of capillary pressure and orders of magnitude increase in computational speed compared to previous works.     

An efficient numerical model for multicomponent compressible flow in fractured porous media

An efficient and accurate numerical model for multicomponent compressible single-phase flow in fractured media is presented. The discrete-fracture approach is used to model the fractures where the fracture entities are described explicitly in the computational domain. We use the concept of cross flow equilibrium in the fractures. This will allow large matrix elements in the neighborhood of the fractures and considerable speed up of the algorithm. We use an implicit finite volume (FV) scheme to solve the species mass balance equation in the fractures. This step avoids the use of Courant–Freidricks–Levy (CFL) condition and contributes to significant speed up of the code. The hybrid mixed finite element method (MFE) is used to solve for the velocity in both the matrix and the fractures coupled with the discontinuous Galerkin (DG) method to solve the species transport equations in the matrix. Four numerical examples are presented to demonstrate the robustness and efficiency of the proposed model. We show that the combination of the fracture cross-flow equilibrium and the implicit composition calculation in the fractures increase the computational speed 20–130 times in 2D. In 3D, one may expect even a higher computational efficiency.     

Comparison of Galerkin-type discretization techniques for two-phase flow in heterogeneous porous media

A comparison of Standard Galerkin, Petrov-Galerkin, and Fully-Upwind Galerkin methods for the simulation of two-phase flow in heterogeneous porous media is presented. On the basis of the coupled pressure-saturation equations, a generalized formulation for all three finite element methods is derived and analysed. For flow in homogeneous media, the Petrov-Galerkin method gives excellent results. But this method fails miserably for problems with heterogeneous media. This is because it is not able to capture correctly processes that take place at interfaces when, for instance, the capillary pressure-saturation relationship after Brooks and Corey is assumed. The Fully-Upwind Galerkin method is superior to the Petrov-Galerkin approach because it is able to give correct results for flow in homogeneous and heterogeneous media for the two models of van Genuchten and Brooks-Corey. The widely used formulation which is correct for the homogeneous case cannot be used for heterogeneous media. Instead the straightforward approach of gradp_c in combination with a chord-slope technique must be utilized.     

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.     

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.     

Bias and uncertainty in regression-calibrated models of groundwater flow in heterogeneous media

Groundwater models need to account for detailed but generally unknown spatial variability (heterogeneity) of the hydrogeologic model inputs. To address this problem we replace the large, m-dimensional stochastic vector β that reflects both small and large scales of heterogeneity in the inputs by a lumped or smoothed m-dimensional approximation γθ_∗, where γ is an interpolation matrix and θ_∗ is a stochastic vector of parameters. Vector θ_∗ has small enough dimension to allow its estimation with the available data. The consequence of the replacement is that model function f(γθ_∗) written in terms of the approximate inputs is in error with respect to the same model function written in terms of β, f(β), which is assumed to be nearly exact. The difference f(β) − f(γθ_∗), termed model error, is spatially correlated, generates prediction biases, and causes standard confidence and prediction intervals to be too small. Model error is accounted for in the weighted nonlinear regression methodology developed to estimate θ_∗ and assess model uncertainties by incorporating the second-moment matrix of the model errors into the weight matrix. Techniques developed by statisticians to analyze classical nonlinear regression methods are extended to analyze the revised method. The analysis develops analytical expressions for bias terms reflecting the interaction of model nonlinearity and model error, for correction factors needed to adjust the sizes of confidence and prediction intervals for this interaction, and for correction factors needed to adjust the sizes of confidence and prediction intervals for possible use of a diagonal weight matrix in place of the correct one. If terms expressing the degree of intrinsic nonlinearity for f(β) and f(γθ_∗) are small, then most of the biases are small and the correction factors are reduced in magnitude. Biases, correction factors, and confidence and prediction intervals were obtained for a test problem for which model error is large to test robustness of the methodology. Numerical results conform with the theoretical analysis.     

Probabilistic collocation and lagrangian sampling for advective tracer transport in randomly heterogeneous porous media

The Karhunen-Loeve (KL) decomposition and the polynomial chaos (PC) expansion are elegant and efficient tools for uncertainty propagation in porous media. Over recent years, KL/PC-based frameworks have successfully been applied in several contributions for the flow problem in the subsurface context. It was also shown, however, that the accurate solution of the transport problem with KL/PC techniques is more challenging. We propose a framework that utilizes KL/PC in combination with sparse Smolyak quadrature for the flow problem only. In a subsequent step, a Lagrangian sampling technique is used for transport. The flow field samples are calculated based on a PC expansion derived from the solutions at relatively few quadrature points. To increase the computational efficiency of the PC-based flow field sampling, a new reduction method is applied. For advection dominated transport scenarios, where a Lagrangian approach is applicable, the proposed PC/Monte Carlo method (PCMCM) is very efficient and avoids accuracy problems that arise when applying KL/PC techniques to both flow and transport. The applicability of PCMCM is demonstrated for transport simulations in multivariate Gaussian log-conductivity fields that are unconditional and conditional on conductivity measurements.