Analytical and numerical analysis of pressure drawdown in a poroelastic reservoir with complete overburden effect considered

Starting from an analytical reservoir model that incorporates full interaction with an elastic overburden, a new hybrid mathematical approach is developed by combining two numerical discretization methods. A tabular reservoir (petroleum reservoir or an aquifer) in an infinite or semi-infinite domain is viewed as a macroscopic displacement discontinuity, allowing use of the efficient displacement discontinuity mathematical method to calculate stresses and displacements that arise because of pressure changes. A 3-D finite element method using a poroelastic formulation is used to discretize the reservoir itself. By coupling the displacement discontinuity and finite element methods, a 3-D large-scale poroelastic reservoir can be simulated within an infinite or semi-infinite domain. The numerical model has been verified through comparison to known solutions, and some time-dependent pressure drawdown problems are analyzed. Results indicate that including the complete overburden (reservoir surroundings) response has a significant effect on pressure drawdown in a poroelastic reservoir during pumping, and should be incorporated in appropriate applications such as well test equations and subsidence analyses.     

数字岩心逆建模理论下的储层参数定量预测方法

数字岩心微观孔隙结构十分复杂,有限元模拟物性参数与弹性参数之间关系是非线性的,直接反演其物性参数准确度低、稳定性差.本文发展了一种数字岩石物理逆建模方法,实现了基于数字岩心的储层参数有效预测.从数字岩心基函数的构建出发,基于有限元方法,计算了一系列具有等间距物性参数值(孔隙度、泥质含量和含水饱和度)的数字岩心弹性参数(体积模量、剪切模量和密度),通过插值算法建立了数字岩心弹性参数三维数据集,从而实现了弹性模量的有限元数值解的快速构建;然后搜索弹性参数的单值等值面,通过等值面的空间交会得到交点,完成储层参数预测.测试结果表明:基于数字岩心逆建模理论的储层参数预测结果与实际模型一致,具有可行性,并且可以通过增加插值点数目提高预测的准确性;孔隙度和泥质含量预测结果稳定性很好,而含水饱和度对噪声的加入较为敏感.     

Asymmetry in the time-lapse seismic response to injection and depletion

Large changes in seismic reflection amplitude have been observed around injectors, and result from the decrease in elastic‐wave velocity due to the increase in pore pressure in the reservoir. In contrast, the velocity change resulting from the decrease in pore pressure in depleting reservoirs is observed to be smaller in magnitude. Elastic‐wave velocities in sandstones vary with stress due to the presence of stress‐sensitive grain boundaries within the rock. Grain‐boundary stiffness increases non‐linearly with increasing compressive stress, due to increased contact between opposing faces of the boundary. This results in a change in velocity due to a decrease in pore pressure that is smaller than the change in velocity caused by an increase in pore pressure, in agreement with time‐lapse seismic observations. The decrease in porosity resulting from depletion is not fully recovered upon re‐pressurization, and this leads to an additional steepening of the velocity vs. effective stress curve for injection relative to depletion. This difference is enhanced by any breakage of cement or weakening of grain contacts that may occur during depletion and by the reopening or formation of fractures or joints and dilation of grain boundaries that may occur during injection.     

Relating 4D seismics to reservoir geomechanical changes using a discrete element approach

A modified discrete element method is briefly introduced and used for modelling reservoir geomechanical response during fluid injection and depletion. The modified approach works as a continuum method until some local failure is initiated, after which it behaves like a discrete element method on a polygonal lattice. The method is advantageous for modelling fracture developments in rocks. It is applied here to synthetic models of two reservoirs taken from the North Sea (Gullfaks and Elgin‐Franklin). For Gullfaks, two cases of water injection were modelled, one with low horizontal effective stress and the other with low vertical effective stress. Vertical fractures are developed in the first case, whereas horizontal fractures are developed in the second case. This would not have been seen using traditional methods. Based on 4D seismics data for the Gullfaks field, one may envision that horizontal fractures could have been formed. The Elgin‐Franklin synthetic model is used to study various scenarios of changing stress field around the depleting reservoir. Based on 4D seismics data from this field, one may see changes that could be interpreted in terms of possible fault reactivation.     

Gravity Changes and Internal Processes: Some Results Obtained from Observations at Three Volcanoes

Temporal gravity changes provide information about mass and/or density variations within and below the volcano edifice. Three active volcanoes have been under investigation; each of them related to a plate boundary: Mayon/Luzon/Philippines, Merapi/Java/Indonesia, and Galeras/Colombia. The observed gravity changes are smaller than previously expected but significant. For the three volcanoes under investigation, and within the observation period, mainly the increase of gravity is observed, ranging from 1,000 nm^–2 to 1,600 nms^–2. Unexpectedly, the gravity increase is confined to a rather small area with radii of 5 to 8 km around the summit. At Mayon and Merapi the parallel GPS measurements yield no significant elevation changes. This is crucial for the interpretation, as the internal pressure variations do not lead to significant deformation at the surface. Thus the classical Mogi-model for a shallow extending magma reservoir cannot apply. To confine the possible models, the attraction due to changes of groundwater level or soil moisture is estimated along the slope of Merapi exemplarily by 2-D modelling. Mass redistribution or density changes were evaluated within the vent as well as deeper fluid processes to explain the gravity variations; the results are compared to the model incorporating the additional effect of elastic deformation.     

Topographic accuracy assessment of bare earth lidar-derived unstructured meshes

This study is focused on the integration of bare earth lidar (Light Detection and Ranging) data into unstructured (triangular) finite element meshes and the implications on simulating storm surge inundation using a shallow water equations model. A methodology is developed to compute root mean square error (RMSE) and the 95th percentile of vertical elevation errors using four different interpolation methods (linear, inverse distance weighted, natural neighbor, and cell averaging) to resample bare earth lidar and lidar-derived digital elevation models (DEMs) onto unstructured meshes at different resolutions. The results are consolidated into a table of optimal interpolation methods that minimize the vertical elevation error of an unstructured mesh for a given mesh node density. The cell area averaging method performed most accurate when DEM grid cells within 0.25 times the ratio of local element size and DEM cell size were averaged. The methodology is applied to simulate inundation extent and maximum water levels in southern Mississippi due to Hurricane Katrina, which illustrates that local changes in topography such as adjusting element size and interpolation method drastically alter simulated storm surge locally and non-locally. The methods and results presented have utility and implications to any modeling application that uses bare earth lidar.     

The scattering of plane P_1 wave by a 3-D sedimentary basin in a fluid-saturated poroelastic half-space using IBEM

The indirect boundary element method(IBEM) is applied to investigate the scattering of elastic waves around a 3-D sedimentary basin filled with fluidsaturated poroelastic medium. Based on this method, the free field and scattered field can be solved according to the boundary conditions. And the numerical accuracy has been verified. The effects of parameters on elastic wave scattering are studied, such as boundary condition, incident frequency,incident angle and porosity of medium. Numerical results illustrate that the amplification effect of surface displacement near poroelastic sedimentary basin is notable. In addition, for the case of large porosity the drainage condition has a significant impact on the response amplitude. Due to the fluid exchange at the interface under the drained condition, the displacement amplitude can be much larger than that under the undrained condition in present study. The study can provide a theoretical basis for the anti-seismic design of engineering structures located in sedimentary basin.     

水库触发地震研究中二维与三维有限元模拟结果比较

紫坪铺水库蓄水载荷引起的库仑应力变化以及是否可以触发汶川M_S8.0地震是地学界关注的一个热点问题。 已有较多研究均建立在二维有限元模拟基础上, 然而二维与三维结果存在较大差异, 本文定量给出了该差异的大小。 首先比较了有限元方法与解析方法求解半无限空间点载荷下弹性问题结果, 确认了三维有限元计算在百万网格下具有足够的计算精度; 然后数值计算并比较了二维、 三维模型中紫坪埔水库蓄水的弹性力学及渗流效应。 结果表明, 总体上弹性加载会造成逆断层更趋安全, 而水压渗透使其趋于危险, 两者叠加仍然没有导致增加汶川地震危险性; 但现有的二维有限元模拟夸大了真实结果, 在汶川震源位置三维计算的应力大小仅为二维计算的1/3到1/4, 因此在严格的讨论中必须使用三维建模。     

Elastic properties of double-porosity rocks using the differential effective medium model

An approach to determining the effective elastic moduli of rocks with double porosity is presented. The double‐porosity medium is considered to be a heterogeneous material composed of a homogeneous matrix with primary pores and inclusions that represent secondary pores. Fluid flows in the primary‐pore system and between primary and secondary pores are neglected because of the low permeability of the primary porosity. The prediction of the effective elastic moduli consists of two steps. Firstly, we calculate the effective elastic properties of the matrix with the primary small‐scale pores (matrix homogenization). The porous matrix is then treated as a homogeneous isotropic host in which the large‐scale secondary pores are embedded. To calculate the effective elastic moduli at each step, we use the differential effective medium (DEM) approach. The constituents of this composite medium – primary pores and secondary pores – are approximated by ellipsoidal or spheroidal inclusions with corresponding aspect ratios. We have applied this technique in order to compute the effective elastic properties for a model with randomly orientated inclusions (an isotropic medium) and aligned inclusions (a transversely isotropic medium). Using the special tensor basis, the solution of the one‐particle problem with transversely isotropic host was obtained in explicit form. The direct application of the DEM method for fluid‐saturated pores does not account for fluid displacement in pore systems, and corresponds to a model with isolated pores or the high‐frequency range of acoustic waves. For the interconnected secondary pores, we have calculated the elastic moduli for the dry inclusions and then applied Gassmann's tensor relationships. The simulation of the effective elastic characteristic demonstrated that the fluid flow between the connected secondary pores has a significant influence only in porous rocks containing cracks (flattened ellipsoids). For pore shapes that are close to spherical, the relative difference between the elastic velocities determined by the DEM method and by the DEM method with Gassmann's corrections does not exceed 2%. Examples of the calculation of elastic moduli for water‐saturated dolomite with both isolated and interconnected secondary pores are presented. The simulations were verified by comparison with published experimental data.     

Parallel 3-D Simulation of a Fault Gouge Using the Lattice Solid Model

Despite the insight gained from 2-D particle models, and given that the dynamics of crustal faults occur in 3-D space, the question remains, how do the 3-D fault gouge dynamics differ from those in 2-D? Traditionally, 2-D modeling has been preferred over 3-D simulations because of the computational cost of solving 3-D problems. However, modern high performance computing architectures, combined with a parallel implementation of the Lattice Solid Model (LSM), provide the opportunity to explore 3-D fault micro-mechanics and to advance understanding of effective constitutive relations of fault gouge layers. In this paper, macroscopic friction values from 2-D and 3-D LSM simulations, performed on an SGI Altix 3700 super-cluster, are compared. Two rectangular elastic blocks of bonded particles, with a rough fault plane and separated by a region of randomly sized non-bonded gouge particles, are sheared in opposite directions by normally-loaded driving plates. The results demonstrate that the gouge particles in the 3-D models undergo significant out-of-plane motion during shear. The 3-D models also exhibit a higher mean macroscopic friction than the 2-D models for varying values of interparticle friction. 2-D LSM gouge models have previously been shown to exhibit accelerating energy release in simulated earthquake cycles, supporting the Critical Point hypothesis. The 3-D models are shown to also display accelerating energy release, and good fits of power law time-to-failure functions to the cumulative energy release are obtained.     

邢台强震区的深部构造对强震孕育影响的三维数值模拟

为了研究地壳结构特征与强震孕育的关系,运用三维有限单元法,计算了邢台地震区壳内应力场的扰动状态。结果表明:地壳内深、浅断裂的存在,造成沿直立断层带的应力相对集中和不同部位的变化。在其上部出现平均应力和水平剪应力的增加,中下部出现平均应力的减小和水平剪应力的增加。前者有利于弹性位能的高度集中和主破裂的发生,后者有利于主破裂前蠕滑的发生。地壳内高速体、低速体、深断裂、莫霍面隆起是孕育强震的主要构造因素,但影响程度,范围、形式不同,高速块体是应力增强的体;直立断层带上及其邻近应力集中程度最高,应力变化梯度最大,是产生震前预滑及主余震发生的有利部位;莫霍面隆起加剧高速体内的应力增强,加大断层带上的应力变化梯度而有利于地震的孕育和发生。     

Breaking Up: Comminution Mechanisms in Sheared Simulated Fault Gouge

The microstructural state and evolution of fault gouge has important implications for the mechanical behaviour, and hence the seismic slip potential of faults. We use 3D discrete element (DEM) simulations to investigate the fragmentation processes operating in fault gouge during shear. Our granular fault gouge models consist of aggregate grains, each composed of several thousand spherical particles stuck together with breakable elastic bonds. The aggregate grains are confined between two blocks of solid material and sheared under a given normal stress. During shear, the grains can fragment in a somewhat realistic way leading to an evolution of grain size, grain shape and overall texture. The ‘breaking up’ of the fault gouge is driven by two distinct comminution mechanisms: grain abrasion and grain splitting. The relative importance of the two mechanisms depends on applied normal stress, boundary wall roughness and accumulated shear strain. If normal stress is sufficiently high, grain splitting contributes significantly to comminution, particularly in the initial stages of the simulations. In contrast, grain abrasion is the dominant mechanism operating in simulations carried out at lower normal stress and is also the main fragmentation mechanism during the later stages of all simulations. Rough boundaries promote relatively more grain splitting whereas smooth boundaries favor grain abrasion. Grain splitting (plus accompanying abrasion) appears to be an efficient mechanism for reducing the mean grain size of the gouge debris and leads rapidly to a power law size distribution with an exponent that increases with strain. Grain abrasion (acting alone) is an effective way to generate excess fine grains and leads to a bimodal distribution of grain sizes. We suggest that these two distinct mechanisms would operate at different stages of a fault’s history. The resulting distributions in grain size and grain shape may significantly affect frictional strength and stability. Our results therefore have implications for the earthquake potential of seismically active faults with accumulations of gouge. They may also be relevant to the susceptibility of rockslides since non-cohesive basal shear zones will evolve in a similar way and potentially control the dynamics of the slide.     

三维有限元数值模拟在宁夏大坝电厂地震安全性评价中的应用

本文综合考虑了宁夏大坝电厂周围２５０ｋｍ范围内的地质、地球物理方面的有关资料，建立了三维有限元模型。利用该模型模拟了该区域的构造应力场，进而对大坝电厂厂址区及附近进行了地震危险性分析。在分析工作中结合安全度的定义提出了危险系数的概念及潜在震源区划分的定量－定性标志。经过分析，作者认为大坝电厂附近构造稳定性较好。     

利用三维孔隙弹性模型探讨紫坪铺水库对汶川地震的影响

本文基于三维孔隙弹性理论,建立了紫坪铺水库及周边地区的有限元模型.根据紫坪铺水库开始蓄水到汶川地震发震时刻的水位变化情况,计算了整个区域的孔隙压力和库仑应力.详细讨论了断层及周围地层的弹性模量和扩散系数对计算结果的影响.计算结果表明:从弹性角度看,断层的弹性模量对汶川地震震源处的库仑应力影响很小;震源处的库仑应力随着断层和周围地层的扩散系数增大而增大.当给定弹性模量和扩散系数代表性值的时候,计算结果表明在汶川地震发震时刻,震源处的库仑应力变化量为+1 kPa左右,这表明紫坪铺水库使得汶川地震发震断层更加危险.是否这个量级的库仑应力就能够触发汶川地震还需要进一步探讨.通过分析库区周边小震的分布,发现小震分布区域均是库仑应力增加的地区,因此紫坪铺水库周边的小震应该与紫坪铺水库蓄水有直接关系.     

弹性随机介质模型的特征频率

为研究弹性随机介质模型中的波场特征,本文使用弹性波动方程正演模拟了平面波在二维弹性随机介质模型中的传播.通过大量正演模拟,我们首先发现若使用不同频率的瑞克子波作为震源函数,计算散射波场的能量相对值ΔE,则对每一个固定的随机介质模型,ΔE都会在某一个与模型对应的震源频率f*处达到最大值.本文由此提出了随机介质模型的特征频率这一新的概念.本文充分说明了特征频率的客观存在性;给出了相应的计算方法;进一步全面研究了随机介质中的各种模型特征（如自相关长度、背景速度、扰动标准差以及模型尺寸等）与模型的特征频率之间的关系;并得到了若干结论和相应的经验关系表达式.     

水库蓄放水对库底岩石介质弹性波速影响的数值模拟

为了解水库水体加、卸载及渗透过程中,库底岩石介质有效应力变化对弹性波速的影响,以及弹性波速变化与孔隙压力、有效应力变化之间的定量关系,以四川省紫坪铺水库为例,在前期已建立库区地质构造与水文地质结构模型,并推导出岩石介质变形与流体渗流耦合数学模型的基础上,进一步建立了有效应力对岩石弹性波速影响的数学模型,利用有限元方法定量计算了库底岩石介质的孔隙压力、附加有效应力和弹性波速的变化。结果表明:在水库水体加、卸载及渗透过程中,地下岩石介质弹性波速增量值的变化并不明显,最大变幅仅为±0.013km/s,变化范围主要集中在库底5km以上的局部区域;弹性波速增量值与附加有效应力变化同步,与水库水位变化形态相似;不同观测点弹性波速变化量的大小与埋藏位置、深度及附加有效应力的变化幅度有关,同一观测点,弹性纵波波速变化幅度大于弹性横波波速。文中数值模拟的结果与一些学者通过环境噪声方法获得的紫坪铺水库附近区域相对波速变化的反演结果在变化形态上大体一致,但变化幅度略有差异。     

卡里巴水库蓄水引起库区应力场变化影响分析

水库地震是近几年国际上地震学和地球动力学领域研究的热点问题。 作为世界上库容最大的卡里巴(Kariba)水库在其蓄水后库区地震活动性不断增加, 并于1963年9月23日发生M6.1地震, 是世界上公认发生过6级以上水库地震的四大水库之一。 对该水库蓄水引起库区应力场变化进行量化估算, 有助于对特大—大型水库蓄水后库区应力场变化及应力触发地震这一基础性科学问题的深入了解。 本文分别应用解析解和数值解方法对卡里巴水库蓄水引起库区弹性应力场、 不排水和排水孔隙压力变化进行了计算分析。 结果表明, 卡里巴水库蓄水引起M6.1地震震源处孔压变化为0.015~0.22 MPa, 库仑应力变化约0.03~0.17 MPa, 触发了此次地震的发生。     

三维竖向断层场地对地面运动的影响

实际工程场地中的断层通常具有三维尺寸,应进行三维分析。本文考虑均匀弹性半空间内存在不同长度、深度的三维断层,由下向上垂直入射一个脉冲波,主要采用三维时域显式有限元结合透射人工边界的方法及相应编写的三维显示有限元波动程序DSI3,计算断层附近地震地面运动的时程反应,研究含有三维断层场地地面运动的一些特点,并与自由场反应对比,以说明三维竖向断层对地震波传播的特点和对断层附近地面运动的影响。     

弹性波边界元法正演模拟

弹性波边界元地震模型方法（ＢＥＥＳＭ），实现了二维和三维问题的纵、横波及转换波的同时模拟，并且能模拟任意复杂构造的地震声波正演模型．根据地震模型的特点，本文发展了数值积分计算与矩阵消元同步进行的块状高斯消元法；用解析法处理奇异积分；用无限元法处理边界吸收问题；采用单元长度随介质速度和计算频率变化的变单元算法，及自动剖分单元等技术，提高了计算精度，节省了内存，缩短了计算时间．     

地震低频伴影的数值模拟与应用

Strong low-frequency energy beneath a hydrocarbon reservoir is called a seismic low-frequency shadow and can be used as a hydrocarbon indicator （Tarter et al., 1979） bu the physical mechanism of the observed low-frequency shadow is still unclear. To stud） the mechanism, we performed seismic numerical simulation of geological models with a hydrocarbon-bearing zone using the 2-D diffusive-viscous wave equation which car effectively model the characteristics of velocity dispersion and transform the seismic dat~ centered in a target layer slice within a time window to the time-frequency domain by usinl time-frequency signal analysis and sort the frequency gathers to common frequency cubes. Then, we observe the characteristics of the seismic low-frequency shadow in the common frequency cubes. The numerical simulations reveal that the main mechanism of seismic lowfrequency shadows is attributed to high attenuation of the medium to high seismic frequency components caused by absorption in the hydrocarbon-filled reservoir. Results from a practical example of seismic low-frequency shadows show that it is possible to identify the reservoir by the low-frequency shadow with high S/N seismic data.