Next Article in Journal
GuEstNBL: The Software for the Guided Estimation of the Natural Background Levels of the Aquifers
Next Article in Special Issue
Probabilistic Model for Real-Time Flood Operation of a Dam Based on a Deterministic Optimization Model
Previous Article in Journal
Introduction of Confidence Interval Based on Probability Limit Method Test into Non-Stationary Hydrological Frequency Analysis
Previous Article in Special Issue
Sustainability Assessment in Water Infrastructure Projects—Existing Schemes and Challenges in Application
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Geomechanical Constraints on Hydro-Seismicity: Tidal Forcing and Reservoir Operation

1
Department of Civil Engineering: Hydraulics, Energy and Environment, Universidad Politécnica de Madrid, 28040 Madrid, Spain
2
Department of Continuum Mechanics and Theory of Structures, Universidad Politécnica de Madrid, 28040 Madrid, Spain
*
Author to whom correspondence should be addressed.
Water 2020, 12(10), 2724; https://doi.org/10.3390/w12102724
Submission received: 30 July 2020 / Revised: 23 September 2020 / Accepted: 24 September 2020 / Published: 29 September 2020
(This article belongs to the Special Issue Planning and Management of Hydraulic Infrastructure)

Abstract

:
Understanding the risk associated with anthropogenic earthquakes is essential in the development and management of engineering processes and hydraulic infrastructure that may alter pore pressures and stresses at depth. The possibility of earthquakes triggered by reservoir impoundment, ocean tides, and hydrological events at the Earth surface (hydro-seismicity) has been extensively debated. The link between induced seismicity and hydrological events is currently based on statistical correlations rather than on physical mechanisms. Here, we explore the geomechanical conditions that could allow for small pore pressure changes due to reservoir management and sea level changes to propagate to depths that are compatible with earthquake triggering at critically-stressed faults (several kilometers). We consider a damaged fault zone that is embedded in a poroelastic rock matrix, and conduct fully coupled hydromechanical simulations of pressure diffusion and rock deformation. We characterize the hydraulic and geomechanical properties of fault zones that could allow for small pressure and loading changes at the ground surface (in the order of tens or hundreds of kPa) to propagate with relatively small attenuation to seismogenic depths (up to 10 km). We find that pressure diffusion to such depths is only possible for highly permeable fault zones and/or strong poroelastic coupling.

1. Introduction

The recent increase in anthropogenic earthquakes [1], especially those that are linked to industrial activities and energy technologies [2,3], has recently attracted industrial, social and scientific interest. The role of pore fluid on the occurrence, spatial distribution, and magnitude of earthquakes has become a central scientific issue. Changes in effective normal stress and shear strength to fault reactivation, due to pore pressure changes, are known to affect fault stability [4,5,6] and explain the earthquakes associated with the impoundment of water in reservoirs [7], the disposal into deep wells of waste water from oil and gas production [8,9,10,11,12], the enhanced geothermal systems [13,14,15], the large-scale CO 2 sequestration [16,17,18,19], or the hydraulic fracturing for unconventional oil and gas extraction [20,21,22]. High pore pressure may also explain slow earthquakes [23,24], aftershock triggering [25,26], or the susceptibility of critical faults to dynamic triggering [11].
Fault reactivation due to pore pressure rise is well established [4,27]. The increase in pore pressure reduces effective stress and the shear strength of the fault, while, at the same time, poroelastic effects may increase shear stress on the fault, weakening the fault and eventually leading to failure and fault reactivation [28,29,30]. However, fault reactivation is not the only condition to trigger a seismic event, since earthquakes are a manifestation of a stick-slip instability [31]. The unstable slip leads to a dynamic rupture phase that is controlled by inertia and characterized by slip velocities about meters per second and the radiation of seismic waves [32].
One of the most intriguing aspects of induced seismicity is the possibility of earthquakes that are triggered by hydrological events on the Earth surface. The increase in pore pressure or stress associated to rainfall-runoff process [33] and monsoon [34], changes in groundwater level of aquifers due to recharge or extraction operations [35,36,37,38,39], oscillations of water level in reservoirs [40,41,42,43], or even tidal cycles [44,45,46,47], have been correlated with the occurrence of earthquakes. It is essential to explore the geomechanical conditions that could allow small pore pressure changes to propagate to depths that are compatible with earthquake triggering at critically-stressed faults (several kilometers) given that the link between induced seismicity and hydrological events is currently based on statistical evidence rather than on physical mechanisms.
New water reservoirs modify the subsurface stress state via two mechanisms: the surface-mass addition due to the self-weight of the impounded water changes both the pore pressure and stress fields; and the water flow from the reservoir to the subsurface increases pore pressure [1,48]. Both theoretical studies [49,50,51] and field observations [27,52,53] have shown that fault reactivation is controlled by the initial stress state and the pore pressure, as well as the fracture network connectivity, rock permeability, fault position, and shear resistance. Earthquakes may be triggered during the first filling, after quick reservoir level changes, or even after several years of reservoir operations [7,40,54,55]. The maximum registered magnitudes range from M 4 to M 6 [42].
One of the most tragic induced earthquakes that was associated to the impoundment of a reservoir was the Koyna earthquake, in India (1967). The epicenter was 5 km away from the Koyna dam. The event was estimated to have a magnitude of M 6.3 [56]. The number of fatalities was about 200 people, and thousands more were injured. Significant structural damages to the dam were reported [57,58]. An earthquake of magnitude 7.9 on the Richter scale was triggered near the Zipingpu dam (China) in 2008, three years after the first filling. The number of fatalities may have reached 80,000 and significant structural damages to many dams in the region were reported [59]. The causal relations between the earthquake and the creation of the reservoir are still subject to debate [60].
Aquifer exploitation may also trigger earthquakes. Water injection in the subsurface increases pore pressure and decreases effective stresses, with the potential for fault reactivation [3]. Likewise, groundwater extraction may increase shear stresses on preferably-oriented faults, eventually triggering earthquakes [61]. Long-term groundwater extraction near Lorca (Spain) may have triggered a Mw 5.1 earthquake on the Richter scale [62,63,64]. Other important events that are associated to groundwater exploitation took place in Gorkha (Nepal) in 2018 with Mw 7,8 [65], the Apennines [38], the San Joaquin valley (California, USA) [66], or Oregon, USA [35].
We present a geomechanical study of the impact of pore pressure changes at the ground surface on the seismic risk at crustal depths. We consider two cases of potential hydro-seismicity: the periodic variation of sea level that is caused by tides and the impounded water level fluctuation due to reservoir management. We consider a poroelastic fault zone embedded in a poroelastic rock with smaller permeability whose hydraulic and mechanical properties may potentially allow for pressure diffusion and mechanical changes towards kilometer-scale-depths. We model the mechanical behavior of the rock as a saturated poroelastic medium that undergoes small deformations, and assume Darcy flow. The coupling between solid and fluid mass conservation is governed by the Biot equations of poroelasticity [67]. We express the formulation in terms of both the Skempton and storage coefficients. We analyze the effect of pore pressure changes on the Coulomb Failure Function (CFF) at some selected depths and assess its impact on the fault stability and, hence, on the eventual occurrence of earthquakes. We perform the computations from this constitutive model through finite element simulations.
The first set of simulations include a two-dimensional (2D) plane-strain model of a strike-slip fault immersed in a homogeneous ground domain in order to explore the hydromechanical conditions that could explain tidal triggering of earthquakes. The system is subjected to periodical tidal actions at the surface (periodic changes in pore pressure and total vertical stress). These cause variations in both the pressure field and effective stresses along the fault, so that the fault permeability may play a significant role in the pressure diffusion. The key question is under what conditions small changes in pressure at the surface can be transmitted to depths in the order of several kilometers, where hypocenters may be located. Given that the tidal periodicity is rather constant, the main parameters that lead to changes in the CFF at depth—say 5 km—are the tidal range, the fault zone permeability and storativity, as well as the Biot coefficient if full poroelastic coupling is considered. We analyze the fault response in order to obtain the evolution and amount of overpressure at such depths, and identify ranges of fault permeability, storativity, and Biot coefficient that would allow for an effective transmission of pressure to depth. We also study the effect of the fault permeability relative to that of the surrounding medium.
We consider a full three-dimensional (3D) model of a strike-slip fault underlying a water reservoir in order to explore the hydromechanical conditions that could explain triggering due to reservoir operations. The fault intersects the ground surface at reservoir site. The domain extends from the ground surface down to 20 km, to minimize boundary effects. The reservoir operations entail water level changes, albeit with either a yearly or a seasonal periodicity. We analyze the effect of such pore pressure variations on the CFF changes at given points on the fault plane. We find that the fault permeability plays a key role on the pressure diffusion and, hence, on the eventual fault instability. Unlike the tidal problem, in this case the analysis spans several years. Periodicity and the action range may both interplay with the fault permeability, which altogether affect the response lag and magnitude in terms of the CFF evolution. The ratio between the fault and the rock permeabilities shows a variety of at-depth responses. Besides, the poroelastic coupling emerges as a fault weakening factor, even when the fault permeability is rather low.

2. Mathematical Model

In the context of geotechnical engineering, the concept of specific storage coefficient arises upon consideration of the deformability of the porous medium. It expresses the changes in the volume of water stored per unit volume of soil in the subsurface that is caused by piezometric head changes, when taking the compressibility of both the water and the soil into account. In this regard, water can be accommodated into or expelled from the porous medium in several ways: (1) the solid frame may expand or compact, or even solid particles may re-accommodate, which is accounted for by the effect of K, the bulk modulus of the porous medium; (2) the solid constituent may expand or contract, accounted for by the effect of K s , the bulk modulus of the solid grains; and, (3) the fluid phase may shrink or expand according to K f , the bulk modulus of the fluid.

2.1. Flow Equation: Single-Phase Darcy Flux

Darcy’s empirical law describes the flow q of a single-phase Newtonian fluid through an undeformable porous medium with porosity ϕ as follows
q = · ( K p ρ f g ) ,
where K is the permeability tensor, p is the fluid pressure or pore pressure, ρ f is the density of the fluid, and g is the gravity vector.
The storativity of the porous medium is incorporated in the Darcy’s law through the specific storage coefficient, S. The coefficient links the increment of fluid mass stored in an elementary control volume with the deformation of the medium. The constraints of the porous medium can be either constant applied stress or constant strain, which results in two definitions of the storage coefficient: (1) the unconstrained specific storage coefficient, S σ , measured under conditions of constant applied stress; and, (2) the constrained specific storage coefficient, S ϵ , measured under conditions of constant strain [68].
The mechanical conditions in real scale geomechanical problems suggest that the constrained specific storage coefficient provides the best mathematical approximation to the real conditions [69]. Under the assumption of undeformable porous medium, the deformation of the solid grains is counteracted by the deformation of the pore space. The mass conservation equation and mechanical equilibrium equation are then uncoupled. The mass conservation equation reads:
S ϵ p t 1 μ f · ( k p ) = Q ( x , t ) ,
where μ f is the dynamic viscosity of the fluid, Q( x ,t) is a distributed source or sink function, and the constrained specific storage coefficient, S ϵ , is given by [69]:
S ϵ = ϕ χ f + ( 1 ϕ ) χ s ,
where χ f is the compressibility of the fluid and χ s is the compressibility of the solid grains.

2.2. Simplified Linear Formulation of Poroelasticity

To model the fluid flow through a deformable porous medium, we draw on the quasi-static linear poroelastic theory, which expresses the mechanical equilibrium and fluid mass conservation as functions of both the solid phase displacement field u and pore fluid pressure p as follows [67,70,71]:
· [ G u ] + ( G 1 2 ν ) · u α B p + f = 0 ,
S ϵ p t · ( k μ f p ) = α B t ( · u ) + Q ( x , t ) .
In these equations, G is the shear modulus of the solid skeleton, ν the drained Poisson ratio, f the body forces vector per unit volume of the bulk, t is time, and α B is the Biot coefficient.
The terms p in Equation (4) —momentum balance equation— and · u in Equation (5)—fluid mass conservation equation—account for the poroelastic coupling between the fluid flow and the solid skeleton deformation in space and time. The poroelastic coupling causes variations in both the pore pressure and deformation fields.
The Biot coefficient represents the ratio of increment of fluid content, ζ , to volumetric strain under drained (constant pore pressure) conditions [69]:
α B = ζ p | p = 0 .
The Biot coefficient is related to the Skempton coefficient, B, which is the ratio of volumetric strain, ϵ , to increment of fluid content under constant mean stress [69]:
B = ϵ ζ | σ = 0 .
The Skempton and Biot coefficients under undrained conditions are related in terms of the drained and undrained Poisson’s ratio, ν and ν u , respectively, by [71]:
α B = 3 ( ν u ν ) B ( 1 2 ν ) ( 1 + ν u ) ,
where ν u in a fully saturated rock with a single solid constituent is given by:
ν u = 3 ν + B α B ( 1 2 ν ) 3 B α B ( 1 2 ν ) .
The constrained specific storage, S ϵ , in a deformable medium, represents the increment of fluid content per unit control bulk volume, ζ , caused by a unit increase in pore pressure under constant control volume, i.e., ϵ = 0 . It reads [69]:
S ϵ = δ ζ δ p | ϵ = 0 .
S ϵ can be expressed in terms of the Biot coefficient, the porosity, the compressibility of the fluid, the compressibility of the solid phase χ s , and the unjacketed pore compressibility χ ϕ , as follows [72]:
S ϵ = α B χ s + ϕ ( χ f χ ϕ ) .
The compressibility of the solid phase is the compressibility of the solid grains if the solid phase consists of a single constituent, i.e., χ s = χ s . Moreover, if the porosity remains constant for equal changes in pore and confining stress, the unjacketed pore compressibility equals the compressibility of the solid phase, i.e., χ ϕ = χ s . Assuming both, S ϵ can be expressed in terms of the Biot coefficient, the porosity, and the compressibilities of the solid and fluid components, as follows [72]:
S ϵ = α B χ s + ϕ ( χ f χ s ) .
The above expression is expressed in terms of the bulk modulus K, as follows [73,74]:
S ϵ = ϕ χ f + 1 K ( α B ϕ ) ( 1 α B ) .
For practical purposes, it is easier to calculate S ϵ by measuring Skempton’s coefficient B and avoid the difficulties of measuring χ n . Hence, it is convenient to express S ϵ in terms of B [71]:
S ϵ = α B 2 ( 1 2 ν ) ( 1 2 ν u ) 2 G ( ν u ν ) .

3. Model Setup

3.1. Model Implementation

We have studied two hydromechanical problems, where variations in the Coulomb Failure Function (CFF) at depth, associated with pore pressure changes on the ground surface, can eventually trigger an earthquake. The CFF characterizes stress changes with regard to stabilizing/destabilizing a fault, based on the Mohr–Coulomb analogy for frictional contact. Hence, it is defined as Δ CFF = μ Δ σ Δ τ , where μ denotes the friction coefficient, Δ σ are the changes in effective compression on the fault, and Δ τ are changes in shear stress acting on the fault.
In the first problem, we study how tides may affect the CFF at depths up to 5 km on a strike-slip fault, depending on the hydromechanical properties of the fault zone and surrounding rock. In the second problem, we analyze the effect of reservoir operations on the CFF at depths up to 20 km. We illustrate both problems and their respective computational domains in Figure 1. We solve both problems using the Finite Element Method. We adopt Taylor–Hood elements with linear discretization for pressure and quadratic one for displacements. We assume the fault permeability is independent of the effective stress evolution, since significant changes in fault permeability require a decrease in effective stress of several MPa [75]. In our problem, the pressure increases are up to several kPa.
The oscillation of the sea level induces pressure changes on the surface, as shown in Figure 1a. The hydraulic properties of faults may make these to work as preferential channels for pressure diffusion towards depth. Our case study is a strike–slip fault that extends from the ground surface down to 5 km deep. We simulate this problem with a 2D vertical domain perpendicular to the fault plane, Figure 1b. We assume plane-strain elasticity for the full poroelastic problem, and mesh the domain through unstructured grids with highly refined elements that are close to the fault and the top boundary, where the tidal load is applied. We use triangular elements with maximum side of 6 m on the faults and the upper boundaries.
The oscillation of the water level in a reservoir due to the regular operation tasks alters pore pressure. Because rivers typically flow in patterns that follow weaknesses in the bedrock, such as faults, reservoirs are usually built over faults, as shown in Figure 1c. We analyze the effect of such water level oscillations on the deep-CFF due to the pressure diffusion along a strike–slip fault that extends from the surface of the ground down to 20 km deep. We simulate this problem with a 3D domain that includes the fault and a 20 km squared portion of the Earth surface. We mesh the domain while using tetrahedral elements with maximum sizes of 100 m and 400 m in the reservoir and the fault, respectively, Figure 1d.

3.2. Tidal Oscillation Model Set-Up

We study the effect of the tidal oscillation on a strike–slip fault with a 2D model, as shown in Figure 1b. We simulate the fault region as a rectangle centered in the squared domain. We list the model parameters in Table 1 and Table 2. As we assume linear flow and linear poroelasticity, the initial pore pressure takes an arbitrary reference value p = p ref . Our calculations address tidal perturbations as a superposition with respect to hydrostatic/lithostatic conditions. We impose no flow on the vertical and bottom boundaries and consider an applied pressure on the upper surface, as given by p r e f + ρ f g h ( t ) , where ρ f is the fluid density and h ( t ) is the tidal sea-level oscillation. We compute h ( t ) with a sinusoidal wave:
h ( t ) = A · sin 2 π T · t ,
where A is the amplitude of the tide and T the period.
The mechanical boundary conditions are null normal displacements, except on the horizontal top boundary where we impose the normal load generated by the tidal oscillation. We assume plane strain conditions and perfect continuity of pore pressures and solid displacements at the contact between the rock matrix and the fault zone.

3.3. Reservoir Oscillation Model Set-Up

We simulate the effect of a reservoir operation on the stability of a strike–slip fault with a 3D model. The fault zone is modeled as an orthohedral subdomain, whose width is W and is centered in a 20 km cubic domain (Figure 1c). Table 3 lists the parameters of the model.
We analyze the effect of reservoir level changes, which cause an additional stress and pressure distribution. Hence, initial pore pressure takes an arbitrary reference value p = p ref . We impose no-flow condition on all boundaries, except at the bottom of the reservoir. The reservoir is modeled as a 1 km side square that is centered on the top boundary, and the water level oscillates following a sinusoidal wave given by Equation (15). The mechanical boundary conditions are null normal displacements on the vertical and bottom horizontal boundaries. On the top horizontal boundary at the reservoir bottom, we impose a normal load equal to the pressure head and, on the remaining upper boundary, we impose a free displacement condition. We assume continuity on the contact between the host rock and the fault.

4. Results and Discussion

The fluid flow and poroelastic effects govern pore pressure transmission in a saturated porous medium. The former is described by Darcy’s law, whereas the Biot equations of poroelasticity encompass both effects. Changes of water level as a consequence of either ocean tides or reservoir oscillations on the ground surface may induce pressure changes in depth as a result of both, the transport of fluid governed by Darcy’s law, and the pore pressure diffusion due to the water load on the surface as a consequence of poroelastic effects. Here, we elucidate how water level oscillations on surface induce pressure changes at great depth, as well as the role of physical properties of faults, Darcy’s law, and poroelastic effects in the pressure transmission and stress changes.

4.1. Tidal Oscillation

We analyze the effect of tidal oscillations on pore pressure changes along a strike-slip fault through two models. The first approach assumes the porous medium to be undeformable and it is modeled through Equation (2). We denote this approach as the Darcy simplification. The accumulation of fluid in the rock varies with pore pressure and it is proportional to the compressibilities of the fluid and solid material of the porous medium, χ f and χ s . Because the porous medium is assumed to be undeformable, i.e., the volumetric deformation is null, so that pressures affect the pore volume available for flow, but otherwise fluid flow is uncoupled with porous matrix deformation. The second approach accounts for the deformation in the porous matrix. Hence, fluid flow is fully coupled with porous matrix deformation by means of the Biot equations of poroelasticity, Equations (4) and (5). The comparison between both model outputs allow for us to assess which part of the pore pressure change at depth is caused by the rock matrix deformation due to the poroelastic effect, and which one is driven by fluid transport and the local storage capacity of the fault zone.

4.1.1. The Effect of the Fault Permeability

The rock and fault permeabilities and their ratio play an important role in the pore pressure transmission. We conduct a set of simulations where, by taking a small constant value of the host rock matrix permeability, we vary the fault permeability and analyze changes in the pore pressure on the fault. We adopt an undeformable porous medium (Darcy model). We depict several vertical contour plots of the pressure field when the change in pressure reaches the deepest regions of the fault presented in Figure 2. The fault permeability is always equal or higher than the rock permeability, so that the fault zone behaves as a preferential path for pressure diffusion and fluid flow.
The equipotential lines are parallel to the ground surface when both permeabilities are equal, as in Figure 2a. However, the increase in fault permeability makes the lines deviate from the ground surface and the rise in pore pressure extends deeper, as in Figure 2b. The higher the fault zone permeability, the deeper propagation of the pore pressure, Figure 2c, reaching even 5000 m depth for reasonably large permeabilities (Figure 2d).
The overpressure on the fault attenuates with depth, as in Figure 3. The compressibility of the fluid and grains make the amplitude of the overpressure wave decrease and delay with depth and elapsed time. The attenuation of the amplitude decreases with permeability: the larger the permeability the lower the attenuation (Figure 3a,b). Fault permeability may also play a significant role in fault stability: the larger the permeability, the larger the mean overpressure along the fault and, consequently, the higher destabilization risk. Moreover, the location of the maximum instantaneous overpressure is deeper as the permeability is lower: it is always on the ground surface for the simulation with the highest fault permeability (Figure 3b), whereas it descends as the fault permeability decreases (Figure 3a).
We quantify the propagation of the overpressure across depth for a given fault permeability through the position of the point, where the overpressure amplitude attenuates 90% (Figure 4a). Low permeability values impede the propagation of the overpressure, but, as permeability increases, the wave penetrates further. Even for high fault permeability, the wave reaches the bottom boundary of our model, which is located at 5000 m depth.

4.1.2. The Poroelastic Coupling

The Biot equations of poroelasticity account for the compressibilities of fluid, solid material of the porous medium, and porous matrix, whereas the uncoupled Darcy formulation only considers the compressibilities of fluid and solid material. Low compressibility values of the solid material enhance the pore pressure transmission (Figure 4b). The maximum attenuation of the pressure at the bottom boundary of our numerical model is approximately 60% for the lowest values of fault permeability, whereas the attenuation becomes total under the uncoupled Darcy formulation. Both of the approaches provide similar results for high values of the fault permeability.
The pressure diffusion coefficient, D = k F / ( μ f S ϵ ) , plays an important role in the pore pressure transmission across depth. We have run five sets of simulations with constant values of fault permeability and varying D (Figure 5). The depth at which the pressure attenuates 90% increases with D for a given k F . The main factor controlling the overpressure advance depends on fault permeability: for low k F , the fault permeability mostly controls the overpressure advance, whereas the diffusion coefficient controls the penetration for those cases with high k F .
The coupling between fluid and mechanical equilibrium allows for us to study the impact of the Young modulus and the Biot coefficient on the attenuation of the overpressure at a depth of 5 km (Figure 6). The Young modulus slightly alters the overpressure wave. The lower the modulus, the higher the overpressure. Nevertheless, the Biot coefficient has great influence: the higher this coefficient the higher the overpressure, except for highly permeable faults, where overpressure mainly diffuses by fluid flow and with poroelastic effects hardly playing any role.
The evolution of the CFF, defined as Δ CFF = μ Δ σ n Δ τ , provides an insight into the likelihood of fault slip that eventually may trigger an earthquake (Figure 7). The tide oscillation induces an increase in pore pressure, which reduces effective stress on the fault, while, at the same time, it alters the shear stress. The change in pore pressure is the main precursor of the CFF evolution, whereas shear stress remains almost constant in time. Consequently, fault permeability is an important property that controls the evolution of CFF (Figure 7a,b).

4.1.3. The Role of the Host Rock Permeability

The ratio of the host rock to fault permeabilities may influence the pressure changes during the tidal cycle. We run a set of simulations by keeping constant the fault permeability and varying the host rock permeability. We simulate the host rock and fault as poroelastic solids, according to Equations (4) and (5). The output of the simulations provides valuable information on the effect of the host rock permeability and poroelastic coupling on the pore pressure diffusion. The results are presented in Figure 8 and Figure 9.
The permeability of the host rock, k, varies from a value that is equal to the fault permeability, k F , up to an almost impermeable value 10 5 · k F , Figure 8. The permeability of the fault is either equal or higher than the host rock, so that the former becomes a preferential way for water flow. The pressure diffusion along the fault generates a horizontal pressure gradient that transfers water from the fault to the host rock.
The extreme values of the host rock permeability produce two opposite behaviors in its pressure field, yet the fault pressure field remains the same. In the first scenario, both the fault and the host rock have the same permeability, Figure 8b. The movement of the fluid front is horizontal throughout the whole domain, and the horizontal pressure gradient is null. There is no fluid transfer between the fault and the host rock. The second scenario involves an almost impermeable host rock. Water does not flow through the rock, but along the fault, as in Figure 8c. The horizontal pressure gradients between the fault and the host rock are high, but the extremely low permeability of the rock does not allow water to flow from the fault to the rock.
The decrease in the permeability of the host rock from k = k F to k = 10 5 · k F implies that less amount of water flows throughout the rock. However, because the fault works as a preferential way for water to flow, the fluid shifts from the fault to the host rock. The equipotential lines are no longer horizontal, as in Figure 8a. The effect of the host rock permeability on the fluid flow along the fault may be classified into two categories. For k-values two or three orders of magnitude lower than k F , the host rock does not allow for the transport of water from the seabed to the inner area of the rock. However, the rock permeability is high enough to let water depart from the fault, as in Figure 8e,f, and the pressure cannot increase in the deepest parts of the fault. However, host rocks with k-values four orders of magnitude lower than k F are impermeable enough to reduce water leaks and allow for pressure to increase in the deepest areas of the fault, as in Figure 8c,d.
The poroelastic coupling may play an important role in the fault pressure build-up. We study this effect by running two sets of simulations with Biot’s coefficient being 0 and 0.8, respectively. The fault permeability is kept constant in all of the simulations, whereas we vary the ratio k / k F by considering several host rock permeabilities. The results are depicted in Figure 9, where we plot the fault pressure build-up at 5000 m depth.
When disregarding the Biot coefficient effect in the simulations, pressure changes are caused exclusively by water transport, whereas, for the simulations with α B = 0.8 , the poroelastic coupling induces pressure changes due to stress alterations in the solid. Both simulation sets, α B = 0 and α B = 0.8 , present a similar V-shape evolution of the pressure build-up at 5000 m depth with the ratio k / k F . The pressure build-up is maximum for k / k F = 1 , it decreases with k / k F up to a value of the ratio from which the pressure starts to increase with the decrease in k / k F . The shape of the curve lies on the phenomenon that was described in the previous paragraphs: the decrease in k promotes the migration of fluid from the fault to the host rock, but from a low enough k-value the host rock works as an impermeable material and water penetrates deeper inside the fault.
The poroelastic coupling describes the interaction between pore pressure and solid deformation. The tidal cycle induces changes in the water load that was applied on the seabed, which may alter pore pressure at large depths. The minimum pressure build-up that is shown in Figure 9 is almost zero and given for the set with α B = 0 , i.e., without poroelastic coupling. However, it is 25% of the pressure on the seabed for the set with α B = 0.8 . Moreover, the poroelastic coupling makes the pressure build-up always higher than simulations with α B = 0 .

4.2. Reservoir Oscillation

We analyze the effect of the reservoir oscillations on the pore pressure build-up through three-dimensional (3D) numerical simulations. Our model couples fluid flow and mechanical equilibrium by adopting the Biot equations of poroelasticity. The dam is 120 m high and the reservoir oscillation is also 100 m. We consider two operation regimes. The first one corresponds to a reservoir with annual regulation, i.e., the water level is zero at the beginning of the year, reaches its maximum value at mid-year, and it comes down to the minimum level at the end of the year. This operation regime is typical of irrigation, water supply, or power reservoirs. The second regime is a daily regulation, which is typical of pumped-storage hydropower reservoirs.
The seismic behavior of a fault depends on both the rate of fluid flow and rate of pore pressure increase, p / t : the higher the fluid flow or p / t , the more likely is the fault to trigger an earthquake [76]. We analyze the evolution of these variables in terms of p and p / t in our case study in the following paragraphs.
Figure 10 illustrates the maximum pressure build-up on a fault crossing an annual regulation reservoir. The maximum pressure is located under the reservoir and attenuates with depth. Moreover, the fault permeability exerts a considerable influence on the extension of the area, where pore pressure changes: the higher the permeability, the larger the affected area.
We analyze the impact of the reservoir operation on the fault pore pressure build-up in Figure 11. The maximum and minimum values of the pressure are quite similar for both reservoir regulations up to approximately 5 km depth. The pressure attenuates about 60% at 1 km depth, 20% at 3 km depth, and 15% at 5 km depth, as in Figure 11a,b. The attenuation is slightly higher for the reservoir with annual operation than the impoundment with daily operation. Nevertheless, differences arise at deepest regions: the attenuation is almost constant in time and equal to 5% at 20 km depth for the daily regulation, whereas it oscillates up to 20% for the annual regulation. The operation of the reservoir has important influence on the pressure diffusion: the longer the oscillation period, the larger the extent of the area of rising pressure, as in Figure 11c,d.
Figure 12 illustrates the evolution of p / t for the two reservoir operations. The faster the reservoir oscillation, the larger the rate p / t , Figure 12a,b. The daily operation is found to produce two orders of magnitude higher rate p/t than the annual operation. The rate p / t attenuates with depth. The attenuation is more pronounced for the daily than for the annual regulation: p / t is almost zero at 20 km depth for the reservoir with daily regulation. Figure 12c,d show the area with appreciable changes in p / t on the fault. Although the daily regulation leads to changes in p / t higher than the annual operation, the affected area with significant variations in p / t is much smaller.
The permeability of the fault as well as the poroelastic coupling play a key role in both the pore pressure attenuation and ratio p / t , as in Figure 13. Poroelastic effects on high permeable faults are almost negligible, the pressure attenuation or the ratio p / t are almost independent of the value of the Biot coefficient, Figure 13a,b. Pressure changes are mostly due to the fluid flow along the fault, and poroelastic effects have little effect on them. Nevertheless, poroelastic effects play a key role on almost impermeable faults, as in Figure 13c,d. The higher the coupling between fluid flow and solid deformation, the higher the pressure build-up and ratio p / t . Moreover, the values of these variables are also higher in impermeable faults than for permeable faults.

5. Conclusions

We have analyzed the geomechanical processes that may explain how water level changes on the ground surface can affect the stability of critically-stressed faults, and that eventually may trigger an earthquake. We have employed two-dimensional and three-dimensional numerical simulations. We have studied two problems: in the first one we analyze the effect of the tidal oscillation on the stability of a strike-slip fault, and in the second one we elucidate the consequences of reservoir oscillations. We adopt a tidal oscillation of 6 m and a reservoir oscillation of 100 m. We simulate two reservoir operations: annual regulation, typical of irrigation, water supply or power reservoirs, and daily regulation, which are common in pumped-storage hydropower reservoirs. Our numerical models include a domain of 5000 m deep for the tidal problem and 20 km deep for the reservoir oscillation simulations.
We quantify fault stability through the change in the Coulomb Failure Function (CFF) with respect to the initial stress state, Δ CFF , defined as Δ CFF = μ Δ σ n Δ τ , where σ n is the normal effective stress on the fault and τ is the shear stress on the fault. The seismic behavior of a fault also depends on both the rate of fluid flow and rate of pore pressure increase, p / t : the higher the fluid flow or p / t , the more likely is the fault to trigger an earthquake. We analyze the evolution of these quantities in terms of pressure change, Δ p , and p / t .
The pore overpressure evolves due to two mechanisms: the fluid flow and the poroelastic effects. Our simulations have shown that the permeability of the fault plays an important role in the overpressure transmission due to the mass flow: the higher the permeability, the deeper the overpressure changes and the lower attenuation of the pressure with depth. Moreover, poroelastic effects enhance the overpressure transmission, especially in faults with low permeability, where the pore overpressure transmission by fluid flow is almost impeded. There is an overpressure attenuation driven by fault permeability, and it is not accurate for neglecting this attenuation.
The ratio of the host rock permeability to the fault permeability considerably affects the pore pressure build-up. Host rocks with permeability that is up to three orders of magnitude lower than fault permeability produce water leakages from the fault to the host rock. Pressure does not increase in the deepest parts of the fault. Nevertheless, host rocks with permeability four orders of magnitude lower than that of the fault work as an impermeable layer and leakages are nearly null. Pressure then changes in the deepest areas of the fault.
The mechanical properties of the fault and the host rock, as well as the Biot coefficient, have a significant influence on the overpressure transmission. We assess the effects of these properties through the Biot equations of poroelasticity. The stiffness of the host rock slightly alters the overpressure transmission: the higher the stiffness, the lower the overpressure. The role of the Biot coefficient is more pronounced: the higher this coefficient, the higher the overpressure.
We have found reservoir operations to be a key parameter on the evolution of the rate p / t . The faster the reservoir oscillation, the larger the rate p / t . The daily operation is found to increase by two orders of magnitude the rate p / t as compared with the annual operation. The rate p / t attenuates with depth. The attenuation is more pronounced for the daily regulation than the annual.
The likelihood of triggering an earthquake due to water level changes on the Earth surface is conditioned by the initial stress state of the fault. Critically stressed faults can be destabilized by small increases in pore pressure, due to, for instance, changes in tidal oscillation. Our simulations have shown that pore pressure is the main mechanism responsible for the fault destabilization, whereas changes in shear stress are negligible. Our methodology arises as a useful tool for engineers, practitioners, and stakeholders to assess fault reactivation for a new project or existing facilities.

Author Contributions

Conceptualization, L.C.-F.; Data curation, P.P. and L.C.-F.; Formal analysis, D.S. and J.C.M.; Investigation, P.P., D.S. and L.C.-F.; Methodology, D.S.; Writing—original draft, P.P., D.S. and J.C.M.; Writing—review & editing, J.C.M. and L.C.-F. All authors have read and agreed to the published version of the manuscript.

Funding

Project supported by a 2019 Leonardo Grant for Researchers and Cultural Creators, BBVA Foundation. P.P. gratefully acknowledges funding from the Spanish Ministry of Science, Innovation and Universities through Beca FPU (Formación de Profesorado Universitario).

Conflicts of Interest

The authors declare no conflict of interest. The BBVA Foundation accepts no responsibility for the opinions, statements and contents included in the project and/or the results thereof, which are entirely the responsibility of the authors.

References

  1. Foulger, G.R.; Wilson, M.P.; Gluyas, J.G.; Julian, B.R.; Davies, R.J. Global review of human-induced earthquakes. Earth Sci. Rev. 2018, 178, 438–514. [Google Scholar] [CrossRef] [Green Version]
  2. National Research Council. Induced Seismicity Potential in Energy Technologies; The National Academies Press: Washington, DC, USA, 2013. [Google Scholar]
  3. Ellsworth, W.L. Injection-Induced Earthquakes. Science 2013, 341, 1225942. [Google Scholar] [CrossRef] [PubMed]
  4. Hubbert, M.K.; Rubey, W.W. Role of fluid pressure in mechanics of overthrust faulting I. Mechanics of fluid-filled porous solids and its application to overthrust faulting. Geol. Soc. Am. Bull. 1959, 70, 115–166. [Google Scholar]
  5. Raleigh, C.; Healy, J.; Bredehoeft, J. An experiment in earthquake control at Rangely, Colorado. Science 1976, 191, 1230–1237. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Sibson, R.H. A note on fault reactivation. J. Struct. Geol. 1985, 7, 751–754. [Google Scholar] [CrossRef]
  7. Simpson, D.; Leith, W.; Scholz, C. Two types of reservoir-induced seismicity. Bull. Seismol. Soc. Am. 1988, 78, 2025–2040. [Google Scholar]
  8. Horton, S. Disposal of hydrofracking waste fluid by injection into subsurface aquifers triggers earthquake swarm in central Arkansas with potential for damaging earthquake. Seismol. Res. Lett. 2012, 83, 250–260. [Google Scholar] [CrossRef] [Green Version]
  9. Keranen, K.M.; Savage, H.M.; Abers, G.A.; Cochran, E.S. Potentially induced earthquakes in Oklahoma, USA: Links between wastewater injection and the 2011 Mw 5.7 earthquake sequence. Geology 2013, 41, 699–702. [Google Scholar] [CrossRef]
  10. Keranen, K.M.; Weingarten, M.; Abers, G.A.; Bekins, B.A.; Ge, S. Sharp increase in central Oklahoma seismicity since 2008 induced by massive wastewater injection. Science 2014, 345, 448–451. [Google Scholar] [CrossRef]
  11. Weingarten, M.; Ge, S.; Godt, J.W.; Bekins, B.A.; Rubinstein, J.L. High-rate injection is associated with the increase in US mid-continent seismicity. Science 2015, 348, 1336–1340. [Google Scholar] [CrossRef] [Green Version]
  12. Shirzaei, M.; Ellsworth, W.L.; Tiampo, K.F.; González, P.J.; Manga, M. Surface uplift and time-dependent seismic hazard due to fluid injection in eastern Texas. Science 2016, 353, 1416–1419. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Majer, E.L.; Baria, R.; Stark, M.; Oates, S.; Bommer, J.; Smith, B.; Asanuma, H. Induced seismicity associated with enhanced geothermal systems. Geothermics 2007, 36, 185–222. [Google Scholar] [CrossRef]
  14. Brodsky, E.; Lajoie, L. Anthropogenic Seismicity Rates and Operational Parameters at the Salton Sea Geothermal Field. Science 2013, 341, 543–546. [Google Scholar] [CrossRef] [PubMed]
  15. Andrés, S.; Santillán, D.; Mosquera, J.C.; Cueto-Felgueroso, L. Thermo-Poroelastic Analysis of Induced Seismicity at the Basel Enhanced Geothermal System. Sustainability 2019, 11, 6904. [Google Scholar]
  16. Zoback, M.D.; Gorelick, S.M. Earthquake triggering and large-scale geologic storage of carbon dioxide. Proc. Natl. Acad. Sci. USA 2012, 109, 10164–10168. [Google Scholar] [CrossRef] [Green Version]
  17. Juanes, R.; Hager, B.H.; Herzog, H.J. No geologic evidence that seismicity causes fault leakage that would render large-scale carbon capture and storage unsuccessful. Proc. Natl. Acad. Sci. USA 2012, 109, E3623. [Google Scholar] [CrossRef] [Green Version]
  18. Vilarrasa, V.; Carrera, J. Geologic carbon storage is unlikely to trigger large earthquakes and reactivate faults through which CO2 could leak. Proc. Natl. Acad. Sci. USA 2015, 112, 5938–5943. [Google Scholar] [CrossRef] [Green Version]
  19. White, J.A.; Foxall, W. Assessing induced seismicity risk at CO2 storage projects: Recent progress and remaining challenges. Int. J. Greenhouse Gas Control 2016, 49, 413–424. [Google Scholar] [CrossRef] [Green Version]
  20. Clarke, H.; Eisner, L.; Styles, P.; Turner, P. Felt seismicity associated with shale gas hydraulic fracturing: The first documented example in Europe. Geophys. Res. Lett. 2014, 41, 8308–8314. [Google Scholar] [CrossRef]
  21. Bao, X.; Eaton, D. Fault activation by hydraulic fracturing in western Canada. Science 2016, 354, 1406–1409. [Google Scholar] [CrossRef]
  22. Deng, K.; Liu, Y.; Harrington, R.M. Poroelastic stress triggering of the December 2013 Crooked Lake, Alberta, induced seismicity sequence. Geophys. Res. Lett. 2016, 43, 8482–8491. [Google Scholar] [CrossRef]
  23. Kodaira, S.; Iidaka, T.; Kato, A.; Park, J.O.; Iwasaki, T.; Kaneda, Y. High pore fluid pressure may cause silent slip in the Nankai Trough. Science 2004, 304, 1295–1298. [Google Scholar] [CrossRef] [PubMed]
  24. Becken, M.; Ritter, O.; Bedrosian, P.A.; Weckmann, U. Correlation between deep fluids, tremor and creep along the central San Andreas fault. Nature 2011, 480, 87. [Google Scholar] [CrossRef] [PubMed]
  25. Nur, A.; Booker, J. Aftershocks Caused by Pore Fluid flow? Science 1972, 175, 885–887. [Google Scholar] [CrossRef]
  26. Miller, S.A.; Collettini, C.; Chiaraluce, L.; Cocco, M.; Barchi, M.; Kaus, B.J. Aftershocks driven by a high-pressure CO2 source at depth. Nature 2004, 427, 724. [Google Scholar] [CrossRef] [PubMed]
  27. Healy, J.; Rubey, W.; Griggs, D.; Raleigh, C. The Denver earthquakes. Science 1968, 161, 1301–1310. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Cueto-Felgueroso, L.; Santillán, D.; Mosquera, J.C. Stick-slip dynamics of flow-induced seismicity on rate and state faults. Geophys. Res. Lett. 2017, 44, 4098–4106. [Google Scholar] [CrossRef]
  29. Cueto-Felgueroso, L.; Vila, C.; Santillan, D.; Mosquera, J.C. Numerical modeling of injection–induced earthquakes using laboratory–derived friction laws. Water Resour. Res. 2018, 54, 1–27. [Google Scholar] [CrossRef]
  30. Andrés, S.; Santillán, D.; Mosquera, J.C.; Cueto-Felgueroso, L. Delayed weakening and reactivation of rate-and-state faults driven by pressure changes due to fluid injection. J. Geophys. Res. Solid Earth 2019, 124, 11917–11937. [Google Scholar] [CrossRef]
  31. Scholz, C.H. Earthquakes and friction laws. Nature 1998, 391, 37–42. [Google Scholar] [CrossRef]
  32. Pampillón, P.; Santillán, D.; Mosquera, J.C.; Cueto-Felgueroso, L. Dynamic and Quasi–Dynamic Modeling of Injection–Induced Earthquakes in Poroelastic Media. J. Geophys. Res. Solid Earth 2018, 123, 5730–5759. [Google Scholar] [CrossRef]
  33. Hainzl, S.; Kraft, T.; Wassermann, J.; Igel, H.; Schmedes, E. Evidence for rainfall-triggered earthquake activity. Geophys. Res. Lett. 2006, 33. [Google Scholar] [CrossRef] [Green Version]
  34. Hainzl, S.; Aggarwal, S.; Khan, P.; Rastogi, B. Monsoon-induced earthquake activity in Talala, Gujarat, India. Geophys. J. Int. 2015, 200, 627–637. [Google Scholar] [CrossRef] [Green Version]
  35. Saar, M.O.; Manga, M. Seismicity induced by seasonal groundwater recharge at Mt. Hood, Oregon. Earth Planet. Sci. Lett. 2003, 214, 605–618. [Google Scholar] [CrossRef]
  36. D’Agostino, N.; Silverii, F.; Amoroso, O.; Convertito, V.; Fiorillo, F.; Ventafridda, G.; Zollo, A. Crustal Deformation and Seismicity Modulated by Groundwater Recharge of Karst Aquifers. Geophys. Res. Lett. 2018, 45, 12253–12262. [Google Scholar] [CrossRef]
  37. Bollinger, L.; Nicolas, M.; Marin, S. Hydrological triggering of the seismicity around a salt diapir in Castellane, France. Earth Planet. Sci. Lett. 2010, 290, 20–29. [Google Scholar] [CrossRef]
  38. Bella, F.; Biagi, P.; Caputo, M.; Cozzi, E.; Monica, G.; Ermini, A.; Plastino, W.; Sgrigna, V. Aquifer-induced Seismicity in the Central Apennines (Italy). Pure Appl. Geophys. 1998, 153, 179–194. [Google Scholar] [CrossRef]
  39. Heinicke, J.; Woith, H.; Alexandrakis, C.; Buske, S.; Telesca, L. Can hydroseismicity explain recurring earthquake swarms in NW-Bohemia? Geophys. J. Int. 2018, 212, 211–228. [Google Scholar] [CrossRef] [Green Version]
  40. Talwani, P. On the nature of reservoir-induced seismicity. Pure Appl. Geophys. 1997, 150, 473–492. [Google Scholar] [CrossRef]
  41. Klose, C. Mechanical and statistical evidence of the causality of human-made mass shifts on the Earth’s upper crust and the occurrence of earthquakes. J. Seismol. 2013, 17, 109–135. [Google Scholar] [CrossRef]
  42. Gupta, H. A review of recent studies of triggered earthquakes by artificial water reservoirs with special emphasis on earthquakes in Koyna, India. Earth Sci. Rev. 2002, 58, 279–310. [Google Scholar] [CrossRef]
  43. Nascimento, A.; Lunn, R.; Cowie, P. Numerical modelling of pore-pressure diffusion in a reservoir-induced seismicity site in northeast Brazil. Geophys. J. Int. 2005, 160, 249–262. [Google Scholar] [CrossRef] [Green Version]
  44. Delorey, A.A.; van der Elst, N.J.; Johnson, P.A. Tidal triggering of earthquakes suggests poroelastic behavior on the San Andreas Fault. Earth Planet. Sci. Lett. 2017, 460, 164–170. [Google Scholar] [CrossRef] [Green Version]
  45. Thomas, A.; Burgmann, R.; Shelly, D.; Beeler, N.; Rudolph, M. Tidal triggering of low frequency earthquakes near Parkfield, California: Implications for fault mechanics within the brittle-ductile transition. J. Geophys. Res. Solid Earth 2012, 117, B05301. [Google Scholar] [CrossRef]
  46. Thomas, A.; Nadeau, R.; Burgmann, R. Tremor-tide correlations and near-lithostatic pore pressure on the deep San Andreas fault. Nature 2009, 462, 1048–1051. [Google Scholar] [CrossRef]
  47. Hainzl, S.; Ben-Zion, Y.; Cattania, C.; Wassermann, J. Testing atmospheric and tidal earthquake triggering at Mt. Hochstaufen, Germany. J. Geophys. Res. Solid Earth 2013, 118, 5442–5452. [Google Scholar] [CrossRef] [Green Version]
  48. El-Hariri, M.; Abercrombie, R.; Rowe, C.; do Nascimento, A. The role of fluids in triggering earthquakes: observations from reservoir induced seismicity in Brazil. Geophys. J. Int. 2010, 181, 1566–1574. [Google Scholar] [CrossRef] [Green Version]
  49. Bell, M.; Nur, A. Strength changes due to reservoir-induced pore pressure and stresses and application to Lake Oroville. J. Geophys. Res. Solid Earth 1978, 83, 4469–4483. [Google Scholar] [CrossRef]
  50. Simpson, D. Triggered Earthquakes. Annu. Rev. Earth Planet. Sci. 1986, 14, 21–42. [Google Scholar] [CrossRef]
  51. Roeloffs, E. Fault stability changes induced beneath a reservoir with cyclic variations in water level. J. Geophys. Res. Solid Earth 1988, 93, 2107–2124. [Google Scholar] [CrossRef]
  52. Zoback, M.; Hickman, S. In situ study of the physical mechanisms controlling induced seismicity at Monticello Reservoir, South Carolina. J. Geophys. Res. Solid Earth 1982, 87, 6959–6974. [Google Scholar] [CrossRef] [Green Version]
  53. Fletcher, J. A comparison between the tectonic stress measured in situ and stress parameters from induced seismicity at Monticello Reservoir, South Carolina. J. Geophys. Res. Solid Earth 1982, 87, 6931–6944. [Google Scholar] [CrossRef]
  54. Matcharashvili, T.; Chelidze, T.; Peinke, J. Increase of order in seismic processes around large reservoir induced by water level periodic variation. Nonlinear Dyn. 2008, 51, 399–407. [Google Scholar] [CrossRef]
  55. Peinke, J.; Matcharashvili, T.; Chelidze, T.; Gogiashvili, J.; Nawroth, A.; Lursmanashvili, O.; Javakhishvili, Z. Influence of periodic variations in water level on regional seismic activity around a large reservoir: Field data and laboratory model. Phys. Earth Planet. Inter. 2006, 156, 130–142. [Google Scholar] [CrossRef]
  56. Gupta, H.K. 50 years of 10 December 1967 M 6.3 Koyna earthquake. J. Geol. Soc. India 2017, 90, 641–644. [Google Scholar] [CrossRef] [Green Version]
  57. Narain, H.; Gupta, H. Koyna earthquake. Nature 1968, 217, 1138–1139. [Google Scholar] [CrossRef]
  58. Santillan, D.; Mosquera, J.; Cueto-Felgueroso, L. Phase-field model for brittle fracture. Validation with experimental results and extension to dam engineering problems. Eng. Fract. Mech. 2017, 178, 109–125. [Google Scholar] [CrossRef]
  59. Ge, S.; Liu, M.; Lu, N.; Godt, J.; Luo, G. Did the Zipingpu Reservoir trigger the 2008 Wenchuan earthquake? Geophys. Res. Lett. 2009, 36, L20315. [Google Scholar] [CrossRef]
  60. Deng, K.; Zhou, S.; Wang, R.; Robinson, R.; Zhao, C.; Cheng, W. Evidence that the 2008 Mw 7.9 Wenchuan Earthquake Could Not Have Been Induced by the Zipingpu Reservoir. Bull. Seismol. Soc. Am. 2010, 100, 2805–2814. [Google Scholar] [CrossRef] [Green Version]
  61. Segall, P. Earthquakes triggered by fluid extraction. Geology 1989, 17, 942–946. [Google Scholar] [CrossRef]
  62. González, P.J.; Tiampo, K.F.; Palano, M.; Cannavó, F.; Fernández, J. The 2011 Lorca earthquake slip distribution controlled by groundwater crustal unloading. Nat. Geosci. 2012, 5, 821. [Google Scholar] [CrossRef] [Green Version]
  63. Avouac, J. Human-induced shaking. Nat. Geosci. 2012, 5, 763–764. [Google Scholar] [CrossRef]
  64. López-Comino, J.A.; Mancilla, F.; Morales, J.; Stich, D. Rupture directivity of the 2011, Mw 5.2 Lorca earthquake (Spain). Geophys. Res. Lett. 2012, 39, L03301. [Google Scholar] [CrossRef]
  65. Kundu, B.; Vissa, N.; Gahalaut, V. Influence of anthropogenic groundwater unloading in Indo-Gangetic plains on the 25 April 2015 Mw 7.8 Gorkha, Nepal earthquake. Geophys. Res. Lett. 2015, 42, 10607–10613. [Google Scholar] [CrossRef] [Green Version]
  66. Amos, C.; Audet, P.; Hammond, W.; Bürgmann, R.; Johanson, I.; Blewitt, G. Uplift and seismicity driven by groundwater depletion in central California. Nature 2014, 509, 483. [Google Scholar] [CrossRef] [PubMed]
  67. Biot, M.A. General theory of three-dimensional consolidation. J. Appl. Phys. 1941, 12, 155–164. [Google Scholar] [CrossRef]
  68. Cheng, A.H. Poroelasticity; Springer International Publishing: Cham, Switzerland, 2016. [Google Scholar]
  69. Wang, H.F. Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology; Princeton University Press: Princeton, NJ, USA, 2000. [Google Scholar]
  70. Rice, J.R.; Cleary, M.P. Some Basic Stress Diffusion Solutions for Fluid-Saturated Elastic Porous Media With Compressible Constituents. Rev. Geophys. Space Phys. 1976, 14, 227–241. [Google Scholar] [CrossRef]
  71. Segall, P. Earthquake and Volcano Deformation; Princeton University Press: Princeton, NJ, USA, 2010. [Google Scholar]
  72. Merxhani, A. An introduction to linear poroelasticity. arXiv 2016, arXiv:physics.geo-ph/1607.04274. [Google Scholar]
  73. Coussy, O. Poromechanics; Wiley: New York, NY, USA, 2004. [Google Scholar]
  74. Verruijt, A. Theory and Problems of Poroelasticity; Delft University of Technology: Delft, The Netherlands, 2016. [Google Scholar]
  75. Ghanimi, M.A.; Leroy, Y.; Kamp, A. Stress-sensitive permeability: Application to fault integrity during gas production. Transp. Porous Media 2017, 118, 345–371. [Google Scholar] [CrossRef]
  76. Talwani, P.; Chen, L.; Gahalaut, K. Seismogenic permeability, ks. J. Geophys. Res. Solid Earth 2007, 112. [Google Scholar] [CrossRef]
Figure 1. Sketch of the model geometry and computational meshes. We study the Coulomb Failure Function (CFF) changes in a vertical strike-slip fault with depth produced by both the tidal effects on a two-dimensional (2D) model (work plane on panel (a)), and reservoir operation (panel (c)). In the 2D case, the fault is embedded in a 5 km side square domain, while, in the three-dimensional (3D) one, it is embedded in a 20 km side cube. We model the fault as a rectangle with a certain width, softer and more permeable than the host rock matrix. The computational mesh of the 2D case is shown in panel (b) and the 3D case in panel (d).
Figure 1. Sketch of the model geometry and computational meshes. We study the Coulomb Failure Function (CFF) changes in a vertical strike-slip fault with depth produced by both the tidal effects on a two-dimensional (2D) model (work plane on panel (a)), and reservoir operation (panel (c)). In the 2D case, the fault is embedded in a 5 km side square domain, while, in the three-dimensional (3D) one, it is embedded in a 20 km side cube. We model the fault as a rectangle with a certain width, softer and more permeable than the host rock matrix. The computational mesh of the 2D case is shown in panel (b) and the 3D case in panel (d).
Water 12 02724 g001
Figure 2. Close-up view (50 m by 150 m) of overpressure near the surface, illustrating the influence of fault permeability on the pore pressure after 63 h (high tide situation), using the Darcy simplification in a 2D model. We represent the increment in the pore pressure for values between 0 and 30 kPa, from warm to cold colors respectively. We have used diverse values of the fault permeability ( k F = 10 15 , 10 13 , 10 11 and 10 9 m 2 ) while keeping the host rock permeability constant ( k = 10 15 m 2 ). The parameters of the simulations related to this Figure are detailed on Table 1. The compressibility of the solid grains, χ s , is 2.33 · 10 10 Pa 1 . The permeability of the fault is: (a) ( k F = 10 15 m 2 ), (b) k F = 10 13 m 2 , (c) k F = 10 11 m 2 , (d) k F = 10 9 m 2 .
Figure 2. Close-up view (50 m by 150 m) of overpressure near the surface, illustrating the influence of fault permeability on the pore pressure after 63 h (high tide situation), using the Darcy simplification in a 2D model. We represent the increment in the pore pressure for values between 0 and 30 kPa, from warm to cold colors respectively. We have used diverse values of the fault permeability ( k F = 10 15 , 10 13 , 10 11 and 10 9 m 2 ) while keeping the host rock permeability constant ( k = 10 15 m 2 ). The parameters of the simulations related to this Figure are detailed on Table 1. The compressibility of the solid grains, χ s , is 2.33 · 10 10 Pa 1 . The permeability of the fault is: (a) ( k F = 10 15 m 2 ), (b) k F = 10 13 m 2 , (c) k F = 10 11 m 2 , (d) k F = 10 9 m 2 .
Water 12 02724 g002
Figure 3. Plot of the evolution of pore pressure increment along the fault zone during a single tidal cycle. We perform 2D models with the Darcy law. We plot the pore pressure changes up to 5 km deep for a 3 m amplitude and 12 h tidal cycle. We consider two values of the fault permeability, k F , ((a), k F = 10 11 m 2 , and (b), k F = 10 9 m 2 ). Pore pressure changes are referred to the maximum value of the surface pore pressure, which corresponds to the high tide situation. We represent the pore pressure increment for several time instants during the tidal cycle. The colors vary from cold (early times) to warm (late times) as the simulation time advances. The color map and range used in panel (b) inset is the same for both panels. The parameters of the simulations related to this Figure are listed in Table 1. The compressibility of the solid grains, χ s , is 2.33 · 10 10 Pa 1 . The permeability of the fault is: (a) k F = 10 11 m 2 , (b) k F = 10 9 m 2 .
Figure 3. Plot of the evolution of pore pressure increment along the fault zone during a single tidal cycle. We perform 2D models with the Darcy law. We plot the pore pressure changes up to 5 km deep for a 3 m amplitude and 12 h tidal cycle. We consider two values of the fault permeability, k F , ((a), k F = 10 11 m 2 , and (b), k F = 10 9 m 2 ). Pore pressure changes are referred to the maximum value of the surface pore pressure, which corresponds to the high tide situation. We represent the pore pressure increment for several time instants during the tidal cycle. The colors vary from cold (early times) to warm (late times) as the simulation time advances. The color map and range used in panel (b) inset is the same for both panels. The parameters of the simulations related to this Figure are listed in Table 1. The compressibility of the solid grains, χ s , is 2.33 · 10 10 Pa 1 . The permeability of the fault is: (a) k F = 10 11 m 2 , (b) k F = 10 9 m 2 .
Water 12 02724 g003
Figure 4. Effect of the poroelastic coupling on pore pressure diffusion along the fault for the tidal problem. We study the relation between the fault permeability and the relative pore pressure increment at 5 km deep. On panel (a), we plot the depth at which arrives the 10% of the surface pressure after 24 tidal cycles for several values of the fault permeability, according to Darcy’s law. As the fault length is 5 km, there is a threshold in the fault permeability for which the pressure saturates. The parameters of the simulations depicted on panel (a) are detailed on Table 1. The compressibility of the solid grains, χ s , is 2.33 · 10 10 Pa 1 . On panel (b) we plot the fault permeability dependence on the pore pressure increment for both a Darcy and a poroelastic case, (dots in blue and red, respectively). In both cases, we refer the pore pressure changes to the maximum value of the surface pore pressure, which corresponds to the high tide situation. The parameters of the simulations depicted with the blue dots on panel (b) are listed on Table 1, and with the red dots on Table 2. The Young modulus of the rock and the fault are 40 and 1 GPa, respectively; the intrinsic rock permeability is 10 15 m 2 ; the intrinsic fault permeability varies between 10 15 and 5 · 10 8 m 2 ; and, the Biot coefficient, α B , is 0.8.
Figure 4. Effect of the poroelastic coupling on pore pressure diffusion along the fault for the tidal problem. We study the relation between the fault permeability and the relative pore pressure increment at 5 km deep. On panel (a), we plot the depth at which arrives the 10% of the surface pressure after 24 tidal cycles for several values of the fault permeability, according to Darcy’s law. As the fault length is 5 km, there is a threshold in the fault permeability for which the pressure saturates. The parameters of the simulations depicted on panel (a) are detailed on Table 1. The compressibility of the solid grains, χ s , is 2.33 · 10 10 Pa 1 . On panel (b) we plot the fault permeability dependence on the pore pressure increment for both a Darcy and a poroelastic case, (dots in blue and red, respectively). In both cases, we refer the pore pressure changes to the maximum value of the surface pore pressure, which corresponds to the high tide situation. The parameters of the simulations depicted with the blue dots on panel (b) are listed on Table 1, and with the red dots on Table 2. The Young modulus of the rock and the fault are 40 and 1 GPa, respectively; the intrinsic rock permeability is 10 15 m 2 ; the intrinsic fault permeability varies between 10 15 and 5 · 10 8 m 2 ; and, the Biot coefficient, α B , is 0.8.
Water 12 02724 g004
Figure 5. Effects of pore pressure diffusicity. Here we plot the depth at which arrives the 10% of the maximum surface pressure for several values of the fault permeability and fault diffusion coefficient. The fault permeability ranges from 10 8 to 10 12 m 2 , while the diffusion coefficient varies between 7 and 15. We represent the results with regard to the square root of D. Table 1 lists the parameters of the simulations.
Figure 5. Effects of pore pressure diffusicity. Here we plot the depth at which arrives the 10% of the maximum surface pressure for several values of the fault permeability and fault diffusion coefficient. The fault permeability ranges from 10 8 to 10 12 m 2 , while the diffusion coefficient varies between 7 and 15. We represent the results with regard to the square root of D. Table 1 lists the parameters of the simulations.
Water 12 02724 g005
Figure 6. Impact of both the poroelastic response and specific storage on pore pressure diffusion through the fault zone. Here we represent the effect of the Biot coefficient ( α B = 0.05 , 0.3 , 0.7 , and 0.95 ) on the relative pressure at 5 km deep (referred to the maximum value of the surface pore pressure). We consider two values of the fault permeability ( k F = 10 12 and 10 9 m 2 ) and diverse values of the Young modulus for the ground (E = 5 , 20 and 40 GPa). The rock stiffness has little effect on the pore pressure diffusion throughout the fault. On the other hand, the more impermeable the fault, the bigger the poroelastic effect. The Young modulus of the fault is 5 GPa and rock intrinsic permeability is 10 16 m 2 . Table 2 details the remaining parameters of the simulations.
Figure 6. Impact of both the poroelastic response and specific storage on pore pressure diffusion through the fault zone. Here we represent the effect of the Biot coefficient ( α B = 0.05 , 0.3 , 0.7 , and 0.95 ) on the relative pressure at 5 km deep (referred to the maximum value of the surface pore pressure). We consider two values of the fault permeability ( k F = 10 12 and 10 9 m 2 ) and diverse values of the Young modulus for the ground (E = 5 , 20 and 40 GPa). The rock stiffness has little effect on the pore pressure diffusion throughout the fault. On the other hand, the more impermeable the fault, the bigger the poroelastic effect. The Young modulus of the fault is 5 GPa and rock intrinsic permeability is 10 16 m 2 . Table 2 details the remaining parameters of the simulations.
Water 12 02724 g006
Figure 7. Time evolution of Coulomb Failure Function for two tidal cases at 5 km depth. Here we represent the CFF changes, as well as the resistance and the shear stress changes. As the latter are almost negligible, Δ CFF and resistance changes are quite similar. Positive values indicate increase in likelihood of failure, while negative values favor fault stability. A higher value of the fault permeability attenuates the Δ CFF changes and increases the phase lag. On panels (a) and (b), the intrinsic fault permeability is k F = 10 11 and 10 9 m 2 , respectively, and the rock intrinsic permeability is 10 15 m 2 . Table 2 lists the parameters of the simulations. The Young modulus of the rock and fault are, respectively, 40 and 1 GPa, and the Biot coefficient, α B , is 0.8, (a) k F = 10 11 m 2 , (b) k F = 10 9 m 2 .
Figure 7. Time evolution of Coulomb Failure Function for two tidal cases at 5 km depth. Here we represent the CFF changes, as well as the resistance and the shear stress changes. As the latter are almost negligible, Δ CFF and resistance changes are quite similar. Positive values indicate increase in likelihood of failure, while negative values favor fault stability. A higher value of the fault permeability attenuates the Δ CFF changes and increases the phase lag. On panels (a) and (b), the intrinsic fault permeability is k F = 10 11 and 10 9 m 2 , respectively, and the rock intrinsic permeability is 10 15 m 2 . Table 2 lists the parameters of the simulations. The Young modulus of the rock and fault are, respectively, 40 and 1 GPa, and the Biot coefficient, α B , is 0.8, (a) k F = 10 11 m 2 , (b) k F = 10 9 m 2 .
Water 12 02724 g007
Figure 8. Influence of the host rock matrix permeability on the pore pressure after 63 hours (high tide). Panels (a) and (b) show a general view (5 km by 5 km) of two models. Panels (cf) show a close-up view (50 m by 150 m) of the upper part of the fault. The permeability of the host rock ranges from k = 10 15 to k = 10 10 m 2 and the permeability of the fault is set to k F = 10 10 m 2 . The remaining hydro-mechanical properties are listed in Table 2. The Young modulus of the host rock and the fault is 40 GPa and the Biot coefficient is 0.8 . (a) k/ k F = 10 3 , (b) k/ k F = 1 , (c) k / k F = 10 5 , (d) k / k F = 10 4 , (e) k / k F = 10 3 , (f) k / k F = 0.1 .
Figure 8. Influence of the host rock matrix permeability on the pore pressure after 63 hours (high tide). Panels (a) and (b) show a general view (5 km by 5 km) of two models. Panels (cf) show a close-up view (50 m by 150 m) of the upper part of the fault. The permeability of the host rock ranges from k = 10 15 to k = 10 10 m 2 and the permeability of the fault is set to k F = 10 10 m 2 . The remaining hydro-mechanical properties are listed in Table 2. The Young modulus of the host rock and the fault is 40 GPa and the Biot coefficient is 0.8 . (a) k/ k F = 10 3 , (b) k/ k F = 1 , (c) k / k F = 10 5 , (d) k / k F = 10 4 , (e) k / k F = 10 3 , (f) k / k F = 0.1 .
Water 12 02724 g008
Figure 9. Effect of the host rock matrix permeability and poroelastic coupling on the pore pressure diffusion. We study the relation between both the rock matrix permeability and the Biot coefficient on the pore pressure build-up at 5 km depth. In order to quantify the effect of the poroelastic coupling, we adopt two values of the fault permeability ( 10 9 and 10 10 m 2 ), and two values of the Biot coefficient ( α B = 0 and 0.8 ). In both cases, we compute the ratio of the maximum pore pressure increment at 5 km depth to the maximum pore pressure increment on the seabed (which corresponds to the high tide). For each permeability ratio, the poroelastic coupling accounts for the difference between the blue dots and green squares. The parameters of these simulations are listed in Table 2. The Young modulus of the host rock and the fault are both 40 GPa.
Figure 9. Effect of the host rock matrix permeability and poroelastic coupling on the pore pressure diffusion. We study the relation between both the rock matrix permeability and the Biot coefficient on the pore pressure build-up at 5 km depth. In order to quantify the effect of the poroelastic coupling, we adopt two values of the fault permeability ( 10 9 and 10 10 m 2 ), and two values of the Biot coefficient ( α B = 0 and 0.8 ). In both cases, we compute the ratio of the maximum pore pressure increment at 5 km depth to the maximum pore pressure increment on the seabed (which corresponds to the high tide). For each permeability ratio, the poroelastic coupling accounts for the difference between the blue dots and green squares. The parameters of these simulations are listed in Table 2. The Young modulus of the host rock and the fault are both 40 GPa.
Water 12 02724 g009
Figure 10. Effect of fault permeability on the pore pressure diffusion after 2.5 years of operation of an annual regulation reservoir. Here we show the contour plots of the pressure build-up on three faults with permeability values k F = 10 11 , 10 10 and 10 9 m 2 on panels (a)–(c), respectively. The resevoir level is at its maximum value after 2.5 years of operation. Table 3 lists the hydro-mechanical properties of the fault and the host rock. (a) k F = 10 11 m 2 , (b) k F = 10 10 m 2 , (c) k F = 10 9 m 2 .
Figure 10. Effect of fault permeability on the pore pressure diffusion after 2.5 years of operation of an annual regulation reservoir. Here we show the contour plots of the pressure build-up on three faults with permeability values k F = 10 11 , 10 10 and 10 9 m 2 on panels (a)–(c), respectively. The resevoir level is at its maximum value after 2.5 years of operation. Table 3 lists the hydro-mechanical properties of the fault and the host rock. (a) k F = 10 11 m 2 , (b) k F = 10 10 m 2 , (c) k F = 10 9 m 2 .
Water 12 02724 g010
Figure 11. Impact of the reservoir operation on the pore pressure diffusion. We plot the results for a reservoir with annual regulation in the left column, and for daily regulation in the right column. On panels (a,b), we represent the ratio of the pore pressure build-up at 1, 3, 5, and 20 km depth to the maximum pressure on the bed of the reservoir. On panels (c,d), we show the distribution of the pore pressure build-up on the fault when the water level in the reservoir is maximum. The hydro-mechanical properties of the solid and water are listed in Table 3.
Figure 11. Impact of the reservoir operation on the pore pressure diffusion. We plot the results for a reservoir with annual regulation in the left column, and for daily regulation in the right column. On panels (a,b), we represent the ratio of the pore pressure build-up at 1, 3, 5, and 20 km depth to the maximum pressure on the bed of the reservoir. On panels (c,d), we show the distribution of the pore pressure build-up on the fault when the water level in the reservoir is maximum. The hydro-mechanical properties of the solid and water are listed in Table 3.
Water 12 02724 g011
Figure 12. Effect of p/t on the pore pressure diffusion. We plot the results for a reservoir with annual regulation in the left column, and for daily regulation in the right column. On panels (a,b), we represent the evolution of p/t at 1, 3, 5, and 20 km depth. On panels (c,d), we depict p/t on the fault when p/t on the bed of the reservoir is maximum, i.e., the reservoir level is half. Table 3 lists the hydro-mechanical properties of the solid and the water.
Figure 12. Effect of p/t on the pore pressure diffusion. We plot the results for a reservoir with annual regulation in the left column, and for daily regulation in the right column. On panels (a,b), we represent the evolution of p/t at 1, 3, 5, and 20 km depth. On panels (c,d), we depict p/t on the fault when p/t on the bed of the reservoir is maximum, i.e., the reservoir level is half. Table 3 lists the hydro-mechanical properties of the solid and the water.
Water 12 02724 g012
Figure 13. Impact of the poroelastic response on the pore pressure attenuation (left column) and p/t (right column) at 3, 5 and 20 km depth. For panels (a) and (b) the fault permeability is k F = 10 11 m 2 , while for panels (c) and (d) is k F = 10 15 m 2 . When the fault is nearly impermeable, the poroelastic coupling is the main mechanism for the pore pressure diffusion in the fault. The remaining parameters of the simulations are detailed in Table 3. The regulation of the reservoir is annual.
Figure 13. Impact of the poroelastic response on the pore pressure attenuation (left column) and p/t (right column) at 3, 5 and 20 km depth. For panels (a) and (b) the fault permeability is k F = 10 11 m 2 , while for panels (c) and (d) is k F = 10 15 m 2 . When the fault is nearly impermeable, the poroelastic coupling is the main mechanism for the pore pressure diffusion in the fault. The remaining parameters of the simulations are detailed in Table 3. The regulation of the reservoir is annual.
Water 12 02724 g013aWater 12 02724 g013b
Table 1. Parameters of the two-dimensional (2D) model with Darcy’s law (Figure 2, Figure 3, Figure 4 and Figure 5).
Table 1. Parameters of the two-dimensional (2D) model with Darcy’s law (Figure 2, Figure 3, Figure 4 and Figure 5).
ParameterValueUnitDescription
L5kmDomain side
W10mFault width
A3mTidal amplitude
T12hTidal period
ρ f 1000kg/m 3 Density of the fluid
μ f 0.001Pa·sDynamic viscosity of the fluid
χ f 4·10 10 Pa 1 Compressibility of the fluid
k10 15 m 2 Intrinsic rock permeability
k F 10 15 –5·10 8 m 2 Intrinsic fault permeability
ϕ 0.1Porosity
Table 2. Parameters of the 2D poroelastic models (Figure 4b, Figure 6, Figure 7, Figure 8 and Figure 9).
Table 2. Parameters of the 2D poroelastic models (Figure 4b, Figure 6, Figure 7, Figure 8 and Figure 9).
ParameterValueUnitDescription
L5kmDomain side
W10mFault width
A3mTidal amplitude
T12hTidal period
ν 0.25Poisson ratio of the rock
ρ s 2500kg/m 3 Dry density of the solid grains of rock
ρ b 2350kg/m 3 Bulk density of the saturated rock
ρ f 1000kg/m 3 Density of the fluid
μ f 0.001Pa·sDynamic viscosity of the fluid
χ f 4·10 10 Pa 1 Compressibility of the fluid
ϕ 0.1Porosity
Table 3. Parameters of the three-dimensional (3D) models (Figure 10, Figure 11, Figure 12, Figure 13).
Table 3. Parameters of the three-dimensional (3D) models (Figure 10, Figure 11, Figure 12, Figure 13).
ParameterValueUnitDescription
L20kmDomain side
W100mFault width
E40GPaYoung modulus of the rock
E F 10GPaYoung modulus of the fault
ν 0.25Poisson ratio of the rock
ρ s 2500kg/m 3 Dry density of the solid grains of rock
ρ b 2350kg/m 3 Bulk density of the saturated rock
ρ f 1000kg/m 3 Density of the fluid
μ f 0.001Pa·sDynamic viscosity of the fluid
χ f 4·10 10 Pa 1 Compressibility of the fluid
k10 16 m 2 Intrinsic rock permeability
k F 10 9 , 10 10 , 10 11 m 2 Intrinsic fault permeability
ϕ 0.1Porosity
α B 0.95Biot coefficient

Share and Cite

MDPI and ACS Style

Pampillón, P.; Santillán, D.; Mosquera, J.C.; Cueto-Felgueroso, L. Geomechanical Constraints on Hydro-Seismicity: Tidal Forcing and Reservoir Operation. Water 2020, 12, 2724. https://doi.org/10.3390/w12102724

AMA Style

Pampillón P, Santillán D, Mosquera JC, Cueto-Felgueroso L. Geomechanical Constraints on Hydro-Seismicity: Tidal Forcing and Reservoir Operation. Water. 2020; 12(10):2724. https://doi.org/10.3390/w12102724

Chicago/Turabian Style

Pampillón, Pedro, David Santillán, Juan Carlos Mosquera, and Luis Cueto-Felgueroso. 2020. "Geomechanical Constraints on Hydro-Seismicity: Tidal Forcing and Reservoir Operation" Water 12, no. 10: 2724. https://doi.org/10.3390/w12102724

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop