Analysis of fracture network connectivity using percolation theory

Connectivity aspects of fracture networks are analyzed in terms of percolation theory. These aspects are of fundamental importance in characterization, exploitation, and management of fractured formations. General connectivity and power law relationships are determined that characterize the density of fractures and average number of intersections per fracture necessary to ensure network connectivity, the likelihood of a fractured formation being hydraulically connected, and the probability that any specific fracture is connected to the conducting portion of the network. Monte Carlo experiments with a two-dimensional fracture network model confirm the percolation theory predictions. These relationships may prove useful in formulating theoretically tractable approximations of fracture nerworks that capture the essential system properties.     

A priori error estimates for the numerical solution of a coupled geomechanics and reservoir flow model with stress-dependent permeability

In this paper we consider the numerical solution of a coupled geomechanics and a stress-sensitive porous media reservoir flow model. We combine mixed finite elements for Darcy flow and Galerkin finite elements for elasticity. This work focuses on deriving convergence results for the numerical solution of this nonlinear partial differential system. We establish convergence with respect to the L ²-norm for the pressure and for the average fluid velocity and with respect to the H ¹-norm for the deformation. Estimates with respect to the L ²-norm for mean stress, which is of special importance since it is used in the computation of permeability for poro-elasticity, can be derived using the estimates in the H ¹-norm for the deformation. We start by deriving error estimates in a continuous-in-time setting. A cut-off operator is introduced in the numerical scheme in order to derive convergence. The spatial grids for the discrete approximations of the pressure and deformation do not need be the same. Theoretical convergence error estimates in a discrete-in-time setting are also derived in the scope of this investigation. A numerical example supports the convergence results.     

含水层层状非均质对地下水流系统的影响

区域尺度上含水层非均质具有复杂的结构性和随机性,难以准确刻画,造成非均质对区域地下水流系统的影响机制研究不够深入。本文以鄂尔多斯盆地白垩系地下水流系统为研究实例,选择典型剖面,采用剖面二维随机数值模拟方法,通过对比不同非均质刻画方法下地下水流场的变化,探讨含水层层状非均质对地下水流系统的影响机制。结果显示,均质条件下模型各向异性（含水层水平和垂向渗透系数比值Kh/Kv）取值为1000时,地下水流场与实际条件较为接近;非均质条件下,渗透系数方差取值0.91,水平相关长度取值5000 m,Kh/Kv取值150时,接近实际条件。研究表明,在大尺度地下水流模拟研究中,采用水平相关长度、渗透系数方差和各向异性值三个变量生成的随机场能很好地刻画含水层的层状非均质特征及其对水流系统的影响控制作用。由于含水层不同尺度层状非均质的叠加效应,采用均质各向异性介质等效概化含水层层状非均质性会造成等效各向异性值偏大失真的效应。     

非饱和透明土优先流迁移规律分析

为实现土的非饱和优先流迁移可视化,设计立柱装置进行非饱和渗流试验,采用透明土与数字图像处理技术,建立归一化像素强度与透明土饱和度的关系,在此基础上,通过室内模型试验,研究优先流路径连通性与相邻优先流路径旋转角对非饱和土优先流迁移的影响。试验结果表明：基于图像灰度像素强度表征非饱和透明土饱和度方法是可行的;全连通优先流（O-O型）与上连通优先流（O-C型）剖面呈现T型,中心轴剖面饱和度与边缘差异明显,下连通优先流路径（C-O型）中土体油压与基质势不足以使流体进入优先流路径形成优先流,入渗趋势与均匀流一致,O-C型优先流的稳定入渗率和湿润锋推移速度分别为C-O型的1.5倍和1.4倍;相邻O-C型优先流之间区域形成新的优先流,土体达到较高饱和度,增长速率随旋转角的增加而减小,优先流转角为90?、60?、30?时稳定入渗速率分别为均匀流的1.5、1.3、1.2倍;流体受重力影响,转角小的优先流仅沿路径一侧水平渗透,但湿润锋推移速度分别为均匀流的1.3、1.4、1.5倍,相邻优先流相互作用减弱,难以形成新的优先流。     

土壤非饱和带中的优先流

优先流是近年来研究土壤水运动所提出的术语。在某些条件下,水分在非饱和带中的运移不是简单的一维垂直向下流动,而是越过大部分土壤的体积沿着优先途径流动。优先流通常分为三类:大孔隙流,漏斗流,指流。本文分别介绍了这三类流动形式,并回顾了有关优先流的理论和实验的研究成果。     

Dispersive Particle Transport: Identification of Macroscale Behavior in Heterogeneous Stratified Subsurface Flows

Dispersive mass transport processes in naturally heterogeneous geological formations (porous media) are investigated based on a particle approach to mass transport and on its numerical implementation using LPT3D, a Lagrangian Particle Tracking 3D code. We are currently using this approach for studying microscale and macroscale space–time behavior (advection, diffusion, dispersion) of tracer plumes, solutes, or miscible fluids, in 1,2,3-dimensional heterogeneous and anisotropic subsurface formations (aquifers, petroleum reservoirs). Our analyses are based on a general advection-diffusion model and numerical scheme where concentrations and fluxes are discretized in terms of particles. The advection-diffusion theory is presented in a probabilistic framework, and in particular, a numerical analysis is developed for the case of advective transport and rotational flows (numerical stability of the explicit Euler scheme). The remainder of the paper is devoted to the behavior of concentration, mass flux density, and statistical moments of the transported tracer plume in the case of heterogeneous steady flow fields, where macroscale dispersion occurs due to geologic heterogeneity and stratification. We focus on the case of perfectly stratified or multilayered media, obtained by generating many horizontal layers with a purely random transverse distribution of permeability and horizontal velocity. In this case, we calculate explicitly the exact mass concentration field C(x, t), mass flux density field f(x, t), and moments. This includes spatial moments and dispersion variance ² _x (t) on a finite domain L, and temporal moments on a finite time scale T, e.g., the mass variance of arrival times ² _T (x). The moments are related to flux concentrations in a way that takes explicitly into account finite space–time scales of analysis (time-dependent tracer mass; spatially variable flow through mass). The multilayered model problem is then used in numerical experiments for testing different ways of recovering information on tracer plume migration, dispersion, concentration and flux fields. Our analyses rely on a probabilistic interpretation that emerges naturally from the particle approach; it is based on spatial moments (particle positions), temporal moments (mass weighted arrival times), and probability densities (both concentrations and fluxes). Finally, as an alternative to direct estimations of the flux and concentration fields, we formulate and study the Moment Inverse Problem. Solving the MIP yields an indirect method for estimating the space–time distribution of flux concentrations based on observed or estimated moments of the plume. The moments may be estimated from field measurements, or numerically computed by particle tracking as we do here.     

Direct numerical simulation of a turbulent stably stratified air flow above a wavy water surface

The influence of the roughness of the underlaying water surface on turbulence is studied in a stably stratified boundary layer (SSBL). Direct numerical simulation (DNS) is conducted at various Reynolds (Re) and Richardson (Ri) numbers and the wave steepness ka. It is shown that, at constant Re, the stationary turbulent regime is set in at Ri below the threshold value Ri_c depending on Re. At Ri > Ri_c, in the absence of turbulent fluctuations near the wave water surface, three-dimensional quasiperiodical structures are identified and their threshold of origin depends on the steepness of the surface wave on the water surface. This regime is called a wave pumping regime. The formation of three-dimensional structures is explained by the development of parametric instability of the disturbances induced by the surface water in the air flow. The DNS results are quite consistent with prediction of the theoretical model of the SSBL flow, in which solutions for the disturbances of the fields of velocity and temperature in the wave pumping regime are found to be a solution of a two-dimensional linearized system with the heterogeneous boundary condition, which is caused by the presence of the surface wave. In addition to the turbulent fluctuations, the three-dimensional structures in the wave pumping regime provide for the transfer of impulse and heat, i.e., the increase in the roughness of the water–air boundary caused by the presence of waves intensifies the exchange in the SSBL.     

含水层渗透系数空间变异性对地下水数值模拟的影响

地下水数值模拟是目前定量研究地下水水量和水质的重要手段。使用基于随机理论的MonteCarlo方法来进行地下水数值模拟。这种方法能较好地考虑水文地质参数的空间变异性。主要将MonteCarlo方法和确定性模型模拟方法的模拟结果在渗透系数场、水头场、速度场和浓度场等方面进行了比较。结果表明:在模拟三维非均质含水层中的溶质运移问题时,充分考虑了含水层渗透系数空间变异性的MonteCarlo法比确定性方法更为有效,模拟精度提高了很多,且对模拟误差及误差来源有合理的数学解释。     

Impact of Data Density and Geostatistical Simulation Technique on the Estimation of Residence Times in a Synthetic Two-dimensional Aquifer

Connectivity patterns of heterogeneous porous media are important in the estimation of groundwater residence time distributions (RTDs). Understanding the connectivity patterns of a hydraulic conductivity ( \(K\) ) field often requires knowledge of the entire aquifer, which is not practical. As such, the method used to estimate unknown \(K\) values using known \(K\) values is important. This study investigates how varying levels of conditioning data and four simulation techniques, one multi-Gaussian and three multi-point, are able to recreate key \(K\) field features and connectivity patterns of a synthetic two-dimensional bimodal distributed ln( \(K\) ) field with highly connected high \(K\) features. These techniques are then assessed in the context of RTD estimation. It was found that the multi-Gaussian technique presented a bias towards earlier travel times with increased conditioning data. This was due to the inability of the method to recreate multiple scales of connecting features. Of the multi-point methods investigated, the facies method was unable to predict early arrival times. The use of a continuous variable training image produced good fits to the observed residence time distribution with a high number of conditioning points. The ability of the methods to predict the shape of residence time distributions appears to be related to their ability to reproduce the connection patterns of higher \(K\) features.     

JHomogenizer: a computational tool for upscaling permeability for flow in heterogeneous porous media

This paper presents an object-oriented programming approach for the design of numerical homogenization programs, called JHomogenizer. It currently includes five functional modules to compute effective permeability and simple codes for computing solutions for flow in porous media. Examples with graphical output are shown to illustrate some functionalities of the program. A series of numerical examples demonstrates the effectiveness of the methodology for two-phase flow in heterogeneous reservoirs. The software is freely available, and the open architecture of the program facilitates further development and can adapt to suit specific needs easily and quickly.     

An efficient approach to generating DEM samples having coincident void ratios with actual sand

The ranges of initial void ratios that can be achieved for numerical samples with the same grading as actual sand under gravity are explored using three-dimensional discrete element method (DEM). A corrected compressible accumulation model is proposed to predict the packing density of sand. Compared with the measured results, the errors of the prediction results are small. For three-dimensional DEM samples, when the ratio of sample size to particle size L/R is equal to or larger than 30, the maximum and the minimum void ratios are minimally affected by the sample size L. The maximum void ratios and the minimum void ratios of numerical samples with spherical particles are much smaller than those of actual sands. The generation method for numerical samples with non-spherical particles is proposed based on the theory of CPM. The variations of void ratios of numerical samples with non-spherical particles are wider than actual sands. So it can meet all the need of DEM simulation on the void ratio of actual sand.     

On the velocity covariance for steady flows in heterogeneous porous formations and its application to contaminants transport

We consider groundwater steady flow in a heterogeneous porous formation of random and stationary log-conductivity Y = ln K, characterized by the mean 〈Y〉, and the two point correlation function C _Y which in turn has finite, and different horizontal and vertical integral scales I and I _v , respectively. The fluid velocity V, driven by a given head drop applied at the boundary, has constant mean value U ≡ (U, 0, 0). Approximate explicit analytical expressions for transverse velocity covariances are derived. The adopted methodology follows the approach developed by Dagan and Cvetkovic (Spatial moments of kinetically sorbing plume in a heterogeneous aquifers, Water Resour. Res. 29 (1993) 4053) to obtain a similar result for the longitudinal velocity covariance. Indeed, the approximate covariances of transverse velocities are determined by requiring that they have the exact first order variances as well as zero integral scale (G. Dagan, Flow and Transport in Porous Formations (Springer, 1989)) , and provide the exact asymptotic limits of the displacement covariance of the fluid particles obtained by Russo (On the velocity covariance and transport modeling in heterogeneous anisotropic porous formations 1. Saturated flow, Water Resour. Res., 31 (1995) 129). Comparisons with numerical results show that the proposed expressions compare quite well in the early time regime, and for Ut/I >100. Since most of the applications, like assessing the effective mobility of contaminants or quantifying the potential hazards of nuclear repositories, require predictions over higher times the proposed approximate expressions provide acceptable results. The main advantage related to such expressions is that they allow obtaining closed analytical forms of spatial moments pertaining to kinetically sorbing contaminant plumes avoiding the very heavy computational effort which is generally demanded. For illustration purposes, we consider the movement of one contaminant species, and show how our approximate spatial moments compare with the numerical simulations.     

Analyzing unsaturated flow patterns in fractured rock using an integrated modeling approach

Characterizing percolation patterns in unsaturated fractured rock has posed a greater challenge to modeling investigations than comparable saturated zone studies due to the heterogeneous nature of unsaturated media and the great number of variables impacting unsaturated flow. An integrated modeling methodology has been developed for quantitatively characterizing percolation patterns in the unsaturated zone of Yucca Mountain, Nevada (USA), a proposed underground repository site for storing high-level radioactive waste. The approach integrates moisture, pneumatic, thermal, and isotopic geochemical field data into a comprehensive three-dimensional numerical model for analyses. It takes into account the coupled processes of fluid and heat flow and chemical isotopic transport in Yucca Mountain’s highly heterogeneous, unsaturated fractured tuffs. Modeling results are examined against different types of field-measured data and then used to evaluate different hydrogeological conceptualizations through analyzing flow patterns in the unsaturated zone. In particular, this model provides clearer understanding of percolation patterns and flow behavior through the unsaturated zone, both crucial issues in assessing repository performance. The integrated approach for quantifying Yucca Mountain’s flow system is demonstrated to provide a practical modeling tool for characterizing flow and transport processes in complex subsurface systems.     

基于转移概率随机模拟的DNAPL运移二维和三维数值模拟对比研究

非均质介质的空间维度变化对重非水相流体（DNAPL） 的运移具有重要的影响。在充分考虑地质体的空间连续 性、不对称性以及各向异性等特征的基础上,采用基于马尔可夫链的转移概率（transition probability） 模型来构建非均质 随机场。该文通过TMVOC-MP软件来模拟DNAPL在非均质介质中的运移规律,探讨非均质随机场的水平空间连续性、 空间维度变化以及侧向运移过程对DNAPL运移的影响。结果表明,介质的水平空间连续性越好,DNAPL在水平方向的 迁移范围越大,在垂向的迁移范围越小;相比于三维模型,二维模型中DNAPL在水平方向的展布更大、在透镜体上的蓄 积量更多,在实际应用中以二维模型代替三维模型会加大模拟结果与实际污染情况之间的误差;侧向运移过程削弱了单 个平面的非均质性对DNAPL运移的控制,当存在侧向运移时,DNAPL绕过透镜体所运移的距离以及在透镜体上的蓄积 量会相应减小。     

Determining the active region model parameter from dye staining experiments for characterizing the preferential flow heterogeneity in unsaturated soils

The active region model (ARM) has been developed as a practical and effective approach for characterizing and representing preferential flow patterns in unsaturated soils. However, studies on methods to determine the ARM parameter (γ) are very limited. The major objective of this work was to refine the methods for determining the ARM parameter (γ) using the data from field-scale dye staining experiments. For this purpose, 13 field-scale dye staining experiments were conducted in silty clay, loam and sand with various initial and boundary conditions. The distributions of soil water content and fraction of stained region were measured to determine the ARM parameter (γ). Three determination methods of ARM parameter (γ) were presented according to the different distribution patterns of the soil water content of stained region. The efficiency of these approaches was demonstrated with the results from 13 field-scale dye staining tests that were directly related to preferential flow patterns. Impacts of soil texture on soil water redistribution, and effects of initial soil water content and infiltration amount on the preferential flow heterogeneity were also discussed in this research.     

Assessment of Discharge through a Dike Breach and Simulation of Flood Wave Propagation

This paper presents a simple and fast method to calculate flow through a dike breach. The approach was based on two-dimensional numerical simulations of idealized dike breakages at straight river-sections. As a result, computation of discharge through a breach can be achieved by use of the new developed formula (denoted as dike break formula). Furthermore, a methodology that combines one-dimensional hydrodynamic modelling, the dike break formula and a simple GIS-based method to estimate inundation areas is described. This fast and easy-to-handle tool can be utilized for near real-time forecasting or evacuation decisions. Detailed predictions were made for a number of flood and dike break scenarios at the River Rhine to prove the accuracy of the new method compared with two-dimensional numerical models.     

Numerical evaluation of multicomponent cation exchange reactive transport in physically and geochemically heterogeneous porous media

Experimental evidence and stochastic studies strongly show that the transport of reactive solutes in porous media is significantly influenced by heterogeneities in hydraulic conductivity, porosity, and sorption parameters. In this paper, we present Monte Carlo numerical simulations of multicomponent reactive transport involving competitive cation exchange reactions in a two-dimensional vertical physically and geochemically heterogeneous medium. Log hydraulic conductivity, log K, and log cation exchange capacity (log CEC) are assumed to be random Gaussian functions with spherical semivariograms. Random realizations of log K and log CEC are used as input data for the numerical simulation of multicomponent reactive transport with CORE^2D, a general purpose reactive transport code. Longitudinal features of the fronts of reactive and conservative species are computed from the temporal and spatial moments of depth-averaged concentrations. Monte Carlo simulations show that: (1) the displacement of reactive fronts increases with increasing variance of log K, while it decreases with the variance of log CEC; (2) second-order spatial moments increase with increasing variances of log K and log CEC; (3) uncertainties in the mean arrival time are largest (smallest) for negatively (positively) correlated log K and Log CEC; (4) cations undergoing competitive cation exchange exhibit different apparent velocities and retardation factors due to both physical and geochemical heterogeneities; and (5) the correlation between log K and log CEC affects significantly apparent cation retardation factors in heterogeneous aquifers.     

Estimation of stream water quality parameter using regime channel theory

For predicting and forecasting fate of non-conservative pollutants downstream from source using advection–dispersion–decay equation (ADDE), estimation of three parameters; mean flow velocity (U), longitudinal dispersion co-efficient (D _L) and decay rate co-efficient (λ), is required a priori. In this three parameters model, estimation of D _L holds difficulties and draws interest towards it. The empirical formulae use the field and experimental data of channel and flow characteristics to estimate D _L. In this paper, an innovative approach has been proposed towards the estimation of D _L using regime channel concept. Having known discharge of flow and silt factor of the riverbed material, hydraulic parameters of the channel can be determined theoretically, which in turn can be used to estimate D _L appropriately making use of a suitable empirical formula.     

The Correlation Bias for Two-Dimensional Simulations by Turning Bands

The turning bands method (TBM) generates realizations of isotropic Gaussian random fields by summing contributions from line processes. We consider two-dimensional simulations and study the correlation bias attributable to the use of only a finite number L of lines. Our analytical and numerical results confirm that the maximal bias is of order 1/L, and that L = 64 lines suffice for excellent covariance reproduction. The notorious banding observed in simulations with an insufficient number of lines is a related but different phenomenon and depends strongly on the choice of the line simulation technique. Clear-cut recommendations for the number of lines necessary to avoid the effect can only be based on practical experience with the specific code at hand.     

Reverse flow in turbidity currents: the role of internal solitons

Pickering & Hiscott, (1985) have demonstrated amply the presence of reverse-flow units within the thick-bedded calcareous wacke (TCW) beds of the turbiditic Cloridorme Formation (Middle Ordovician, Gaspé Peninsula, Quebec, Canada). These reverse-flow units are underlain and overlain by units which reveal flow in the primary (obverse) direction. In this paper, a model is proposed for this reverse flow, based on the probable nature of the primary turbidity flow. It appears that the initial flow was highly elongated (thickness h? length L), with h～ 500 m, velocity U～ 2 m s^-1 and sediment concentration C～ 1·25%_o. The rate of momentum loss of the flow is estimated by means of a useful parameter which we call the ‘drag distance’, symbol d_D, defined by where h and L are the thickness and length of the flow, respectively; _cC_d is a combined drag coefficient representing friction on the bottom and at the upper interface; and _fC_d is a form-drag coefficient related to the shape and size of the head. d_D is the distance travelled by a current of constant h and L, flowing over a horizontal bottom and obeying a quadratic friction law, for an e-fold reduction in velocity. Simple considerations, confirmed by our own experiments (described in this paper), show that such an elongated turbidity current cannot be reflected as a whole from an adverse slope: when the nose of the current reaches the slope, it forms a hump, which surges backwards and sooner or later breaks up into a series of internal solitons. The latter, probably numbering 4–7, will cause reverse flow at a given point as they pass by, provided that the residual velocity in the tail is not too great. Flow in the original (obverse) direction will be re-established after the passage of the solitons. Quiescent periods in front of, between and behind the solitons, when soliton-associated currents cancelled out the residual obverse flow, would allow the deposition of thin mud-drapes. Additional flow reversals observed in a few of the TCW beds cannot be explained readily by the re-passage of solitons, since wave breaking at the ends of the basin would cause massive energy loss; internal seiches are the preferred explanation for these later reversals.