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

: 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 ﬁnd that pressure diffusion to such depths is only possible for highly permeable fault zones and/or strong poroelastic coupling. Data Formal analysis, 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. L.C.-F.


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 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.

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.

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 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: 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]: where χ f is the compressibility of the fluid and χ s is the compressibility of the solid grains.

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]: 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]: 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]: 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]: where ν u in a fully saturated rock with a single solid constituent is given by: 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 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]: 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]: The above expression is expressed in terms of the bulk modulus K, as follows [73,74]: 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]:

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.   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).
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.

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 Tables 1 and 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 , where ρ f is the fluid density and h(t) is the tidal sea-level oscillation. We compute h(t) with a sinusoidal wave: where A is the amplitude of the tide and T the period.     Table 1. The compressibility of the solid grains, 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. 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.  Table 2 details the remaining parameters of the simulations.   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.
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.

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.    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 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. The regulation of the reservoir is annual.
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.

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.

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.

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.

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).

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 Figures 8 and 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.

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.

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.