Modeling nonlinear response of fractured rocks and reservoirs

Since the early work of Athey (1930), there have been many attempts to describe the various nonlinear behaviors of rocks and soils in terms of functionals having only a few parameters, while nevertheless being able to fit the complicated available data with satisfactory accuracy. Such approaches have not been universally applied however, and the present analyses are intended to draw attention to the possibility of using such nonlinear fitting methods on old as well as new data sets. In particular, some special emphasis is placed here on re‐examining the well‐known laboratory data of Coyner (1984) on rocks in light of such modeling tools, and we find that the nonlinear approach again has several clear advantages – especially in terms of reducing the number of variables needed to describe the observed behavior of both bulk modulus and porosity of rocks undergoing large changes in pressure. Copyright © 2016 John Wiley & Sons, Ltd.     

Modeling of subsidence and stress-dependent hydraulic conductivity for intact and fractured porous media

Summary This study investigates the changes in deformation and stress dependent hydraulic conductivities that occur as a result of underground mining in intact and fractured porous media. The intact porous medium is assumed to be comprised of regularly packed spherical grains of uniform size. The variation in grain size or pore space due to the effect of changing intergranular stresses results in a change in rock hydraulic conductivity. A model is developed to describe the sensitivity of hydraulic conductivity to effective stresses through Hertzian contact of spherical grains. The fractured porous medium is approximated as an equivalent fracture network in which a single fracture is idealized as a planar opening having a constant equivalent thickness or aperture. Changes in fracture aperture as a result of changes in elastic deformation control the variation of hydraulic conductivity. A model is presented to illustrate the coupling between strain and hydraulic conductivity. Subsidence induced deformations that result from mining induced changes in hydraulic conductivity in both intact and fractured media. These changes are examined and compared with results from a mining case study.     

Simulation of multiphase flow in fractured reservoirs using a fracture-only model with transfer functions

We present a fracture-only reservoir simulator for multiphase flow: the fracture geometry is modeled explicitly, while fluid movement between fracture and matrix is accommodated using empirical transfer functions. This is a hybrid between discrete fracture discrete matrix modeling where both the fracture and matrix are gridded and dual-porosity or dual-permeability simulation where both fracture and matrix continua are upscaled. The advantage of this approach is that the complex fracture geometry that controls the main flow paths is retained. The use of transfer functions, however, simplifies meshing and makes the simulation method considerably more efficient than discrete fracture discrete matrix models. The transfer functions accommodate capillary- and gravity-mediated flow between fracture and matrix and have been shown to be accurate for simple fracture geometries, capturing both the early- and late-time average behavior. We verify our simulator by comparing its predictions with simulation results where the fracture and matrix are explicitly modeled. We then show the utility of the approach by simulating multiphase flow in a geologically realistic fracture network. Waterflooding runs reveal the fraction of the fracture–matrix interface area that is infiltrated by water so that matrix imbibition can occur. The evolving fraction of the fracture–matrix interface area turns out to be an important characteristic of any particular fracture system to be used as a scaling parameter for capillary driven fracture–matrix transfer.     

3-D geometric description of fractured reservoirs

Traditional stochastic modeling of fracture networks usually failed because it required unaccessible statistics and may not be able to honor available local data. This paper presents an algorithm for the 3D geometric simulation of fractured reservoirs. It is based on geological rules of fracture propagation and interaction. It is part of a methodology which aims at integrating diverse data about the fracture system in the subsurface. This information can come from well cores and logs, analog outcrops, geomechanical stress studies, seismic surveys; it may be quantitative or qualitative, and have different degrees of reliability.     

An evolutionary system for exploitation of fractured geothermal reservoirs

In this work, we present a novel methodology to integrate one of the most advanced technique for modeling fractured media for underground systems with a semantics-based genetic programming technique. The objective of the study is to develop a global framework to forecast the temperature of fractured reservoirs. The numerical method used to solve the physical equations is able to handle different fracture distributions without changing the background computational grid, i.e. the mesh of the rock matrix, as well as letting geometrically uncoupled the one co-dimensional fracture meshes. In the context of temperature forecasting, the use of a recently defined variant of genetic programming is taken into account for finding (quasi-)perfect solutions with high probability and for generating models able to produce near optimal predictions also on unseen data. The proposed computational intelligence technique integrates, in a recently developed version of genetic programming that uses semantic genetic operators, a “greedy” crossover and a self tuning algorithm. Experimental results confirm the suitability of the proposed method in predicting the correct temperature distribution in probes inside the domain.     

Heterogeneity preserving upscaling for heat transport in fractured geothermal reservoirs

In simulation of fluid injection in fractured geothermal reservoirs, the characteristics of the physical processes are severely affected by the local occurence of connected fractures. To resolve these structurally dominated processes, there is a need to develop discretization strategies that also limit computational effort. In this paper, we present an upscaling methodology for geothermal heat transport with fractures represented explicitly in the computational grid. The heat transport is modeled by an advection-conduction equation for the temperature, and solved on a highly irregular coarse grid that preserves the fracture heterogeneity. The upscaling is based on different strategies for the advective term and the conductive term. The coarse scale advective term is constructed from sums of fine scale fluxes, whereas the coarse scale conductive term is constructed based on numerically computed basis functions. The method naturally incorporates the coupling between solution variables in the matrix and in the fractures, respectively, via the discretization. In this way, explicit transfer terms that couple fracture and matrix solution variables are avoided. Numerical results show that the upscaling methodology performs well, in particular for large upscaling ratios, and that it is applicable also to highly complex fracture networks.     

水平井压裂储层应力场分布模拟方法

在研究分析水力压裂对储层岩石力学特性参数影响的基础上,提出一种压力储层应力场分布模拟计算方法。通过建立水平井储层原地应力场模型和水力压裂产生人工裂缝诱导应力场模型,并且利用实际的水力压裂测井参数对储层原地应力场和压裂产生裂缝诱导应力场分布进行了模拟计算。模拟计算结果表明,压裂产生人工裂缝会对储层应力场分布造成很大影响;压裂后储层应力主要在裂缝周围得到积累,并且距离裂缝越远,应力值积累越少;压裂生成裂缝长度也会影响储层应力场分布,裂缝越长,裂缝诱导应力场减小越慢。     

Boundary element analysis of contaminant transport in fractured porous media

A two-dimensional boundary integral method to analyse the flow of contaminant in fractured media having a two- or three-dimensional orthogonal fracture network is presented. The method assumes that the fractures provide the paths of least resistance for transport of contaminants while the matrix, because of its low permeability, acts as ‘storage blocks’ into which the contaminant diffuses. Laplace transform is used to eliminate the time variable in the governing equation in order to facilitate the formulation of a boundary integral equation in the Laplace transform space. Conventional boundary element techniques are applied to solve for the contaminant concentrations at specified locations in the spatial domain. The concentration in the time domain is then obtained by using an efficient inversion technique developed by Talbot. The method is able to analyse the behaviour of waste repositories which have diminishing concentration due to the mass transport of the contaminant into the surrounding fractured media.     

Ensemble clustering for efficient robust optimization of naturally fractured reservoirs

In the development of naturally fractured reservoirs (NFRs), the existence of natural fractures induces severe fingering and breakthrough. To manage the flooding process and improve the ultimate recovery, we propose a numerical workflow to generate optimal production schedules for smart wells, in which the inflow control valve (ICV) settings can be controlled individually. To properly consider the uncertainty introduced by randomly distributed natural fractures, the robust optimization would require a large ensemble size and it would be computationally demanding. In this work, a hierarchical clustering method is proposed to select representative models for the robust optimization in order to avoid redundant simulation runs and improve the efficiency of the robust optimization. By reducing the full ensemble of models into a small subset ensemble, the efficiency of the robust optimization algorithm is significantly improved. The robust optimization is performed using the StoSAG scheme to find the optimal well controls that maximize the net-present-value (NPV) of the NFR’s development. Due to the discrete property of a natural fracture field, traditional feature extraction methods such as model-parameter-based clustering may not be directly applicable. Therefore, two different kinds of clustering-based optimization methods, a state-based (e.g., s _w profiles) clustering and a response-based (e.g., production rates) clustering, are proposed and compared. The computational results show that the robust clustering optimization could increase the computational efficiency significantly without sacrificing much expected NPV of the robust optimization. Moreover, the performance of different clustering algorithms varies widely in correspondence to different selections of clustering features. By properly extracting model features, the clustered subset could adequately represent the uncertainty of the full ensemble.     

Shear behaviour of fractured rock as a function of size and shear displacement

This paper presents laboratory results regarding the shear behaviour of an artificial tensile fracture generated in granite. We used a direct shear rig to test fractures of different sizes (from 100 mm to 200 mm) under various shear displacements up to 20 mm and cyclic shear stresses with constant normal stress of 10 MPa. To determine the evolution of surface damage and aperture during shear, cyclic loading was performed at designated shear displacements. These changes in the surfaces topography were measured with a laser profilometer ‘non-contact surface profile measurement system’. In addition, changes were also measured directly by using pressure-sensitive film.

The results showed that the standard deviation (SD) of the initial aperture of the sheared fracture significantly increases with both shear displacement and size, which result in an increase in the non-linearity of the closure curve (since the matedness of the fracture surfaces decreases with shear displacement). Therefore, we concluded that shear dilation is not only governed by the surfaces sliding over each other, but is also strongly influenced by the non-linearity of closure with shear displacement. Furthermore, while the shear stiffness of the fracture during the initial stage decreases with fracture size, it increases with fracture size in the residual stage. This can be attributed to the fact that only small asperities with short wavelengths were mainly damaged by shearing. Moreover the result showed that the damaged zones enlarge and localise with shear displacement, and eventually tend to form perpendicular to the shear displacement.     

Double-porosity modelling of oscillatory gas motion and contaminant transport in a fractured porous medium

A double-porosity model is used to describe the oscillatory gas motion and associated contaminant transport induced by cyclical variations in the barometric pressure at the surface of a fractured porous medium. Flow along the fractures and within the permeable matrix blocks is locally one-dimensional. The interaction between fractures and blocks includes seepage of fluid as well as diffusion of contaminant. To guard against artificial numerical diffusion, the FRAM filtering remedy and methodology of Chapman is used in calculating the advective fluxes along fractures and within blocks. The entire system of equations, including the fracture-matrix interaction terms, is solved by a largely implicit non-iterative algorithm which remains stable and conservative even when the computational time step is large compared to the cross-block transit time of pressure waves. The numerical accuracy is tested by comparison with exact solutions for oscillatory and unidirectional flows, some of which include diffusion interaction between the fracture and the matrix. The method is used to estimate the rate of vertical transport of radioactive gases through the rubblized chimney produced by an underground nuclear explosion.     

Upscaling permeability for three-dimensional fractured porous rocks with the multiple boundary method

Upscaling permeability of grid blocks is crucial for groundwater models. A novel upscaling method for three-dimensional fractured porous rocks is presented. The objective of the study was to compare this method with the commonly used Oda upscaling method and the volume averaging method. First, the multiple boundary method and its computational framework were defined for three-dimensional stochastic fracture networks. Then, the different upscaling methods were compared for a set of rotated fractures, for tortuous fractures, and for two discrete fracture networks. The results computed by the multiple boundary method are comparable with those of the other two methods and fit best the analytical solution for a set of rotated fractures. The errors in flow rate of the equivalent fracture model decrease when using the multiple boundary method. Furthermore, the errors of the equivalent fracture models increase from well-connected fracture networks to poorly connected ones. Finally, the diagonal components of the equivalent permeability tensors tend to follow a normal or log-normal distribution for the well-connected fracture network model with infinite fracture size. By contrast, they exhibit a power-law distribution for the poorly connected fracture network with multiple scale fractures. The study demonstrates the accuracy and the flexibility of the multiple boundary upscaling concept. This makes it attractive for being incorporated into any existing flow-based upscaling procedures, which helps in reducing the uncertainty of groundwater models.     

Modeling spatial fracture intensity as a control on flow in fractured rock

Spatial fracture intensity (P ₃₂, fracture area by volume) is an important characteristic of a jointed rock mass. Although it can hardly ever be measured, P ₃₂ can be modeled based on available geological information such as spatial data of the fracture network. Flow in a mass composed of low-permeability hard rock is controlled by joints and fractures. In this article, models were developed from a geological data set of fractured andesite in LanYu Island (Taiwan) where a site is investigated for possible disposal of low-level and intermediate-level radionuclide waste. Three different types of conceptual models of spatial fracture intensity distribution were generated, an Enhanced Baecher’s model (EBM), a Levy–Lee Fractal model (LLFM) and a Nearest Neighborhood model (NNM). Modeling was conducted on a 10 × 10 × 10 m synthetic fractured block. Simulated flow was forced by a 1% hydraulic gradient between two vertical x–z faces of the cube (from North to South) with other boundaries set to no-flow conditions. Resulting flow vectors are very sensitive to spatial fracture intensity (P ₃₂). Flow velocity increases with higher fracture intensity (P ₃₂). R-squared values of regression analysis for the variables velocity (V/V _max) and fracture intensity (P ₃₂) are 0.293, 0.353, and 0.408 in linear fit and 0.028, 0.08, and 0.084 in power fit. Higher R ² values are positively linked with structural features but the relation between velocity and fracture intensity is non-linear. Possible flow channels are identified by stream-traces in the Levy–LeeFractal model.     

断裂岩石蠕剪中的细观接触损伤试验研究

为深刻理解构成断裂岩体长期抗剪强度的细观机制,对岩石断裂面的细观接触和损伤演化进行了试验研究。采用压剪方式和巴西劈裂方式制作出两类断裂岩石,在恒定法向力作用下对破坏岩石试件进行分级施加剪切力的蠕变试验。在加载前、中、后对断裂岩石分别进行CT和激光扫描,观察到不同蠕变阶段断裂岩石的细观接触和损伤状态,获得不同性质的断裂岩石抗剪强度特征。试验研究表明：构成断裂岩石长期抗剪强度的机制主要有两个：一是细观凹凸啮合体的抗剪断能力;二是表观凹凸接触体的抗摩擦能力。拉破裂岩石表面局部粗糙度相对较大,抗剪强度以第1种破坏机制为主,剪切破坏岩石表面宏观起伏度较大,抗剪强度是以第2种破坏机制为主。在蠕变剪切过程中,两种机制交织在一起,并随时间或剪位移的增加而相互转换。     

Semi-empirical equations for the systematic decrease in permeability with depth in porous and fractured media

Permeability loss with depth is a general trend in geological media and plays an essential role in subsurface fluid flow and solute transport. In the near surface zone where groundwater movement is active, the decrease in permeability with depth is dominated by the mechanical compaction of deformable media caused by the increase in lithostatic stress with depth. Instead of using empirical equations from statistical analysis, by considering the well-defined relationships among permeability, porosity, fracture aperture and effective stress under lithostatic conditions, new semi-empirical equations for the systematic depth-dependent permeability are derived, as well as the equations for the depth-dependent porosity in a porous medium and the depth-dependent fracture aperture in a fractured medium. The existing empirical equations can be included in the new equations as special cases under some simplification. These new semi-empirical equations perform better than previous equations to interpret the depth-dependent permeability of the Pierre Shale (with a maximum depth of approximately 4,500 m) and the granite at Stripa, Sweden (with a maximum depth of about 2,500 m).     

Architectural characteristics and petrophysical properties evolution of a strike-slip fault zone in a fractured porous carbonate reservoir

This paper describes the structural, petrophysical and hydromechanical properties relationships between a small fault zone and the porous layered carbonate series which host it. In a gallery located at 250-m depth, the deformation of a 22-m thick section of layered carbonates-, affected by a strike slip-fault have been characterized by means of structural (Q-value), acoustic velocities (V_p), porosity and uniaxial compressive strength (σ_c) measurements conducted in situ at the meter scale, and on laboratory samples at the infra-centimeter scale. A clear influence of the layers initial properties on fault architecture and properties evolution is underlined. In the porous layers with a low σ_c, there is an important accommodation of the deformation by micro-mechanisms resulting in a progressive decrease in the porosity toward the fault core. In the low-porosity layers with a high σ_c, deformations are accommodated toward the fault core by: an increase in the fracture porosity, in the micro-cracks porosity, and by displacements along pre-existing fractures resulting from a joint roughness decrease. The fault zone appears as relatively stiff and low permeable zones intercalated with low stiffness and high fracture permeability zones that extend one to tens of meters from the fault following the initial properties contrasts and geometry of the sedimentary layers.     

Performance of multi-stage fractured horizontal wells with stimulated reservoir volume in tight gas reservoirs considering anomalous diffusion

Horizontal well combined with volume fracturing technology has been extensively employed in the development of tight gas reservoirs. The disordered distribution of the induced and natural fractures in the reservoirs leads to the existence of the anomalous diffusion, so the conventional Darcy law has some limitations in describing the fluid flow under this circumstance. This paper introduces the fractional Darcy law to take into account the effect of the anomalous diffusion and then extends the conventional model of the multi-stage fractured horizontal (MSFH) well with the presence of the stimulated reservoir volume (SRV). The generated point source model for dual-porosity composite system includes the fractional calculus and its solution in Laplace space is derived. The superposition principle and the numerical discrete method are applied to obtain the solution for the MSFH well with SRV. Stehfest inversion method is used to transform the pseudo-pressure and production rate from Laplace space to real space. Type curves for pseudo-pressure and production rate are presented and analyzed. The influence of the relevant parameters on pseudo-pressure behavior and production rate decline is discussed in detail. The proposed model enriches the flow models of the MSFH well with SRV and can be used to more accurately interpret and forecast the transient pressure and transient rate.