Coupled Process Models of Fluid Flow and Heat Transfer in Hydrothermal Systems in Three Dimensions

Geothermal fields and hydrothermal mineral deposits are manifestations of the interaction between heat transfer and fluid flow in the Earth’s crust. Understanding the factors that drive fluid flow is essential for managing geothermal energy production and for understanding the genesis of hydrothermal mineral systems. We provide an overview of fluid flow drivers with a focus on flow driven by heat and hydraulic head. We show how numerical simulations can be used to compare the effect of different flow drivers on hydrothermal mineralisation. We explore the concepts of laminar flow in porous media (Darcy’s law) and the non-dimensional Rayleigh number (Ra) for free thermal convection in the context of fluid flow in hydrothermal systems in three dimensions. We compare models of free thermal convection to hydraulic head driven flow in relation to hydrothermal copper mineralisation at Mount Isa, Australia. Free thermal convection occurs if the permeability of the fault system results in Ra above the critical threshold, whereas a vertical head gradient results in an upward flow field.     

Coupled Process Models as a Tool for Analysing Hydrothermal Systems

Hydrothermal systems are characterised by complex interactions between heat transfer, fluid flow, deformation, species transport and chemical reactions. Numerical models can provide quantitatively constrained information in regions where acquisition of new data is difficult or expensive thus providing a means for reducing risks, costs, and effort during targeting, production, and management of resources linked to hydrothermal systems. Here we show how numerical simulations of hydrothermal processes can be used to better understand coupled reactive transport in modern geothermal systems and in ancient hydrothermal ore deposits. We give examples based on the Enhanced Geothermal System at Soultz-sous-Forêts in France, hydrothermal mineralisation at Mount Isa in Australia, and the geothermal resource at Hamburg-Allermöhe in Germany.     

洋壳渗透率结构对其内部热液对流的影响研究——基于热力学数值模拟

洋底特殊环境的限制使数值模拟成为研究海底热液对流与成矿机制的有效方法.本文通过数值模拟的方法,研究洋壳渗透率单因素变化对热液对流系统的形态和输出参数(热液喷发温度和热流值)进行研究.模拟结果表明,洋壳平均渗透率分别与热液喷发温度和热量输出呈反相关和正相关关系,符合达西流体热对流的基本解析规律.另外,洋壳渗透率的垂向变化不会使洋壳内部的对流系统产生明显的横向偏移,只使热液羽规模在垂向上产生规律性变化,且渗透率越大,热液羽规模越小.渗透率在水平方向的变化则是引发热液羽和热液喷发位置横向偏移的主要因素,且只有在特定的渗透率水平变化率范围内,上涌热液羽才会发生向高渗区的明显偏移,这也从一定程度上解释了现今拆离断层相关的热液系统均未发育于断层终止线上的事实.以大西洋Trans-Atlantic Geotraverse(TAG)热液区为例,断裂带高渗区必然影响相邻洋壳内部热液的上涌路径,但受区内拆离断裂带周缘的渗透率水平变化规律的影响,上涌热液羽不至于完全偏移至断裂带,从而产生拆离断层上盘的热液活动.     

Modelling Feedback of Chemical Reactions on Flow Fields in Hydrothermal Systems

Coupled reactive transport models of hydrothermal systems provide new insights and deeper understanding of the processes occurring due to fluid flow, heat transfer, solute transport, and chemical reactions. Basic concepts of species transport (diffusion, dispersion, and advection) and chemical precipitation and dissolution reactions are discussed, and five end-member types of reactive transport environments are introduced. One of these reactive transport environments, named ‘reactions within thermal gradients’, is used to demonstrate how free thermal convection can lead to redeposition of minerals and, due to the feedback of reaction on the flow field, a change of the convection pattern. The direct consequence of changing the flow field is a significant variation of the temperature distribution within the modelled area. With the example it is shown how reactive transport simulation can be applied for the detailed study of fossil and recent hydrothermal systems.     

Influence of permeability on hydrothermal circulation in the sediment-buried oceanic crust

1 Introduction Hydrothermal circulation is the key process of hydrothermal activity. Modern seafloor hydrothermal circulation can be divided into three parts: convective cells in the oceanic curst, interface between seafloor and ocean and hydrothermal plume. Hydrothermal convection in the crust is the dominant part of the whole seafloor hydrothermal circulation. The distribu-tion and nature of hydrothermal system in the oceanic crust are controlled by crust thermal structures and permeability …     

Modeling hydrothermal fluid circulation and gravity signals at the Phlegraean Fields caldera

The Phlegraean Fields caldera is an active volcanic system where episodes of ground deformation are accompanied by significant changes in geochemical and geophysical parameters monitored at the surface. These changes derive from a complex interaction between magmatic system and hydrothermal fluid circulation. We calculate the gravity changes associated with the variable density of hydrothermal fluids. We simulate the multi-phase and multi-component fluid circulation triggered by a pulsating magma degassing, periodically increasing the discharge of CO₂-enriched fluids into the shallow hydrothermal system. The simulated evolution of the hydrothermal system successfully reproduces the observed composition of gas discharged at the surface. At the same time, results indicate that changes in average fluid density generate a detectable gravity signal that is of the same order of magnitude of the observed changes. This contribution to gravity changes can explain the peculiar behavior of gravity data collected at Solfatara, where surface hydrothermal phenomena are present. Simultaneous fitting of two independent sets of monitoring data (gas composition and gravity changes) confirms the conceptual model proposed for the hydrothermal system at Solfatara, and it provides new insights for the interpretation of gravity data.     

Geothermal constraints on the hydrological regime of the TAG active hydrothermal mound, inferred from long-term monitoring

During August 1994 to March 1995, a period that included ODP Leg 158 drilling, bottom-water and sub-bottom temperatures were continuously logged by a long-term temperature monitoring system ‘Daibutsu’ at the base of the central black-smoker complex (CBC) and within the low heat flow zone at the TAG hydrothermal mound on the Mid-Atlantic Ridge. The temperature of hydrothermal fluid at CBC was also measured with a small high-temperature probe ‘Hobo’. Bottom-water temperature variations measured with Daibutsu at both sites have predominant semi-diurnal periods, causing the sub-bottom temperatures to fluctuate at these periods with reduced amplitudes and phase delays at sub-bottom depths. Seawater entrainment into the mound has been previously suggested at the low heat flow zone. We quantitatively evaluate the seawater entrainment rate at both sites from a one-dimensional numerical model, combined with a heat conduction model for the semi-diurnal variations. The entrainment rate of seawater at the base of CBC is estimated as 1.3±0.5×10⁻⁵ m/s, at least from August 17 to 30, 1994. On the other hand, the seawater entrainment rate at the low heat flow zone was undetected by long-term temperature monitoring at shallow sub-bottom depth. Nevertheless an increase in heat flow observed at the low heat flow zone during ODP drilling can be interpreted as a decrease in the entrainment rate of seawater. Before ODP Leg 158, Daibutsu measured three sub-bottom temperature anomalies at the base of CBC not derived from bottom-water temperature variations and Hobo also detected a CBC fluid temperature anomaly, indicating some natural changes in fluid flow within the mound. Daibutsu and Hobo also measured temperature anomalies during and after drilling at the ODP TAG-1 area. The Hobo temperature anomalies are inferred to have occurred when the cold fluid entrained through the drill holes at TAG-1 site reached or cooled the main fluid path to CBC. The entrained seawater through the drill holes appears to have contributed to dissolution and precipitation of anhydrite within the mound and perhaps affected the local permeability structure inside the mound. The temperature anomalies measured with Daibutsu at the base of CBC may have been induced by the change in the fluid flow pattern as a result of such permeability changes within the mound.     

Geochemical and isotopic evidence for seawater contamination of the hydrothermal system of Taal Volcano, Luzon, the Philippines

 The hydrologic structure of Taal Volcano has favored development of an extensive hydrothermal system whose prominent feature is the acidic Main Crater Lake (pH<3) lying in the center of an active vent complex, which is surrounded by a slightly alkaline caldera lake (Lake Taal). This peculiar situation makes Taal prone to frequent, and sometimes catastrophic, hydrovolcanic eruptions. Fumaroles, hot springs, and lake waters were sampled in 1991, 1992, and 1995 in order to develop a geochemical model for the hydrothermal system. The low-temperature fumarole compositions indicate strong interaction of magmatic vapors with the hydrothermal system under relatively oxidizing conditions. The thermal waters consist of highly, moderately, and weakly mineralized solutions, but none of them corresponds to either water–rock equilibrium or rock dissolution. The concentrated discharges have high Na contents (>3500 mg/kg) and low SO₄/Cl ratios (<0.3). The Br/Cl ratio of most samples suggests incorporation of seawater into the hydrothermal system. Water and dissolved sulfate isotopic compositions reveal that the Main Crater Lake and spring discharges are derived from a deep parent fluid (T≈300  °C), which is a mixture of seawater, volcanic water, and Lake Taal water. The volcanic end member is probably produced in the magmatic-hydrothermal environment during absorption of high-temperature gases into groundwater. Boiling and mixing of the parent water give rise to the range of chemical and isotopic characteristics observed in the thermal discharges. Incursion of seawater from the coastal region to the central part of the volcano is supported by the low water levels of the lakes and by the fact that Lake Taal was directly connected to the China sea until the sixteenth century. The depth to the seawater-meteoric water interface is calculated to be 80 and 160 m for the Main Crater Lake and Lake Taal, respectively. Additional data are required to infer the hydrologic structure of Taal. Geochemical surveillance of the Main Crater Lake using the SO₄/Cl, Na/K, or Mg/Cl ratio cannot be applied straightforwardly due to the presence of seawater in the hydrothermal system. Received: 12 February 1997 / Accepted: 26 January 1998     

Numerical simulation of groundwater flow patterns using flux as upper boundary

Since the 1960s, most of the studies on groundwater flow systems by analytical and numerical modelling have been based on given‐head upper boundaries. The disadvantage of the given‐head approach is that the recharge into and discharge from a basin vary with changes in hydraulic conductivity and/or basin geometry. Consequently, flow patterns simulated with given‐head boundaries but with different hydraulic conductivities and/or basin geometry may not reflect the effects of these variables. We conducted, therefore, numerical simulations of groundwater flow in theoretical drainage basins using flux as the upper boundary and realistically positioned fluid‐potential sinks while changing the infiltration intensity, hydraulic conductivities, and geometric configuration of the basin. The simulated results demonstrate that these variables are dominant factors controlling the flow pattern in a laterally closed drainage basin. The ratio of infiltration intensity to hydraulic conductivity (R_ic) has been shown to be an integrated pattern‐parameter in a basin with a given geometric configuration and possible fluid‐potential‐sink distribution. Successively, the changes in flow patterns induced by stepwise reductions in R_ic are identical, regardless of whether the reductions are due to a decrease in infiltration intensity or an increase in hydraulic conductivity. The calculated examples show five sequential flow patterns containing (i) only local, (ii) local–intermediate, (iii) local–intermediate–regional, (iv) local–regional, and (v) just regional flow systems. The R_ic was found to determine also whether a particular sink is active or not as a site of discharge. Flux upper boundary is preferable for numerical simulation when discussing the flow patterns affected by a change of infiltration, the hydraulic conductivity, or the geometry of a basin. Copyright © 2012 John Wiley & Sons, Ltd.     

Nonholonomic model of deformation of highly porous sandstone under its internal crushing

A mathematical model of elastoplastic rock deformation, under which brittle failure and shear slippages occur on the microscopic level and cataclastic flow on the macrolevel, is formulated. Since the model allows for the motion of a fluid in the pore space with its own velocity, the only applicable representation is the Eulerian representation with the Cauchy stress tensor. Naturally, the data of the rheological tests of rock samples correspond to their Lagrange representation. The problem is solved in terms of the current deformation of the material, i.e., the methodology of Mark Wilkins [Wilkins et al., 1973] applied in the adequate numerical solutions is used. The constructed model agrees with the general rules of the nonholonomic theory of plasticity, which indicates that the current parameters of the geological material depend on the path of loading. The yield surface consists of the generalized Drucker-Prager section, which corresponds to macro-slip, and the closing (at elevated pressure) cap-curve. The crushing effect corresponds to the moment when the figurative point of Mohr’s plane attains the cap-curve. The kinematics is based on nonassociative flow rules and is constructed by the summation of the known experimental data for sands and sandstones. The computer calculations of homogeneous states closely correspond to the basic types of the tests of fluid-saturated sandstone. It turned out that the rotational movements in the problems considered are so small that the corrections according to Oldroyd (which are necessary for making the formulas consistent with the theory of finite deformations) can be ignored. The changes in the coefficient of side thrust (the side pressure) appearing in the course of irreversible rock deformation are found, and the levels of crushing are evaluated.     

A theoretical study of hydrothermal convection and the origin of the ophiolitic sulphide ore deposits of Cyprus

Hydrothermal circulation of seawater has been suggested as a mass transport mechanism for the formation of sulphide ore deposits in the ophiolitic rocks of Cyprus. Since ophiolitic sequences are generally regarded as fragments of oceanic crust and upper mantle, hydrothermal circulation of a form inferred from geological observations on Cyprus may be analogous to that thought to occur in oceanic crust at spreading ridges. The hypothesis that ore deposits were formed in ascending plumes of hot, buoyant fluid is examined by considering thermal convection in a permeable medium. To match the inferred pattern of circulation, finite amplitude convection in a cylindrical geometry is studied using finite difference approximations. These results combined with available geological and geochemical data are applied to understand better the physical controls on mineralisation.A simple model for the formation of the hydrothermal ore deposits of Cyprus is discussed. The model is semi-quantitatively reasonable in terms of vertical fluid flow rate, thermal structure, permeability and basal heat flow, and predicts volumes of maximum mineralisation similar to those observed. Three factors are identified which were important in confining mineralisation to a small volume immediately beneath the sea water/rock boundary: (1) hot fluid was confined to a narrow core zone of a rising plume, (2) the upward fluid flux was greatest in this same core zone, and (3) significant temperature decrease occurred within a thin surface boundary layer.     

Pore pressure profiles in fractured and compliant rocks1

Fluid permeability in fractured rocks is sensitive to pore-pressure changes. This dependence can have large effects on the flow of fluids through rocks. We define the permeability compliance γ= 1/k(k/δp_p)_pc, which is the sensitivity of the permeability k to the pore pressure p_p at a constant confining pressure p_c, and solve the specific problems of constant pressure at the boundary of a half-space, a cylindrical cavity and a spherical cavity. The results show that when the magnitude of permeability compliance is large relative to other compliances, diffusion is masked by a piston-like pressure profile. We expect this phenomenon to occur in highly fractured and compliant rock systems where γ may be large. The pressure profile moves rapidly when fluids are pumped into the rock and very slowly when fluids are pumped out. Consequently, fluid pressure, its history and distribution around injection and production wells may be significantly different from pressures predicted by the linear diffusion equation. The propagation speed of the pressure profile, marked by the point where δp_p/δx is a maximum, decreases with time approximately as and the amplitude of the profile also dissipates with time (or distance). The effect of permeability compliance can be important for fluid injection into and withdrawal from reservoirs. For example, excessive drawdown could cause near-wellbore flow suffocation. Also, estimates of the storage capacity of reservoirs may be greatly modified when γ is large. The large near-wellbore pressure gradients caused during withdrawal by large γ can cause sanding and wellbore collapse due to excessive production rates.     

Fluid involvement in normal faulting

Evidence of fluid interaction with normal faults comes from their varied role as flow barriers or conduits in hydrocarbon basins and as hosting structures for hydrothermal mineralisation, and from fault-rock assemblages in exhumed footwalls of steep active normal faults and metamorphic core complexes. These last suggest involvement of predominantly aqueous fluids over a broad depth range, with implications for fault shear resistance and the mechanics of normal fault reactivation. A general downwards progression in fault rock assemblages (high-level breccia-gouge (often clay-rich) → cataclasites → phyllonites → mylonite → mylonitic gneiss with the onset of greenschist phyllonites occurring near the base of the seismogenic crust) is inferred for normal fault zones developed in quartzo-feldspathic continental crust. Fluid inclusion studies in hydrothermal veining from some footwall assemblages suggest a transition from hydrostatic to suprahydrostatic fluid pressures over the depth range 3–5 km, with some evidence for near-lithostatic to hydrostatic pressure cycling towards the base of the seismogenic zone in the phyllonitic assemblages. Development of fault-fracture meshes through mixed-mode brittle failure in rock-masses with strong competence layering is promoted by low effective stress in the absence of thoroughgoing cohesionless faults that are favourably oriented for reactivation. Meshes may develop around normal faults in the near-surface under hydrostatic fluid pressures to depths determined by rock tensile strength, and at greater depths in overpressured portions of normal fault zones and at stress heterogeneities, especially dilational jogs. Overpressures localised within developing normal fault zones also determine the extent to which they may reutilise existing discontinuities (for example, low-angle thrust faults). Brittle failure mode plots demonstrate that reactivation of existing low-angle faults under vertical σ₁ trajectories is only likely if fluid overpressures are localised within the fault zone and the surrounding rock retains significant tensile strength. Migrating pore fluids interact both statically and dynamically with normal faults. Static effects include consideration of the relative permeability of the faults with respect to the country rock, and juxtaposition effects which determine whether a fault is transmissive to flow or acts as an impermeable barrier. Strong directional permeability is expected in the subhorizontal σ₂ direction parallel to intersections between minor faults, extension fractures, and stylolites. Three dynamic mechanisms tied to the seismic stress cycle may contribute to fluid redistribution: (i) cycling of mean stress coupled to shear stress, sometimes leading to postfailure expulsion of fluid from vertical fractures; (ii) suction pump action at dilational fault jogs; and, (iii) fault-valve action when a normal fault transects a seal capping either uniformly overpressured crust or overpressures localised to the immediate vicinity of the fault zone at depth. The combination of σ₂ directional permeability with fluid redistribution from mean stress cycling may lead to hydraulic communication along strike, contributing to the protracted earthquake sequences that characterise normal fault systems.     

Refractive Flow and Treatment Systems: Conceptual, Analytical, and Numerical Modeling

Refractive flow and treatment (RFT) systems are designed for passive or low-maintenance in situ ground water remediation for rock or soil of low to moderate permeability. An RFT system captures and refracts contaminated ground water and conveys it to an in situ permeable treatment zone without the need for pumping. Flow to the treatment zone is through one or more high-permeability collection cells, and flow from the treatment zone back into the adjacent native media is through one or more high-permeability dispersal cells.
Conceptual, analytical, and numerical modeling demonstrates the potential for RFT systems to be successful. Analytical modeling shows that the most important factor for this success is that RFT system components be engineered to have comparatively high hydraulic conductivities. A numerical model, capable of representing site-specific conditions, is required for actual RFT system design.     

The influence of local fluid flow and the microstructure on elastic and anelastic rock properties

The frequency dependent mechanism of local fluid flow was found to be the decisive absorption and dispersion mechanism in fluid containing sandstones. In the ultrasonic frequency range local fluid flow and grain surface effects control the behaviour of highly porous and highly permeable rock if a pore fluid is present. Both mechanisms depend less on macroscopic rock parameters like porosity and permeability than essentially on microscopic parameters like crack size, crack density and grain contact properties. To demonstrate directly the important influence of the microstructure on the rock elastic and anelastic properties the microstructure of a sandstone was artificially changed by thermal cracking. The cracked rock exhibits a clearly changed behaviour at low uniaxial as well as at high hydrostatic pressure despite small changes of porosity and permeability. Fluid effects increase due to cracking. The experimental results are explained by means of a rock, model and local fluid flow. These results emphasize that it is the microstructure which controls the elastic and anelastic rock behaviour, even at high hydrostatic pressure.     

Porosity, permeability, and fluid flow in the Yellowstone geothermal system, Wyoming

Cores from two of 13 U.S. Geological Survey research holes at Yellowstone National Park (Y-5 and Y-8) were evaluated to characterize lithology, texture, alteration, and the degree and nature of fracturing and veining. Porosity and matrix permeability measurements and petrographic examination of the cores were used to evaluate the effects of lithology and hydrothermal alteration on porosity and permeability. The intervals studied in these two core holes span the conductive zone and the upper portion of the convective geothermal reservoir. Variations in porosity and matrix permeability observed in the Y-5 and Y-8 cores are primarily controlled by lithology. Y-8 intersects three distinct lithologies: volcaniclastic sandstone, perlitic rhyolitic lava, and non-welded pumiceous ash-flow tuff. The sandstone typically has high permeability and porosity, and the tuff has very high porosity and moderate permeability, while the perlitic lava has very low porosity and is essentially impermeable. Hydrothermal self-sealing appears to have generated localized permeability barriers within the reservoir. Changes in pressure and temperature in Y-8 correspond to a zone of silicification in the volcaniclastic sandstone just above the contact with the perlitic rhyolite; this silicification has significantly reduced porosity and permeability. In rocks with inherently low matrix permeability (such as densely welded ash-flow tuff), fluid flow is controlled by the fracture network. The Y-5 core hole penetrates a thick intracaldera section of the 0.6-Ma Lava Creek ash-flow tuff. In this core, the degree of welding appears to be responsible for most of the variations in porosity, matrix permeability, and the frequency of fractures and veins. Fractures are most abundant within the more densely welded sections of the tuff. However, the most prominent zones of fracturing and mineralization are associated with hydrothermal breccias within densely welded portions of the tuff. These breccia zones represent transient conduits of high fluid flow that formed by the explosive release of overpressure in the underlying geothermal reservoir and that were subsequently sealed by supersaturated geothermal fluids. In addition to this fracture sealing, hydrothermal alteration at Yellowstone appears generally to reduce matrix permeability and focus flow along fractures, where multiple pulses of fluid flow and self-sealing have occurred.     

A Poroelastic Description of Permeability Evolution

Pore pressure changes in a geothermal reservoir, as a result of injection and/or production of water, result in changes of stress acting on the reservoir rock and, consequently, changes in the mechanical and transport properties of the rock. Bulk modulus and permeability were measured at different pressures and temperatures. An outcropping equivalent of Rotliegend reservoir rock in the North German Basin (Flechtinger sandstone) was used to perform hydrostatic tests and steady state fluid flow tests. Permeability measurements were conducted while cycling confining pressure; the dependence of permeability on stress was determined at a constant downstream pressure of 1 MPa. Also, temperature was increased stepwise from 30 to 140 °C and crack porosity was calculated at different temperatures. Although changes in the volumes of cracks are not significant, the cracks control fluid flow pathways and, consequently, the permeability of the rock. A new model was derived which relates microstructure of porosity, the stress–strain curve, and permeability. Porosity change was described by the first derivative of the stress–strain curve. Permeability evolution was ascribed to crack closure and was related to the second derivative of the stress–strain curve. The porosity and permeability of Flechtinger sandstone were reduced by increasing the effective pressure and decreased after each pressure cycle.     

Fracturing and hydrothermal alteration in normal fault zones

Large normal fault zones are characterized by intense fracturing and hydrothermal alteration. Displacement is localized in a slip zone of cataclasite, breccia and phyllonite surrounding corrugated and striated fault surfaces. Slip zone rock grades into fractured, but less comminuted and hydrothermally altered rock in the transition zone, which in turn grades abruptly into the wall rock. Fracturing and fluid flow is episodic, because permeability generated during earthquakes is destroyed by hydrothermal processes during the time between earthquakes.Fracture networks are described by a fracture fabric tensor (F). The permeability tensor (k) is used to estimate fluid transport properties if the trace of F is sufficiently large. Variations in elastic moduli and seismic velocities between fault zone and wall rock are estimated as a function of fracture density (). Fracturing decreases elastic moduli in the transition zone by 50–100% relative to the country rock, and similar or even greater changes presumably occur in the slip zone.P-andS-wave velocity decrease, andV _p /V _s increases in the fault zone relative to the wall rock. Fracture permeability is highly variable, ranging between 10^–13 m² and 10^–19 m² at depths near 10 km. Changes in permeability arise from variations in effective stress and fracture sealing and healing.Hydrothermal alteration of quartzo-feldspathic rock atT>300°C creates mica, chlorite, epidote and alters the quartz content. Alteration changes elastic moduli, but the changes are much less than those caused by fracturing.P-andS-wave velocities also decrease in the hydrothermally altered fault rock relative to the country rock, and there is a slight decrease inV _p /V _s , which partially offsets the increase inV _p /V _s caused by fracturing.Fracturing and hydrothermal alteration affect fault mechanics. Low modulus rock surrounding fault surfaces increases the probability of exceeding the critical slip distance required for the onset of unstable slip during rupture initiation. Boundaries between low modulus fault rock and higher modulus wall rock also act as rupture guides and enhance rupture acceleration to dynamic velocity. Hydrothermal alteration at temperatures in excess of 300°C weakens the deeper parts of the fault zone by producingphyllitic mineral assemblages. Sealing of fracture in time periods between large earthquakes generates pods of abnormally pressured fluid which may play a fundamental role in the initiation of large earthquakes.     

Simulation of two-phase fluid flow through compactible reservoirs

Fluid-flow simulators used in the oil industry model the movement of fluids through a porous reservoir rock. These simulators either ignore coupling between the flow and concurring deformation of the solid rock frame or take it into account approximately, in the so-called loose or staggered-in-time mode. In contrast to existing simulators, the one we describe here fully couples two-phase (oil and water) flow to subsurface deformation and simultaneously accounts for all relevant physical phenomena. As such, our flow simulator inherently links time-dependent fluid pressures, saturations, permeabilities and flow velocities to stresses in the whole subsurface. These stresses relate to strains through the non-linear theory of elasticity, allowing us to model time-lapse changes in seismic velocities and anisotropy. The velocity variations manifest themselves in time shifts and reflection amplitudes that are conventionally measured from 4D seismic data. Changes in anisotropy produce time-dependent shear-wave splitting that can be used for monitoring the horizontal stresses.