**1. Introduction**

Experimental and numerical simulations show that disorder plays a key role in driving stress accumulation in the crust and energy nucleation during earthquakes [1], nevertheless it was not clarified ye<sup>t</sup> how stress variations trigger breaking processes in such heterogeneous media. Indeed, earthquakes can be due to several stress sources, such as magmatic

**Citation:** Zaccagnino, D.; Telesca, L.; Doglioni, C. Different fault response to stress during the seismic cycle. *Appl. Sci.* **2021**, *11*, 9596. https:// doi.org/10.3390/app11209596

Academic Editor: Igal M. Shohet

Received: 6 September 2021 Accepted: 11 October 2021 Published: 14 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

intrusion or overpressured liquids; moreover, faulting is also affected by temperature, confining and pore pressure, and rock brittleness. This is why the comprehension of the response of faulting to additional stress was so actively investigated for 50 years. There are a few exogenous stress sources useful for this purpose: fluid injection is a widespread technique in stimulating production from oil and natural gas wells and improve geothermal energy generation. Although it is usually associated with microseismicity, several events with moderate magnitudes were also related to this practice [2]. This is why it is of paramount importance to improve our knowledge about the conditions under which intermediate magnitude events might occur. Also, it is crucial to note that injection and depletion generally happen at different wells, leading to a complex underground liquid circulation. Thus, it is not so easy to model how fluid injection drives spatial variations in pore pressure. An additional source of complexity is due to the variability of the time interval between the beginning of fluid injections and the onset of seismic activity. Therefore, fluid injection cannot be considered an efficient way of monitoring fault response to stress perturbation, at least over the time interval we are interested in. For this reason, we do not focus on this kind of stress source. A second possibility may be the controlled use of explosive, a tested tool in engineering of rock blasting, drilling, and mining, and moreover it is at the base of field reflection and refraction seismology. Unfortunately, this method is useful only for stress pulses simulations.

On the contrary, lunar and solar tides continuously induce periodic deformations in solid earth. Tidal harmonics are featured by different frequencies, so that their periods range from 10<sup>4</sup> to 10<sup>9</sup> s. The displacement of the tidal bulge can be decomposed into its vertical and horizontal components that depend on latitude, and may amount respectively up to 40 cm and 20 cm in the case of the semidiurnal M2 tide. Solid Tides depend on depth reaching their highest intensity around 1000 km below the surface [3]. Despite the fact that solid and ocean tidal stress (0.1 kPa–100 kPa) is fairly smaller than the earthquake stress drops (1 MPa–30 MPa, [4]) it was sufficiently proven that tides can trigger earthquakes (e.g., [5–9]) even though global seismicity weakly correlates with the Moon and Sun's distances from Earth and long catalogs ( ≥10<sup>4</sup> events, [10]) are needed to detect this effect accurately. For these reasons, in this article we investigate how the response of faults to stress modulations changes during the "seismic cycle" using tidal perturbations. In particular, we focus on recent Greek seismicity.

Since the tidal bulge is misplaced relative to the gravitational Earth-Moon alignment, being about 0.3–2.4 degrees eastward of it [11,12] due to the delay in reaction for the anelastic component of the Earth as a response to the tidal pull, a westerly-directed horizontal drag acts on the lithosphere and the continuously slows down Earth rotation [13]. The crux of the matter is that, besides the symmetric oscillatory motions, the nonlinear response of the low velocity zone (LVZ) breaks the symmetry of the tidal traveling wave on the Earth surface: the small asymmetry produces a net drift motion of any material point interacting with the gravitational wave force, in the direction following the rotation of the Moon. Therefore, detailed properties of the combined tidal oscillation and tidal drift depend on the degree of deformability of the lithosphere (a.k.a. Love and Shida numbers) due to its local temperature and geochemistry. These modulations act through two different mechanisms on the outer layers of the Solid Earth: the brittle and outermost part of the crust is elastically affected by tidal waves, which induce a stress variation in rocks that, depending on the local geodynamics, promotes or, on the contrary, can prevent the achievement of the critical breaking point, so that it plays a statistically significant role in fault activation [14–16], while solid tides actively modulate plate motion at low frequencies over geological periods [17,18] due to a heterogeneous dissipation of the tidal torque at the LVZ level.

The issue of the tidal triggering of seismicity is not trivial, since it is necessary to take into account not only of the effect of Earth tides, but also variations in pore and confinement pressure, the orientation of the DC (Double Couple contribution) in relation with focal mechanism, and CLVD (Compensated Linear Vector Dipole) of earthquakes and local geophysical heterogeneities.

To make matters even more complicated, a further difficulty must be considered: the observed seismicity spans up to 11 orders of magnitude over time (from 10 s typical of microseismicity up to 10<sup>12</sup> s for earthquakes recurrence time intervals), and 14 concerning to energy (if one considers events ranging from M*w* 0 to 9.5) and seven in space (a M*w* 9.0 can breaks crust up to 1000 km); this creates an enormous obstacle of both resolution and saturation, catalog incompleteness [19], and unreliability of statistical data.

Finally, to the technical complexity of measuring real stress, we add the estimation of uncertainties of seismic parameters such as magnitude and depth.

Beyond the triggering mechanism, in this work we focus on the possibility of highlighting the growth of critical states in the crust induced by stress accumulation in rocks through the measurement of correlations between some features of the tidal perturbation and seismic activity.

There is a flurry of scientific articles devoted to understand how tides influence seismic activity, but only a few of them (e.g., [20]) so far extensively studied whether tidal perturbation might somehow provide information about the stability of faults close to rupture, with a few exceptions regarding some particular case studies (e.g., [21,22]). The novelty of this article with respect to previous scientific literature consists in highlighting not only the variations in responsiveness of faults to stress while reaching the critical point, i.e., before a large earthquake, but it also shows how seismic response of faults changes over time as a function of their instability.

#### **2. Materials and Methods**

The lunisolar tides deform the Earth up to 60 cm twice a day, moreover the weight of the ocean tides gives a periodic load on the Earth's surface strongly dependent on bathymetry [23]. Although the displacements are relatively large, the associated changes in strain at the Earth's surface are tiny, extremely difficult to measure accurately, and even more tricky to model. The main difference between liquid tides and solid tides is in the phase: rocks react quickly to solicitations, while fluid masses need a characteristic time span to move, so they are affected by the tide with a certain phase shift.

We model tidal stress according to the following method.

Considering two massive celestial bodies with a spherical distributed mass density one realizes that the gravitational and centrifugal forces are balanced whenever their volumes are deformed by tides.

$$
\vec{F}\_{\mathcal{M}}(P) = \vec{F}\_{\mathcal{G}}(P) + \vec{F}\_{\mathcal{C}}(P) \tag{1}
$$

from which the tidal acceleration is immediately obtained

$$|\vec{a}\_{\mathcal{M}}(R)| = \frac{GM}{(R-r)^2} - \frac{GM}{R^2} \simeq \frac{2MGr}{R^3} \tag{2}$$

in the last step *R r* is assumed. The gravitational potential due to celestial body with mass *M*, at a point *P* at a distance *r* from the center of the planet, can be expanded in a series of powers of *<sup>r</sup>*/*R*, where *R* is the distance from the center of the Earth to the celestial body.

$$\mathcal{W}(P) = V(P) + \vec{F}\_G(P) \cdot \vec{r}(P) + const \tag{3}$$

$$\mathcal{W}(\vec{\mathcal{R}}) = -\frac{GM}{|\vec{\mathcal{R}} + \vec{r}|} + \frac{GM}{R^3}\vec{\mathcal{R}} \cdot \vec{r} - \frac{GM}{R} \tag{4}$$

if *l* = |- *R* +*r*|

$$\mathcal{W}(\vec{\mathcal{R}}) = -\frac{GM}{R} \left( \frac{R}{l} - \frac{\vec{\mathcal{R}} \cdot \vec{r}}{R^2} - 1 \right) \tag{5}$$

by expanding in series of Legendre polynomials we obtain

$$W(R,\Psi) = \frac{GM}{R} \sum\_{n=2}^{\infty} \left(\frac{r}{R}\right)^n P\_n(\cos\Psi) \tag{6}$$

Since *r*/*R* for the Moon is ∼1/60 while for the Sun is ∼1/23,000, the contributions of successive terms to the potential rapidly decrease. For the Moon, ∼98% of the total tidal potential, and for the Sun, the higher orders are completely negligible for our purpose. Moreover, from the ratios of masses and their mean distances, it follows that the solar tidal perturbation is ∼0.459 times the lunar tides. So, we can write that the total tidal potential is given by the sum of lunar and solar perturbations in the following form

$$\mathcal{W}(R,\Psi) \simeq \frac{GM}{R} \left(\frac{r}{R}\right)^2 P\_2(\cos\Psi) \tag{7}$$

where Ψ is the zenith of the body with respect to P and *P*2 is the second degree Legendre polynomial.

$$
\cos\Psi = \cos\theta\cos\delta + \sin\theta\sin\delta\cos(\phi-\alpha) \tag{8}
$$

*θ* is the colatitude and *φ* is the easterly longitude of P, *δ* is the codeclination, and *α*is the right-ascension of the body. We can write the potential so that three different contributions are highlighted

$$\mathcal{W}(R,\Psi\_1,\Psi\_2,\Psi\_3) = \frac{3GMr^2}{4R^3} [\Psi\_1 + \Psi\_2 + \Psi\_3] \tag{9}$$

with

$$\Psi\_1 = 3\left(\sin^2\left(\frac{\pi}{2} - \theta\right) - \frac{1}{3}\right)\left(\sin^2\left(\frac{\pi}{2} - \delta\right) - \frac{1}{3}\right) \tag{10}$$

$$\Psi\_2 = \sin\left(2\left(\frac{\pi}{2} - \theta\right)\right) \sin\left(2\left(\frac{\pi}{2} - \delta\right)\right) \cos\left(\phi - \alpha\right) \tag{11}$$

$$\Psi\_3 = \cos^2\left(\frac{\pi}{2} - \theta\right) \cos^2\left(\frac{\pi}{2} - \delta\right) \cos\left(2\left(\phi - \alpha\right)\right) \tag{12}$$

For the sake of simplicity

$$D = \frac{3GMr^2}{4R^3}$$

is called Doodson's parameter.

The vertical displacement is obtained by dividing *W* for the local value of the gravitational acceleration

$$\delta\mathcal{H}(\mathbf{R}, \Psi\_1, \Psi\_2, \Psi\_3) = \frac{\mathcal{W}(\mathbf{R}, \Psi\_1, \Psi\_2, \Psi\_3)}{\mathcal{S}} = \frac{3GMr^2}{4gR^3} \left(\Psi\_1 + \Psi\_2 + \Psi\_3\right) \tag{13}$$

The three terms represent the zonal, tesseral, and sectoral tides respectively. The zonal and sectoral contributes are responsible for tides with a half the Moon's revolution period, while the tesseral one for planet's rotation period tides. A paramount research work in this field was conducted by A. T. Doodson (1890–1968), who identified 378 tidal harmonics that were collected in a celebrated catalog in 1921 [24]. The largest tidal harmonics [25] are shown in Table 1. The elastic deformation of the Earth has two modes: spheroidal and toroidal, but tidal forces excite only the spheroidal modes [26].



For the sake of simplicity, we assume that the Earth is spherically symmetric, nonrotating, elastic, and isotropic. Nevertheless, the effect of sphericity of the Earth's layering cannot be neglected if one wishes to calculate surface waves of long wavelength. Spheroidal deformations of the SNREI model as a function of depth, Lamé's coefficients and the gravitational acceleration can be evaluated by using a set of functions *yi* with *i* = 1, . . . , 6 which satisfy a set of six ordinary differential equations [27]

$$\frac{dy\_i(r,n)}{dr} = \sum\_{j=1}^6 f(\rho(r), \lambda(r), \mu(r), \lg(r), n) y\_j(r,n),\tag{14}$$

where *n* is the *n*-th mode. For our purpose, it is enough to consider *n* = 2. Displacements can be written in spherical coordinates as follows:

$$\begin{cases} u\_r = \frac{l\_2(r)}{\mathcal{S}(r)} \mathcal{W}(r, \theta, \phi) \\ u\_\theta = \frac{l\_2'(r)}{\mathcal{S}(r)} \frac{\partial \mathcal{W}(r, \theta, \phi)}{\partial \theta} \\ u\_\phi = \frac{l\_2(r)}{\mathcal{S}(r)\sin\theta} \frac{\partial \mathcal{W}(r, \theta, \phi)}{\partial \phi} \end{cases} \tag{15}$$

The value of *g*(*r*) rises from the surface up to a depth of about 700 km to a maximum of 9.99 m/s2. In the lower mantle, *g*(*r*) lingers stable and increases abruptly near the Gutenberg discontinuity, reaching 10.16 m/s2. Gravity continuously falls in the core with a rate that depends on density till it reaches a zero-value at the center of the Earth.

In this work, tidal perturbations are considered with respect to seismicity and, above all, crustal seismicity (up to ∼100 km at depth), so that if we assume *g*(*r*) *g*(*r* = 6371 km) an error ∼0.1% is introduced, which is negligible. The strain components in spherical coordinates are obtained from (15) by derivation

$$\begin{cases} \varepsilon\_{rr} = \frac{\partial u\_r}{\partial r} \\ \varepsilon\_{\theta\theta} = \frac{1}{r} \left( \frac{\partial u\_\theta}{\partial \theta} + u\_r \right) \\ \varepsilon\_{\theta\phi} = \frac{1}{r \sin \theta} \left( \frac{\partial u\_\theta}{\partial \phi} + u\_r \sin \theta + u\_\theta \cos \theta \right) \\ \varepsilon\_{\theta\phi} = \frac{1}{r} \left( \frac{1}{\sin \theta} \frac{\partial u\_\theta}{\partial \phi} + \frac{\partial u\_\phi}{\partial \theta} - u\_\phi \cot \theta \right) \end{cases} \tag{16}$$

For computational simplicity, the radial component of strain can also be evaluated by using

$$
\varepsilon\_{\mathcal{II}} = -\frac{\nu}{1-\nu} (\varepsilon\_{\theta\theta} + \varepsilon\_{\phi\phi}) \tag{17}
$$

where *ν* is the Poisson's coefficient that can be computed at depth with

$$\nu(r) = \frac{\lambda(r)}{2(\lambda(r) + \mu(r))}\tag{18}$$

in turn, the Lame's coefficients can be obtained starting from the speed of the seismic P and S waves as a function of depth (PREM, [28]).

$$\begin{cases} \lambda(r) = \rho(r) \left( v\_P^2(r) - v\_S^2(r) \right) \\ \mu(r) = \rho(r) v\_S^2(r) \end{cases} \tag{19}$$

So, the needed components of stress in spherical coordinates are

$$\begin{cases} \sigma\_{\theta\theta}(r) = \lambda(r) \left( \varepsilon\_{rr} + \varepsilon\_{\phi\phi} + \varepsilon\_{\theta\theta} \right) + 2\mu(r)\varepsilon\_{\theta\theta} \\ \sigma\_{\phi\phi}(r) = \lambda(r) \left( \varepsilon\_{rr} + \varepsilon\_{\phi\phi} + \varepsilon\_{\theta\theta} \right) + 2\mu(r)\varepsilon\_{\phi\phi} \\ \sigma\_{\theta\phi}(r) = \mu(r)\varepsilon\_{\theta\phi} . \end{cases} \tag{20}$$

At last, it is necessary to take into account the spatial orientation of faults, which provides information about the tectonic stress tensor. Given the strike *α* of the seismological source, the tangential stress is

$$
\sigma\_a^{(\pm)} = \sigma\_{\theta\theta}(r) \cos^2 \alpha + \sigma\_{\phi\phi}(r) \sin^2 \alpha \pm 2\sigma\_{\theta\phi}(r) \sin \alpha \cos \alpha. \tag{21}
$$

Then, the quantities of geophysical interest are the following

$$\begin{cases} \sigma\_{\delta} = \sigma\_{\alpha}^{(+)} \sin \delta \cos \delta \\ \sigma\_{\hbar} = \sigma\_{\alpha}^{(+)} \sin^{2} \delta \\ \sigma\_{\delta} = \frac{1}{3} \left( \sigma\_{\alpha}^{(+)} + \sigma\_{\alpha}^{(-)} \right) \end{cases} \tag{22}$$

where *δ* is the dip angle of the fault. *α* and *δ* are inferred by comparing the focal mechanisms of local seismicity with the maps of actual faults. Since focal mechanisms are calculated only for earthquakes with significant magnitude, usually larger than 3.5–4.0, routinely recorded small events are assumed to occur on fault planes whose angles of strike and dip are given by the averages of the available ones.

*σs* is the so-called shear stress, which is positive in extensional tectonics, *σn* is the normal stress acting orthogonally to the fault plane, and *σc* is the confining stress due to the weight of the overlying rocks and fluids [29]. *σrr*, *σrθ*, and *<sup>σ</sup>rφ* are not considered since they are negligible up to 300 km deep [3] and ∼95% of the seismic energy is nucleated within the depth 0–50 km. For these reasons, only the horizontal shear stresses *σθθ* and *σφφ* can effectively play a role in triggering earthquakes.

A still open problem concerns the functions *yi*: for the calculations above, only *y*1 and *y*3 are needed since

$$\begin{cases} y\_1(r,2) = \frac{h\_2(r)}{\mathcal{S}(r)}\\ y\_3(r,2) = \frac{f\_2'(r)}{\mathcal{S}(r)} \end{cases} \tag{23}$$

they can be obtained by integrating with the fourth order Runge–Kutta method, a system of six coupled ordinary differential equations starting from a set of suitable boundary conditions. In turn, these can be calculated by considering the solutions in the case of an isotropic and homogeneous sphere on the Gutenberg discontinuity. The result is a combination of spherical Bessel function of the first kind that can be computed with a power series expansion.

As regards ocean tides, a foreword is necessary.

Despite ocean loading being able to induce stress up to 100 kPa, which is much larger than the stress due to solid tides (0.1–3 kPa), it is locally generated and usually focused over small surfaces (-10<sup>4</sup> km2), with some exceptions such as off the coast of New Zealand, the Madagascar Channel, the Java–Timor Sea and offshore Alaska. In practice, the main contribution of oceanic tides derives, unlike solid tides, from vertical stress [30]

$$
\sigma\_{zz} = -\rho g h \tag{24}
$$

where *h* is the amplitude of the tide and *ρ* 1030 kg/m3; the radial stress spread horizontally through the Poisson's coefficient. So, working in local Cartesian coordinates, if I assume that the vertical stress acts symmetrically *<sup>σ</sup>xx*<sup>=</sup>*σyy* [31]

$$
\sigma\_{\text{xx}} = \frac{\nu}{1-\nu} \sigma\_{zz} \tag{25}
$$

therefore, comparing with Equation (22) we ge<sup>t</sup>

$$\begin{cases} \sigma\_s \approx 0\\ \sigma\_n \approx \sigma\_{xx} \cos^2 \delta + \sigma\_{yy} \sin^2 \delta\\ \sigma\_c = \frac{1}{3} \left( \sigma\_{xx} + \sigma\_{yy} + \sigma\_{zz} \right) \end{cases} \tag{26}$$

the predicted tidal height *H* can be provided by the NAO.99b software [32].

At last, it is convenient to introduce the Coulomb Failure Stress (CFS) [33]

$$\text{CFS} = |\sigma\_{ts}| - \mu(\sigma\_{tn} - p) - S\_0 \tag{27}$$

where *σts* and *σns* are, respectively, the shear and normal stress, *p* is the pore pressure and *S*0 stands for the coesion of rocks. It is often assumed that changes in *p* are proportional to the normal stress change across the fault plane, so that rescaling *μ*

$$\text{CFS} = \sigma\_{\text{fs}} - \mu \sigma\_{\text{tr}} - \text{S}\_0 \tag{28}$$

with *μ* ∼ 0.4–0.8.

> Since the original state of stress is unknown, the ΔCFS [34] is usually studied

$$
\Delta \text{CFS} = \sigma\_s + \mu \sigma\_n \tag{29}
$$

where *σs* is the change in shear stress on the fault plane induced by the tidal perturbation in the slip direction and *σn* is the tidal normal stress. Positive ΔCFS is associated with encouraged seismicity, while negative values produce stress shadow effects which inhibit slip [35].

Since tidal stress upon the fault is known in principle, we can perform a correlation analysis according to the following steps:


rameters. The dominant contribution comes from the strike and dip angle errors, so that

$$\begin{cases} \epsilon\_{shar} \simeq \sqrt{\epsilon\_+^2 \sin^2 \delta \cos^2 \delta + \epsilon\_\delta^2 \sigma\_+^2} \cos^2 \delta\\ \epsilon\_{normal} \simeq \sqrt{\epsilon\_+^2 \sin^4 \delta + \epsilon\_\delta^2 \sin^2 2\delta} \\ \epsilon\_{confinement} \simeq \frac{1}{3} \sqrt{\epsilon\_+^2 + \epsilon\_-^2} \approx 0.47 \epsilon\_+ \end{cases} \tag{30}$$

where replaces the usual symbol for standard deviation to avoid misunderstanding with stress components *σs*, *σc*, *σn* and *σ*(±) *α* .

+ and − are the uncertainties of the positive and negative tangential stresses given by

$$\varepsilon\_{+} = \varepsilon\_{-} \simeq \varepsilon\_{a} \sqrt{(2r\_{\theta\theta}\cos\theta\sin\theta)^{2} + (2r\_{\phi\theta}\cos\theta\sin\theta)^{2} + (2r\_{\theta\phi}\cos2\theta)^{2}} \tag{31}$$

• The correlation between the magnitude of the seismic events and the intensities of the tidal stress components acting on the fault is calculated over fixed time intervals Δ*t* according to the following formula:

$$\rho\_{t\_w}^{(j)} = \frac{\sum\_{i=1}^{N\_{l\_n}} \left(M\_{wi} - \overline{M\_{w}}\right) \left(\sigma\_i^{(j)} - \overline{\sigma^{(j)}}\right)}{\sqrt{\sum\_{i=1}^{N\_{l\_n}} \left(M\_{wi} - \overline{M\_{w}}\right)^2 \sum\_{k=1}^{N\_{l\_n}} \left(\sigma\_k^{(j)} - \overline{\sigma^{(j)}}\right)^2}} \tag{32}$$

where *j* = *s*, *n*, *c*, ΔCFS meaning respectively the shear, normal, confinement and ΔCFS components of stress; *Ntn* is the number of failures occurred during the *n*-th time step.

ΔCFS is calculated for each event occurred within the selected region whose magnitude is above the completeness magnitude. Since we are interested in understanding how sensitivity of faults to additional stress modulations changes during the seismic cycle, we calculate the correlation between ΔCFS and the nucleated seismic energy for several time intervals (about 20 in Figures 1–3). To do this, the number of events used for the calculation of the correlation must be large enough to suppress stochastic fluctuations (the number of earthquakes for each point in Figures 1–3 is >200). On the other hand, short time intervals are not suitable for our goal because tides with not negligible amplitudes have semiannual and yearly frequencies. Therefore, Δ*t* < 1 yr can affect the correlation value. Moreover, we are looking for slow processes of progressive destabilization of crustal volumes, then averaging does not cause information loss, but only noise attenuation.

• The uncertainty on the correlation index is obtained by propagation of tidal stress errors and magnitudes.
