Variable-density groundwater flow and solute transport in porous media containing nonuniform discrete fractures

Variations in fluid density can greatly affect fluid flow and solute transport in the subsurface. Heterogeneities such as fractures play a major role for the migration of variable-density fluids. Earlier modeling studies of density effects in fractured media were restricted to orthogonal fracture networks, consisting of only vertical and horizontal fractures. The present study addresses the phenomenon of 3D variable-density flow and transport in fractured porous media, where fractures of an arbitrary incline can occur. A general formulation of the body force vector is derived, which accounts for variable-density flow and transport in fractures of any orientation. Simulation results are presented that show the verification of the new model formulation, for the porous matrix and for inclined fractures. Simulations of variable-density flow and solute transport are then conducted for a single fracture, embedded in a porous matrix. The simulations show that density-driven flow in the fracture causes convective flow within the porous matrix and that the high-permeability fracture acts as a barrier for convection. Other simulations were run to investigate the influence of fracture incline on plume migration. Finally, tabular data of the tracer breakthrough curve in the inclined fracture is given to facilitate the verification of other codes.     

Crossover from capillary fingering to compact invasion for two-phase drainage with stable viscosity ratios

Motivated by a wide range of applications from enhanced oil recovery to carbon dioxide sequestration, we have developed a two-dimensional, pore-level model of immiscible drainage, incorporating viscous, capillary, and gravitational effects. This model has been validated quantitatively, in the very different limits of zero viscosity ratio and zero capillary number; flow patterns from modeling agree well with experiment. For a range of stable viscosity ratios (μ_injected/μ_displaced ? 1), we have increased the capillary number, N_c, and studied the way in which the flows deviate from capillary fingering (the fractal flow of invasion percolation) and become compact for realistic capillary numbers. Results exhibiting this crossover from capillary fingering to compact invasion are presented for the average position of the injected fluid, the fluid–fluid interface, the saturation and fractional flow profiles, and the relative permeabilities. The agreement between our results and earlier theoretical predictions [Blunt M, King MJ, Scher H. Simulation and theory of two-phase flow in porous media. Phys Rev A 1992;46:7680–99; Lenormand R. Flow through porous media: limits of fractal patterns. Proc Roy Soc A 1989;423:159–68; Wilkinson D. Percolation effects in immiscible displacement. Phys Rev A 1986;34:1380–90; Xu B, Yortsos YC, Salin D. Invasion Percolation with viscous forces. Phys Rev E 1998;57:739–51] supports the validity of these general theoretical arguments, which were independent of the details of the porous media in both two and three dimensions.     

A test case for the simulation of three-dimensional variable-density flow and solute transport in discretely-fractured porous media

A test case has been developed for three-dimensional simulations of variable-density flow and solute transport in discretely-fractured porous media. The simulation domain is a low-permeability porous matrix cube containing a single non-planar fracture. The initial solute concentration is zero everywhere. A constant solute concentration is assigned to the top of the domain, which increases near-top fluid density and induces downward density-driven flow. The test case is therefore comparable to downwelling of a dense brine below a saline disposal basin or a waste repository. Numerous fingers and distinct convection cells develop early in the fracture but the fingers later coalesce and convection becomes less apparent. To help test other variable-density flow and transport models, results of the test case are presented both qualitatively (concentration contours and velocity fields) and quantitatively (penetration depth, mass flux, total mass stored, maximum fracture and matrix velocity).     

Exhalation of radon and thoron: the question of the effect of thermal gradients in soil

Experimental measurements show no evidence for diurnal variation in exhalation of radon and thoron from soil due to convection induced by thermal gradients in the top few decimeters of soil as suggested by several authors. Estimates based on convective calculations indicate that even in the unlikely event of either vertical temperature gradients large enough to cause vertical instability, or horizontal gradients sufficient to cause significant convection, any effect would be too small to be detected. These same calculations suggest that it is difficult to conceive of cases involving typical thermal gradients in unfractured porous media such as soil where thermally induced convection would play an important role in transport of radon or thoron.     

The influence of precipitate formation on the chemical oxidation of TCE DNAPL with potassium permanganate

A three-dimensional two-phase flow model is coupled to a non-linear reactive transport model to study the efficacy of potassium permanganate treatment on dense, non-aqueous phase liquid (DNAPL) source removal in porous media. A linear relationship between the soil permeability (k) and concentration of manganese dioxide precipitate ([MnO_2(s)]), k = k_o + S_rind [MnO_2(s)], is utilized to simulate nodal permeability reductions due to precipitate formation. Using published experimental column studies, an S_rind = −5.5 × 10⁻¹⁶ m² L/mg was determined for trichloroethylene (TCE) DNAPL. This S_rind was then applied to treatment simulations on three-dimensional TCE DNAPL source zones comprising either DNAPL at residual saturations, or DNAPL at pooled saturations.     

Modeling changes in the permeability in a gypsified fractured-porous rock massif

A method is proposed for evaluating rock hydraulic conductivity (k _x;y , LT^?1) in grid model blocks to be used in finite-difference simulation of mass exchange in gypcified fractured-porous rocks containing groundwater. The values of k _x and k _y were determined taking into account the dominating form of transport (advection or transverse dispersion), contributing most to gypsum dissolution. The method was used to simulate groundwater flow in a vertical section of gypsified rocks in the base of the dam at the Lower-Kafirniganskii Hydropower System. The changes in rock permeability in the base of the dam were used to evaluate whether there is a need to take seepage-control measures (cement-grout curtain, blanket, and a gypsum fill over the blanket). A leaching zone was found to form near the curtain and blanket during the simulation period (100 years) under background seepage conditions. In this period, the total water flow in the lower pool can increase by no more than 10–15%. Under such conditions, the “seepage-control” measures can enhance the reliability of the dam at the Lower-Kafirniganskii Hydropower System.     

Hot springs mediate spatial exchange of heat and mass in the enclosed sediment domain: A stability perspective

In many natural environments, such as in underwater hot springs and hydrothermal vents, thermal gradients are accompanied with changes in the concentration of chemical compounds transported to the seawater, causing the so-called double-diffusive, mixed convection. To study the physical scenarios in such systems, a vertical channel filled with a porous medium saturated with saline water is considered. The motion in the sediment-filled channel is induced by two buoyancy forces and an external pressure gradient, similar to the situation in a vent with an upward flow direction. The fluid flow has been modeled by an extended Darcy model, and the flow instability mechanisms have been studied numerically. The linear stability analysis is performed considering a wide range of Darcy number (Da = 10⁻⁵ -10⁻⁸). The instability boundary curve showed three distinct dynamic regimes: (i) Rayleigh-Taylor (R-T), (ii) log-log non-linear variation, and (iii) log-log linear variation. The domain of different regimes were sensitive to external pressure gradient as well as permeability. Similar to cross-diffusive natural convection in pure viscous fluids, a linear relationship between logarithmic absolute values of critical thermal Rayleigh number (∣Ra_T∣) and solute Rayleigh number (Ra_C) is found in the third regime. Based on the permeability, for any solute Rayleigh number (Ra_C), there existed a minimum value of Reynolds number (Re), below which R-T type of instability appeared. Above this minimum value, the instability was due to two buoyancy forces, known as buoyant instability. Simulations of secondary flow via energy analysis demonstrated the development of complex dynamics at the critical state in all three regimes characterized by transition of multi to uni-cellular structures and vice verse.     

Analytical modeling of one-dimensional diffusion in layered systems with position-dependent diffusion coefficients

Diffusion in stratified porous media is common in the natural environment. The objective of this study is to develop analytical solutions for describing the diffusion in layered porous media with a position-dependent diffusion coefficient within each layer. The orthogonal expansion technique was used to solve a one-dimensional multi-layer diffusion equation in which the diffusion coefficient is expressed as a segmented linear function of positions in the porous media. The behavior of the solutions is illustrated using several examples of a three-layer system, with constant diffusion coefficient α₁ in layer 1 (0 < x < l₁), α₃ in layer 3(l₂ < x < l₃), and a linearly position-dependent diffusion coefficient α₁(1 + Δ(x − l₁)/(l₂ − l₁)) in the center layer (Δ = α₃/α₁ − 1). Because of the asymmetry of the layered system, the diffusion and related concentration distributions are also asymmetrical. For a given Δ value, the smaller the value of (l₂ − l₁)/l₃, the more significant the accumulation of concentration in the middle transition zone (l₁ < x < l₂), the sharper the change in the concentration profile of spatial distribution. Therefore, transition between two layers has significant effects on diffusion.     

Stability of buoyancy opposed mixed convection in a vertical channel and its dependence on permeability

Non-Darcy mixed convective flow of water due to external pressure gradient and buoyancy opposed forces are considered in a vertical channel filled with porous medium, which can be either isotropic or anisotropic. The linear theory of stability analysis has been used to numerically investigate the dependence of the transition behavior of the fully developed basic flow on the permeability of the medium. Numerical experiments indicate that mainly two main instability modes appear: Rayleigh–Taylor (R–T) and buoyant instability. For Darcy numbers (Da) ?10⁻⁹, R–T instability dominates within the entire Reynolds number (Re) range considered here. It was also found that for the same Re, the fully developed base flow is highly unstable (stable) for porous media with high (low) permeability. Further, it was seen that the disturbance isotherm cells migrate from the channel walls toward the centerline when permeability is reduced. Reducing the permeability by one order of magnitude (corresponding to a decrease of Darcy number from 10⁻⁶ to 10⁻⁷) increases base flow stability approximately 20-fold. For higher Reynolds numbers, buoyant, mixed and shear instability of the basic flow were found when Da was increased from 10⁻⁷ to 10⁻³. However, for cases in which permeability and porosity behaved as suggested by Carman–Kozeny relation (CKR), buoyant stability was the only mode of instability. Critical values of the Rayleigh (Ra) and Darcy (Da) numbers in the R–T mode of instability were related to each other by the hyperbolic function RaDa = −2.465.     

An examination of τc-Pd earthquake early warning method using a strong-motion building array

The use of characteristic period τ_c and peak displacement amplitude Pd of the initial P wave in earthquake early warning (EEW) was proposed by Wu and Kanamori 1, 2, 3 and 4. Here we apply this approach to strong-motion records from a building sensor array installed in Taitung County, Taiwan. This building was damaged during the 2006 M_w=6.1 Taitung earthquake with a peak ground velocity (PGV) of up to 38.4 cm/s at an epicentral distance of 14.5 km. According to our analysis, the peak displacement amplitude Pd is a better indicator for the destructiveness of an earthquake than τ_c because τ_c is more sensitive to the signal-to-noise ratio (SNR) than Pd. In accordance with previous studies, only the structurally damaging Taitung earthquake generated a Pd value larger than 0.5 cm (a threshold for identifying damaging events). Using Pd as an indicator for destructive earthquakes does not lead to missing or false alarms for EEW purposes.     

Testing density-dependent groundwater models: two-dimensional steady state unstable convection in infinite,finite and inclined porous layers

This study proposes the use of several problems of unstable steady state convection with variable fluid density in a porous layer of infinite horizontal extent as two-dimensional (2-D) test cases for density-dependent groundwater flow and solute transport simulators. Unlike existing density-dependent model benchmarks, these problems have well-defined stability criteria that are determined analytically. These analytical stability indicators can be compared with numerical model results to test the ability of a code to accurately simulate buoyancy driven flow and diffusion. The basic analytical solution is for a horizontally infinite fluid-filled porous layer in which fluid density decreases with depth. The proposed test problems include unstable convection in an infinite horizontal box, in a finite horizontal box, and in an infinite inclined box. A dimensionless Rayleigh number incorporating properties of the fluid and the porous media determines the stability of the layer in each case. Testing the ability of numerical codes to match both the critical Rayleigh number at which convection occurs and the wavelength of convection cells is an addition to the benchmark problems currently in use. The proposed test problems are modelled in 2-D using the SUTRA [SUTRA––A model for saturated–unsaturated variable-density ground-water flow with solute or energy transport. US Geological Survey Water-Resources Investigations Report, 02-4231, 2002. 250 p] density-dependent groundwater flow and solute transport code. For the case of an infinite horizontal box, SUTRA results show a distinct change from stable to unstable behaviour around the theoretical critical Rayleigh number of 4π² and the simulated wavelength of unstable convection agrees with that predicted by the analytical solution. The effects of finite layer aspect ratio and inclination on stability indicators are also tested and numerical results are in excellent agreement with theoretical stability criteria and with numerical results previously reported in traditional fluid mechanics literature.     

Monotonic solution of flow and transport problems in heterogeneous media using Delaunay unstructured triangular meshes

Transport problems occurring in porous media and including convection, diffusion and chemical reactions, can be well represented by systems of Partial Differential Equations. In this paper, a numerical procedure is proposed for the fast and robust solution of flow and transport problems in 2D heterogeneous saturated media. The governing equations are spatially discretized with unstructured triangular meshes that must satisfy the Delaunay condition. The solution of the flow problem is split from the solution of the transport problem and it is obtained with an approach similar to the Mixed Hybrid Finite Elements method, that always guarantees the M-property of the resulting linear system. The transport problem is solved applying a prediction/correction procedure. The prediction step analytically solves the convective/reactive components in the context of a MAST Finite Volume scheme. The correction step computes the anisotropic diffusive components in the context of a recently proposed Finite Elements scheme. Massa balance is locally and globally satisfied in all the solution steps. Convergence order and computational costs are investigated and model results are compared with literature ones.     

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.     

Exploring the nature of non-Fickian transport in laboratory experiments

Observations of non-Fickian transport in sandbox experiments [Levy M, Berkowitz B. Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. J Contam Hydrol 2003;64:203–26] were analyzed previously using a power law tail ψ(t) ∼ t^−1−β with 0 < β < 2 for the spectrum of transition times comprising a tracer plume migration. For each sandbox medium a choice of β resulted in an excellent fit to the breakthrough curve (BTC) data, and the value of β decreased slowly with increasing flow velocity. Here, the data are reanalyzed with the full spectrum of ψ(t) gleaned from analytical calculations [Cortis A, Chen Y, Scher H, Berkowitz B. Quantitative characterization of pore-scale disorder effects on transport in “homogeneous” granular media. Phys Rev E 2004;10(70):041108. doi: 10.1103/PhysRevE.70.041108], numerical simulations [Bijeljic B, Blunt MJ. Pore-scale modeling and continuous time random walk analysis of dispersion in porous media. Water Resour Res 2006;42:W01202. doi: 10.1029/2005WR004578] and permeability fields [Di Donato G, Obi E-O, Blunt MJ. Anomalous transport in heterogeneous media demonstrated by streamline-based simulation. Geophys Res Lett 2003;30:1608–12s. doi: 10.1029/2003GL017196]. We represent the main features of the full spectrum of transition times with a truncated power law (TPL), ψ(t) ∼ (t₁ + t)^−1−βexp(−t/t₂), where t₁ and t₂ are the limits of the power law spectrum. An excellent fit to the entire BTC data set, including the changes in flow velocity, for each sandbox medium is obtained with a single set of values of t₁, β, t₂. The influence of the cutoff time t₂ is apparent even in the regime t < t₂. Significantly, we demonstrate that the previous apparent velocity dependence of β is a result of choosing a pure power law tail for ψ(t). The key is the change in the log–log slope of the TPL form of ψ(t) with a shifting observational time window caused by the change in the mean velocity. Hence, the use of the full spectrum of ψ(t) is not only necessary for the transition to Fickian behavior, but also to account for the dynamics of these laboratory observations of non-Fickian transport.     

Dimensional analysis of two-phase flow including a rate-dependent capillary pressure–saturation relationship

The macroscopic modelling of two-phase flow processes in subsurface hydrosystems or industrial applications on the Darcy scale usually requires a constitutive relationship between capillary pressure and saturation, the P_c(S_w) relationship. Traditionally, it is assumed that a unique relation between P_c and S_w exists independently of the flow conditions as long as hysteretic effects can be neglected. Recently, this assumption has been questioned and alternative formulations have been suggested. For example, the extended P_c(S_w) relationship by Hassanizadeh and Gray [Hassanizadeh SM, Gray WG. Mechanics and thermodynamics of multiphase flow in porous media including interphase boundaries. Adv Water Resources 1990;13(4):169–86] proposes that the difference between the phase pressures to the equilibrium capillary pressure is a linear function of the rate of change of saturation, thereby introducing a constant of proportionality, the coefficient τ. It is desirable to identify cases where the extended relationship needs to be considered. Consequently, a dimensional analysis is performed on the basis of the two-phase balance equations. In addition to the well-known capillary and gravitational number, the dimensional analysis yields a new dimensionless number. The dynamic number Dy quantifies the ratio of dynamic capillary to viscous forces. Relating the dynamic to the capillary as well as the gravitational number gives the new numbers DyC and DyG, respectively. For given sets of fluid and porous medium parameters, the dimensionless numbers Dy and DyC are interpreted as functions of the characteristic length and flow velocity. The simulation of an imbibition process provides insight into the interpretation of the characteristic length scale. The most promising choice for this length scale seems to be the front width. We conclude that consideration of the extended P_c(S_w) relationship may be important for porous media with high permeability, small entry pressure and high coefficient τ when systems with a small characteristic length (e.g. steep front) and small characteristic time scale are under investigation.     

An alternative approach for the Istanbul earthquake early warning system

Two recent catastrophic earthquakes that struck the Marmara Region on 17 August 1999 (M_w=7.4) and 12 November 1999 (M_w=7.2) caused major concern about future earthquake occurrences in Istanbul and the Marmara Region. As a result of the preparations for an expected earthquake may occur around Istanbul region, an earthquake early warning system has been established in 2002 with a simple and robust algorithm, based on the exceedance of specified thresholds of time domain amplitudes and the cumulative absolute velocity (CAV) levels (Erdik et al., 2003 [1]). In order to improve the capability of Istanbul earthquake early warning system (IEEWS) for giving early warning of a damaging earthquake in the Marmara Region, we explored an alternative approach with the use of a period parameter (τ_c) and a high-pass filtered vertical displacement amplitude parameter (Pd) from the initial 3 s of the P waveforms as proposed by Kanamori (2005) [2] and Wu and Kanamori (2005) 3 and 4. The empirical relationships both between τ_c and moment magnitude (M_w), and between Pd and peak ground velocity (PGV) for the Marmara Region are presented. These relationships can be used to detect a damaging earthquake within seconds after the arrival of P waves, and can provide on-site warning in the Marmara Region.     

Finite element techniques for convective-diffusive transport in porous media

Operator-splitting techniques are applied to convective-diffusive transport problems in porous media. The convection is treated by applying a modified method of characteristics to time-step along the characteristics of the convective part of the flow. The nonsymmetry in the spatial operator is addressed via a Petrov-Galerkin method which uses a test function to achieve stability through a balancing of the remaining convection, the diffusion, and any possible reaction terms. The use of time-stepping along characteristics allows the use of large time-steps in a stable but accurate fashion. If local phenomena are important, self-adaptive local grid refinement techniques can be coupled with the operator splitting.     

Forward induced seismic hazard assessment: application to a synthetic seismicity catalogue from hydraulic stimulation modelling

The M _w 3.2-induced seismic event in 2006 due to fluid injection at the Basel geothermal site in Switzerland was the starting point for an ongoing discussion in Europe on the potential risk of hydraulic stimulation in general. In particular, further development of mitigation strategies of induced seismic events of economic concern became a hot topic in geosciences and geoengineering. Here, we present a workflow to assess the hazard of induced seismicity in terms of occurrence rate of induced seismic events. The workflow is called Forward Induced Seismic Hazard Assessment (FISHA) as it combines the results of forward hydromechanical-numerical models with methods of time-dependent probabilistic seismic hazard assessment. To exemplify FISHA, we use simulations of four different fluid injection types with various injection parameters, i.e. injection rate, duration and style of injection. The hydromechanical-numerical model applied in this study represents a geothermal reservoir with preexisting fractures where a routine of viscous fluid flow in porous media is implemented from which flow and pressure driven failures of rock matrix and preexisting fractures are simulated, and corresponding seismic moment magnitudes are computed. The resulting synthetic catalogues of induced seismicity, including event location, occurrence time and magnitude, are used to calibrate the magnitude completeness M _c and the parameters a and b of the frequency-magnitude relation. These are used to estimate the time-dependent occurrence rate of induced seismic events for each fluid injection scenario. In contrast to other mitigation strategies that rely on real-time data or already obtained catalogues, we can perform various synthetic experiments with the same initial conditions. Thus, the advantage of FISHA is that it can quantify hazard from numerical experiments and recommend a priori a stimulation type that lowers the occurrence rate of induced seismic events. The FISHA workflow is rather general and not limited to the hydromechanical-numerical model used in this study and can therefore be applied to other fluid injection models.     

A computational model for coupled multiphysics processes of CO2 sequestration in fractured porous media

In this paper, a computational model for the simulation of coupled hydromechanical and electrokinetic flow in fractured porous media is introduced. Particular emphasis is placed on modeling CO₂ flow in a deformed, fractured geological formation and the associated electrokinetic flow. The governing field equations are derived based on the averaging theory and the double porosity model. They are solved numerically with a mixed discretization scheme, formulated on the basis of the standard Galerkin finite element method, the extended finite element method, the level-set method and the Petrov–Galerkin method. The standard Galerkin method is utilized to discretize the equilibrium and the diffusive dominant field equations, and the extended finite element method, together with the level-set method and the Petrov–Galerkin method, are utilized to discretize the advective dominant field equations. The level-set method is employed to trace the CO₂ plume front, and the extended finite element method is employed to model the high gradient in the saturation field front. The proposed mixed discretization scheme leads to a convergent system, giving a stable and effectively mesh-independent model. The accuracy and computational efficiency of the proposed model is evaluated by verification and numerical examples. Effects of the fracture spacing on the CO₂ flow and the streaming potential are discussed.     

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.