Stress heterogeneities and failure mechanisms induced by temperature and pore-pressure increase in volcanic regions

Earth and Planetary Science Letters 525 (2019) 115765

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

Stress heterogeneities and failure mechanisms induced by temperature and pore-pressure increase in volcanic regions M.E. Belardinelli ∗ , M. Bonafede, M. Nespoli a r t i c l e

i n f o

Article history: Received 11 February 2019 Received in revised form 21 May 2019 Accepted 7 August 2019 Available online xxxx Editor: R. Bendick Dedicated to Enzo, mentor supporter and friend Keywords: thermo-poro elasticity degassing magma chamber mixed mode thrust-tensile fractures Coulomb failure stress pore pressure change

a b s t r a c t We study the strain and stress fields produced by temperature and pore pressure increases within and outside a Thermo-Poro-Elastic (TPE) inclusion (the source region), embedded within a medium (the matrix) in isothermal drained conditions. This model is suitable to describe a crustal region in a volcanic environment pervaded by hot pressurized fluids released by an underlying magma chamber. After introducing the pertinent constitutive relations, a formal solution for the displacement field is provided in terms of the Green’s function for an elastic medium with drained isothermal elastic moduli, employing a generalization of Eshelby (1961) procedure. If an unbounded medium is considered, a displacement potential can be introduced, obeying the Laplace equation within the source region and the Poisson equation within the matrix. If a spherically symmetric source region is considered, simple analytical solutions are obtained for the displacement, the strain and the stress fields, showing that thrust faulting mechanisms are promoted within the source region while normal faulting mechanisms prevail in the embedding matrix. Employing reasonable numerical values for the thermo-poro-elastic parameters, suitable to describe highly porous sedimentary rock, strain and stress variations are found to be significant even for moderate changes of temperature and pore pressure. Variations of the Coulomb failure function are high in the TPE region and are strongly dependent on the friction coefficient and pore pressure. Application of these results to the 1982-84 and 2011-13 unrest episodes at Campi Flegrei caldera (Italy) suggests that an oblique dike intrusion across a previously unfaulted TPE region took place with a mixed tensile-thrust dislocation mechanism in both events, as previously inferred from accurate inversion of geodetic data. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Stresses and strains within and around magma chambers are generally overlooked in mechanical models. Sometimes they are modeled by plastic or viscoelastic stresses (see e.g. Segall, 2016; Bonafede and Ferrari, 2009). In volcanic regions magma intrusion at shallow depth is generally accompanied by exsolution of volatiles (mostly H2 O and CO2 ) at high temperature and pressure which can efficiently migrate across fissures and cracks of the embedding rock, bringing sharp and significant increases of temperature  T and of pore pressure  p. The increasing pore pressure, in turn, is expected to open pre-existing fissures and cracks wider, thus providing an enhanced effective permeability, as modeled e.g. by Zencher et al. (2006). Temperature and pore pressure changes provide significant stress and deformation and it is then necessary to include thermoporo-elasticity in the set of governing equations in order to eval-

*

Corresponding author. E-mail address: [email protected] (M.E. Belardinelli).

https://doi.org/10.1016/j.epsl.2019.115765 0012-821X/© 2019 Elsevier B.V. All rights reserved.

uate the impact of exsolved fluids and possibly discriminate their role from that due to magma intrusion. The ground deformation field observed at Campi Flegrei, a restless Caldera near Naples (Italy) during two major uplift episodes in 1970 and 1982-84 was interpreted in terms of the joint effects of overpressure in a shallow magma chamber and migration of hot pressurized fluids (Bonafede, 1990). After these major inflation episodes, totaling more than 3.5 m uplift, the caldera entered a phase of subsidence accompanied by minor unrest episodes, generally ascribed to fluid migration episodes (D’Auria et al., 2011). Chiodini et al. (2003) and Todesco and Berrino (2005) model geochemical anomalies, ground deformation and gravity changes in the Campi Flegrei caldera in terms of effects induced by heated and pressurized fluids released by a deeper magmatic source and circulating in the shallow hydrothermal system. Lima et al. (2009) suggest that short time-scale uplift and deflation episodes occurring at Campi Flegrei Caldera since 1970 can be modeled in terms of the episodic expulsion of aqueous fluids (released during magma solidification) from a deep lithostatically-pressurized reservoir toward an overlying hydrostatic reservoir.

2

M.E. Belardinelli et al. / Earth and Planetary Science Letters 525 (2019) 115765

2. Constitutive relations The strain field i j in a thermo-poro-elastic (TPE) medium which undergoes a change of stress τi j , a change of pore pressure  p and a change of temperature  T is provided by the following constitutive relation (e.g. McTigue, 1986):

i j =

Fig. 1. CGPS baseline length variation (cm) between two benchmarks encompassing the maximum uplift area and geochemical signal (arbitrary units) at two fumaroles in the Solfatara crater. Redrawn from Chiodini et al. (2015). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)



1 2μ



ν

τi j −

1+ν

τkk δi j +

1 3H

1

 p δi j + α  T δi j

(1)

3

where μ is the rigidity of the medium, ν the drained isothermal Poisson ratio, H is Biot’s constant and α is the coefficient of thermal expansion. As shown by Rice and Cleary (1976) H may 2μ(1+ν ) be computed from the bulk modulus K = 3(1−2ν ) of the porous medium in drained isothermal conditions and the bulk modulus K s pertinent to the solid phase in absence of interconnected porosity, according to the relation

1 H

=

1 K

−

1 K s

Separating the strain and stress tensors into their isotropic and deviatoric components, designated respectively with superfix ( I ) and ( D ), we have

1

1

3

3

i(jI ) ≡ kk δi j = Fig. 2. Conceptual model of the Thermo-Poro-Elastic source (yellow), affected by temperature and pore-pressure changes due to release of volatiles from a deeper magmatic reservoir (orange).

More recently, during 2005-2014, complex inflation episodes (Fig. 1) took place: a longer time, accelerating uplift signal (∼ 20 cm) was interpreted in terms of rock heating caused by enrichment in water of the magmatic fluids and by an increment in their flux, as inferred by temperature-pressure geo-indicators (Chiodini et al., 2015). Onto this long-term geodetic signal, short term pulses overlap, interpreted in terms of the injection of magmatic fluids into the hydrothermal system, as inferred from geochemical analyses of fumarolic gases, which follow the geodetic signal with a time delay ∼ 200 days. The purpose of this paper is to provide a mathematical model of the strain and stress field accompanying the intrusion of hot and pressurized fluids within a porous or fractured crustal volume. The source region is modeled as a thermo-poro elastic (TPE) medium, undergoing a sudden increase of temperature  T and pore pressure  p (Fig. 2). The source region is assumed to be surrounded by a medium in isothermal drained conditions, where both  T and  p vanish. Similar problems have been already considered (e.g. Bonafede, 1990, 1991; Rinaldi et al., 2010) with the aim of describing the surface displacement provided by a TPE source. Here we intend to model the deformation and the stress in the near field and even in the interior of the TPE source region, since intruding dikes (possibly leading to an eruption) are expected to be emplaced close to the inner caldera boundary or even within it. We shall see that a strongly heterogeneous stress field is generated within and outside the TPE source, which explains the extreme variability of seismic mechanisms determined at Campi Flegrei (47% normal, 23% inverse and 30% strike-slip) according to D’Auria et al. (2015) and provides support to the mixed tensile-thrust mechanism inferred from the inversion of geodetic data employing a realistic description of the heterogeneous caldera subsoil (Trasatti et al., 2011, 2015). Even if the proposed model was inspired by several observations made in the Campi Flegrei caldera, it is not specifically designed for this area and, thanks to its structural and geometrical simplicity, it may apply, at least in principle, to many other geothermal and volcanic regions of the world.



1 3K

τkk +

1

1

3

2μ

i(jD ) ≡ i j − kk δi j =



1 1 − 2ν

H

 + α  T δi j ,

τi(jD )

The inverse relation providing simply obtained from (1) as

τ i j = 2μ  i j +

p

(2)

τi j in terms of i j ,  p and  T is





νkk − (1 + ν )

p 3H

1

+ αT 3



 δi j (3)

The previous relations show that the effective bulk modulus ¯

K e f f = − V ∂∂ VP  13 τkk /kk (where P¯ = − 13 τkk denotes the average confining pressure, not to be confused with pore pressure p) depends on temperature and pore pressure changes, while the effective rigidity (the ratio of shear stress to twice the shear strain) is insensitive to them. From the relation (3) the isotropic and deviatoric stress components are:



1

τi(jI ) ≡ τkk δi j = K kk − 3

K H

  p − K α  T δi j ,

1

τi(jD ) ≡ τi j − τkk δi j = 2μi(jD ) 3

(4)

3. Stress free strain due to  T and  p Consider a bounded TPE region V s , undergoing a change of temperature and pore pressure, surrounded by an unbounded medium in drained and isothermal conditions (Fig. 2). In order to compute the stress and strain field induced by such changes within and outside V s we generalize the conceptual method devised by Eshelby (1961) and described in Aki and Richards (1980) for a purely elastic medium: we first remove the source volume V s (the inclusion) from the embedding material (the matrix) by cutting along the surface enclosing V s while keeping both the (0) inclusion and the matrix subjected to the same tractions τi j n j which were acting over before the cut (step 1 in Fig. 3); to be more specific, we denote the surface delimiting the matrix as + and the surface containing the inclusion, with opposite orientation, as − .

M.E. Belardinelli et al. / Earth and Planetary Science Letters 525 (2019) 115765

3

Fig. 3. Graphical illustration of the four conceptual steps devised by Eshelby (1961). The TPE region is yellow. The outward normal to the surface bounding the matrix is n+ , the outward normal to the surface bounding the inclusion is n− .

According to Kirchoff uniqueness theorem the initial shape and volume is preserved for both the inclusion and the matrix after step 1. Second, we increase the temperature and pressure within V s by  T and  p, leaving the tractions over + and − unchanged (step 2 in Fig. 3). From (1) we have the following “stressfree strain” within V s :

 i j =

1 3H

The displacement field produced by the traction discontinuity

Tk = −3K 0 nk acting on the surface elements d  located in x is given by (e.g. Aki and Richards, 1980)



u i (x) =

1

 p δi j + α  T δi j

with 0 =

1 p 3H

+ 13 α  T

(5)

The inclusion volume of course changes from V s to, say, V s ; if  p and  T are constant within the inclusion, V s = V s (1 + kk ) = V s (1 + 30 ). In step 3 we apply over − additional artificial tractions τi j n− keeping the pore pressure and the temperaj ture unchanged (isothermal drained conditions), which restore the initial shape and volume V s of the inclusion; clearly τi j n− = j

−C i jkl kl n−j = −3K 0 n− i , where C i jkl denote general elastic coef-

ficients and K = λ + 23 μ is the drained isothermal bulk modulus in an isotropic medium (for which C i jkl = λδi j δkl + μ(δik δ jl + δil δ jk ) with λ, μ Lamé constants in drained isothermal conditions). These tractions generate within the inclusion the artificial stress

τi j = −C i jkl kl = −3K 0 δi j Finally (step 4), the inclusion is inserted back into its original place: tractions over − differ by τi j n− with respect to + . j This amounts to a traction discontinuity +3K 0 n− between + i and − . When we remove this artificial discontinuity of tractions, (i.e. we add the traction discontinuity −3K 0 n− in isothermal and i drained conditions), the inclusion and the matrix adjust into a new equilibrium configuration (Fig. 3) with continuous tractions across . In the following sections we compute displacements, strain and stress fields caused by the application of this traction discontinuity describing the effects due to  T and  p changes within the region V S . Outside the source region V S , the medium is assumed in isothermal free drainage conditions; its constitutive relation is then perfectly elastic, employing drained, isothermal elastic parameters. According to the previous considerations, stresses in the source region V S and in the matrix are given, respectively, by

⎧

⎪ τi(jS ) ≡ 2μ i j + 1−12ν [νkk − (1 + ν ) 0 ] δi j ⎪ ⎪ ⎪ ⎪ ⎨ if x ∈ V s τi j (x) =

⎪ (M ) ν  δ ⎪ τ ≡ 2 μ  + ⎪ i j i j kk ⎪ ij 1 −2 ν ⎪ ⎩ otherwise

3K 0 nk G ik (x, x )d 

(7)



3

or, shortly

i j = 0 δi j

4. The deformation field provided by  T and  p

(6)

where G ik (x, x ) is the drained isothermal Green’s function i.e. the i-th component of displacement in x produced by a unit force applied in x in direction xk while keeping constant the pressure and the temperature. Employing Gauss theorem and assuming that 0 is constant, we may write



u i (x) = 3K 0 VS

and

∂ G ik (x, x )  dV ∂ xk

(8)

∂ G ik (x,x ) may be interpreted as the displacement in x pro∂ xk

duced by three orthogonal force dipoles with unit moment centered in x . It is worth mentioning that G ik (x, x ) is the same Green’s function as for a purely elastic medium, provided that the drained isothermal elastic constants are employed. Due to the singularity of the Green’s function within V S , the actual computation of the displacement field from eqs. (7)-(8) is generally awkward, and the computation of strain and stress components is even more so. 5. Unbounded medium

Significant simplifications are obtained if an unbounded medium is considered. We shall obtain equations formally similar to those governing the gravity field within and around an assigned density distribution, with the displacement field u replacing the gravity acceleration g. The Green’s function in a homogeneous isotropic elastic space is (e.g. Landau and Lifschitz, 1967):

G ik (x, x ) =

1 16πμ(1 − ν )



(3 − 4ν )

δik r

+

(xi − xi )(xk − xk )



r3

where r = [(x1 − x1 )2 + (x2 − x2 )2 + (x3 − x3 )2 ]1/2 = |x − x | is the distance between x and x . By straightforward calculation we obtain

xi − xi ∂ ∂ G ik (x, x ) 1 − 2ν 1 − 2ν = =  ∂ xk 8πμ(1 − ν ) r 3 8πμ(1 − ν ) ∂ xi

  1

(9)

r

Exploiting the property G ik (x, x ) = G ik (x − x ) (valid only for an unbounded medium) we may rewrite derivatives with respect to xi as derivatives with respect to xi with a change of sign, so that

∂ ∂ G ik (x, x ) 1 − 2ν =−  ∂ xk 8πμ(1 − ν ) ∂ xi

  1 r

=−

1 − 2ν 8πμ(1 − ν )

 

∇

1 r

(10)

4

M.E. Belardinelli et al. / Earth and Planetary Science Letters 525 (2019) 115765

and (8) may be rewritten in vector form as

u(x) = −0

(1 + ν ) ∇ 4π (1 − ν )

1

(x) = VS

|x − x  |

VS

1

|x − x  |

2μ(1+ν ) 3(1−2ν )

where the relation K =





dV  ,

dV  = −

so that

1 ∇ 4π

(1 + ν ) (1 − ν )



δ(x − x )dV  =



1 if x ∈ V S 0

otherwise

(12)

(13)

so that 1 is the dilatation within the TPE source. Accordingly, the displacement potential obeys the Poisson equation ∇ 2 S = −4π within the TPE source region V S and the Laplace equation ∇ 2 M = 0 in the matrix surrounding V S . Since displacements and tractions must be continuous across , the solutions S and M must match to each other properly over , as described below. In the embedding matrix, where kk = ∇ · u = 0 according to (13), and

(M )

τi j

∂ 2 M ∂ xi ∂ x j , we have from (6)

μ1 ∂ 2 M =− . 2π ∂ x i ∂ x j

Within the TPE source,



1 ∂ 2 S

(S)

τi j = −2μ1

(14)

kk = 1 and i j = − 4π1

4π ∂ x i ∂ x j

∂ 2 S ∂ 2 M − ∂ xi ∂ x j ∂ xi ∂ x j

 n j = 4π n i

(16)



This condition derives from the constraint of traction continuity across the inclusion surface and has no obvious counterpart in the gravity analogy mentioned at the beginning of the present section; it is equivalent to requiring that the linear dilatation/contraction i j ni n j in the normal direction to the boundary of the TPE region is discontinuous across the boundary, with jump amplitude −1 (contraction outside and dilatation inside along the normal direction). 6. Spherical source The geometric shape of the TPE domain might be very complex since it is strictly related to the spatial distribution of permeability and to the presence of fracture zones for which we usually have little - if any - “a priori” knowledge. Since our preliminary aim is mainly speculative, we shall assume spherically symmetric configurations. If the TPE inclusion V S is a sphere with radius a and center at the origin of an unbounded medium, invariance arguments suggest that u = u r (r ) rr is directed radially and its radial component u r depends only on the distance r from the origin. In spherical coordinates, inside V S we have 2

∇ S =

1 ∂ r2 ∂ r



r

2 ∂ S

∂r



= −4π

and the stress field is, from (15):

4

τi(jS ) = − μ1 δi j Outside V S , from ∇ 2 M = 0 we get similarly

M = −

C

+D

r

(M )

=⇒

ur

(r ) = −

1 C 4π r 2

and the inessential constant D may be set to zero in order that M vanishes at infinity; the constraint that u be continuous in r = a yields C = − 43 π a3 and continuity of is obtained if B = 2π a2 , so that we have finally

= 2π a 2 ×

(15)

Accordingly, the continuity of tractions τi j n j applied on the boundary of the source region V S requires a discontinuity in the second derivatives of the potential



1 ∂ S 1 = 1 r 4π ∂ r 3

⎧ ⎨ 1− ⎩

1 r2 3 a2

 if r < a

2a 3r

if r ≥ a

The previous equation and the following ones are written extracting the appropriate dimensional scale factor. The displacement field from (11) is

∂ 2 S ∂ xi ∂ x j , so that

 + δi j

(S)

u r (r ) = −

+B

3

VS

i j = − 4π1

r

where the integration constant A must vanish because the displacement must be bounded in r = 0 and B does not contribute to the displacement, but must be employed to connect continuously S with M in r = a. Within the source region, the displacement field from eq. (11) is

From this expression we immediately see that ∇ × u = 0 while, taking into account that ∇ 2 ( 1r ) = −4π δ(x − x ), we get

1 2 ∇ ·u=− ∇ = 1 4π

A

3

S (r ) = − π r 2 −

was employed and

1 = 0

2

(11)

1

u r (r ) =

3

1a ×

⎧ ⎨ ⎩

r a

if r < a 2

a r2

if r ≥ a

We observe in particular that, according to eq. (5), the stress-free strain i j (step 2 in Fig. 3) would inflate the sphere from the initial radius a to a0 = a(1 + 0 ), but the elastic constraint of the medium surrounding V S shrinks the sphere to the final radius  a1 = a + u (a) = a 1 + 13 1 ∼ a(1 + 59 0 ), the final estimate being obtained in the Poisson approximation ν = 1/4. The stress field outside the source domain is obtained from (14)

τi(jM ) =

2μ 3



1a3

δi j r3

−3

xi x j



r5

.

The stress components are more conveniently rewritten in spherical coordinates r , θ, ϕ ; exploiting the spherical symmetry of the problem we may anticipate that the tangential stresses τr θ , τr ϕ and τθ ϕ must vanish while τrr , τθθ and τϕϕ must be independent of θ and ϕ ; accordingly we may compute τrr (r ) as τxx (x = r , y = z = 0), τθθ (r ) = τzz (x = r , y = z = 0) and τϕϕ (r ) = τ y y (x = r , y = z = 0). We obtain explicitly



τrr (r ) =

− 43 μ1 − 43

if r < a

μ1 ×

and

τθ θ (r ) = τϕϕ (r ) =



 a 3 r

if r ≥ a

− 43 μ1 + 23 μ1 ×

if r < a

 a 3 r

if r > a

The previous solutions are illustrated graphically in Fig. 4. We observe in particular that τrr is everywhere compressive and is continuous across the surface of the sphere (as imposed according to

M.E. Belardinelli et al. / Earth and Planetary Science Letters 525 (2019) 115765

5

and from (14), (15):

4

τrr( S ) = − μ1 + 3

μ1 C ; π r3

τrr(M ) =

μ1 A π r3

Imposing the boundary condition τrr = − P 0 in r = a, the continuity of displacement in r = b and the continuity of we get



4

C=

3



P0

−

π a3 ,

μ1

A=−

4 3

π (b3 − a3 ) −

P0

μ1

π a3 ,

D = 2π b 2 .

Fig. 4. Normalized displacement potential × and non-vanishing stress components



1 2π a2



, radial displacement u ×



3



a1

3 τrr , τθ θ , τϕϕ × [ 4μ ] due to an increase of 1

temperature  T and pore pressure  p within a spherical TPE inclusion (yellow zone) embedded in an isothermal and drained unbounded matrix (gray). Arrows at the top describe the principal stress axes and the normal fault mechanisms expected within the matrix.

Eshelby method described in section 2) while τθθ and τϕϕ suffer a discontinuity from the TPE domain (compressional) to the matrix domain (extensional). We may easily check explicitly that the discontinuity condition (16) is fulfilled. The stress field is purely isotropic within the TPE domain and purely deviatoric in the embedding matrix. No shear stress is induced within the TPE region but, if the ground surface is envisaged laying horizontally above the sphere, normal faulting processes are favored above the TPE region where rˆ is vertical upwards. It is worthwhile pointing out that this solution is identical to the solution for a fluid sphere embedded within an elastic medium if an overpressure  P = −τrr (r = a) = 43 μ1 acts from within. 6.1. Spherical shell A formally similar problem possessing closed form solutions is the spherical TPE shell enclosing a spherical fluid-like reservoir with radius a at pressure P 0 , bounded by a spherical surface with radius b > a. The region r ≤ a may be considered as a magmatic fluid reservoir at pressure P 0 or, more generally, as a viscoelastic region surrounding a magmatic reservoir where, in the long term (after a few Maxwell times), an isotropic stress state τi j = − P 0 δi j is attained (see e.g. Dragoni and Magnanensi, 1989; Bonafede and Ferrari, 2009). The problem is again spherically symmetric and the equations to be solved are



1 ∂ r2

∂r

r2

∂ ∂r





if a ≤ r ≤ b if r > b

−4π

=

0

with boundary conditions prescribed over the spherical surface r = a and over r = b

τrr = − P 0 in r = a,

Continuity of τrr across r = b is automatically satisfied. In conclusion we obtain the following displacement (in the radial direction):

τrr and ur continuous across r = b

ur =

⎧   P 0 a3 1 a3 ⎪ (S) ⎪ u = +  r − ⎨ r 1 2 2

in a ≤ r ≤ b

⎪ ⎪ ⎩ ur( M ) =

in r > b

2

(r ) = S (r ) = − π r − 3 A

(r ) = M (r ) = −

r

C r

+ D (if a ≤ r ≤ b),

+ B (if r > b)

where A , B , C , D are constants of integration. According to eq. (11) the constants B and D do not contribute to the displacement field; B may be put equal to zero but D must be chosen in order to make continuous across r = b. We obtain from (11): (S)

ur

=

1 3

r−

1 C , 4π r 2

(M )

ur

=−

1 A 4π r 2

3 1

r (b3 − a3 )

3

r2

+ 1

4μr 2

(17)

and the following (non-vanishing) stress components:

τrr =

τθ θ

⎧   a3 4 a3 ⎪ ⎪ ⎨ τrr( S ) = − P 0 3 − μ1 1 − 3

in a ≤ r ≤ b

⎪ ⎪ ⎩ τrr( M ) = − P 0

in r > b

r a3 r3

3 4

r (b3 − a3 )

3

r3

− μ1

⎧   a3 4 a3 ⎪ (S) ⎪ ⎪ τθ θ = P 0 3 − μ1 1 + 3 ⎪ ⎪ 3 2r 2r ⎪ ⎪ ⎨ in a < r < b = τϕϕ = ⎪ ⎪ a3 2 (b3 − a3 ) ⎪ (M ) ⎪ τ = P + μ ⎪ 0 1 θθ ⎪ ⎪ 3 2r 3 r3 ⎩ in r > b

(18)

m The displacement um i and stress τi j which were present before degassing processes started are simply obtained posing 1 = 0 and P 0 = P m in the previous solutions (degassing processes plausibly are accompanied by a pressure drop from P m to P 0 < P m within the magmatic reservoir r < a). The incremental displacement and ( S /M ) ( S /M ) stress are u i − um and τi j − τimj ; they are simply obtained i substituting P 0 with ( P 0 − P m ) in the previous equations. In the simplest case when P 0 = P m the incremental fields are

ur =

⎧   1 r a2 ⎪ ⎪ ⎨ 1 a − 2 3 ⎪ ⎪1

a r (b3 − a3 ) a2

in a ≤ r ≤ b

⎩ 1 a in r > b 3 a3 r2 ⎧   3 ⎪ ⎪ − 4 μ1 1 − a ⎨ in a < r < b 3 3 r τrr = 3 3 3 ⎪ ⎪ ⎩ − 4 μ1 (b − a ) a in r > b 3 3 3

We obtain the following general solutions

2

4μr P 0 a3

τθ θ =

a

r

⎧   4 a3 ⎪ ⎪ ⎨ − μ1 1 + 3

in a < r < b

⎪ ⎪ ⎩ + μ1

in r > b

3 4 3

2r (b − a3 ) a3

(19)

3

a3

2r 3

which are shown in Fig. 6. Remarkable features of the solution are

• In the embedding matrix r > b the isotropic stress

1 3 kk

τ vanishes so that the incremental stress field is purely deviatoric: moreover, measurements taken in this region, outside the TPE source, cannot provide any evidence for the existence of a TPE region in a < r < b since the solutions are identical to those provided by a pressurized sphere, with overpressure

6

M.E. Belardinelli et al. / Earth and Planetary Science Letters 525 (2019) 115765

Fig. 5. (a) TPE region (yellow) enclosing a pressurized fluid-like reservoir (orange) and surrounded by an unbounded matrix (gray) in drained isothermal conditions. (b) Local reference system employed to compute the Coulomb failure function in the neighborhood of points vertically above the reservoir, like point S (in the TPE source) or point M (in the embedding matrix).

Fig. 6. Normalized displacement potential × and non-vanishing stress components



1 2π a2



, radial displacement u ×



3 a1



3 τrr , τθ θ , τϕϕ × [ 4μ ] due to a uniform in1

crease of temperature  T and pore pressure  p within a TPE shell (yellow) with inner radius a, outer radius b = 1.5 a, embedded within a drained isothermal matrix (gray). The isotropic stress 13 τkk and maximum shear stress τmax = 12 (τθ θ − τrr ) are also shown (normalized). All quantities are incremental values with respect to the initial condition before degassing and the overpressure within the magma chamber (orange domain) is assumed constant P 0 = P m during degassing. Arrows at the top describe the principal stress axes and the expected fault mechanisms, inverse within the TPE shell (as at point S in Fig. 5) and normal in the matrix (as at point M).

 P = 43 μ1 (b a−3a ) , surrounded by an elastic medium. In par3

Fig. 7. Solid lines show the normalized radial displacement u × vanishing stress components



3 a1

 and non-

3 τrr , τθ θ , τϕϕ × [ 4μ ] due to a uniform increase of 1

temperature  T and pore pressure  p within a spherical shell with b = 1.5 a, embedded within a matrix in isothermal drained conditions. All quantities are incremental values with respect to the initial condition before degassing and the overpressure within the magma chamber drops from P m to P 0 = P m − 0.2  p during degassing. Dashed lines refer to the case P m = P 0 (considered in Fig. 6) and are shown for ease of comparison. The isotropic stress 13 τkk (green line) is unchanged if

r > a but the magnitude of the maximum shear stress τmax = 12 (τθ θ − τrr ) (red) increases significantly within the TPE shell and decreases outside. Arrows at the top describe the expected fault mechanisms, inverse within the TPE shell (as at point S in Fig. 5) and normal in the matrix (as at point M).

3

ticular, if the ground surface is envisaged laying horizontally above the shell, normal faulting processes are expected to take place in the matrix domain around point M (in Fig. 5). • Within the TPE shell all incremental stress components are negative (compressive) as in the previous case of the spherical TPE source shown in Fig. 4; the isotropic stress 13 τkk is constant within the TPE shell but a significant deviatoric component is present, with opposite sign w.r. to the embedding matrix, promoting thrust faulting processes in the region around point S (in Fig. 5). • Employing parameter values μ = 6 GPa, ν = 0.2, H = 10 GPa, α = 3 × 10−5 K−1 pertinent to highly porous sedimentary rocks (e.g. Berea sandstone, Rice and Cleary, 1976), a temperature increase  T = 10 K and a pore pressure increase  p = 1 MPa in a spherical shell with a = 1 km and b = 1.5 km provide the following numerical values of the normalization constants: 43 μ1 = 1.6 × 106 Pa for stress components, 13 1 a = 20 cm for the displacement; the thermo-poro-elastic effects appear to be significant. Assuming the previous parameter values, the thermal contribution is 3 times the poro-elastic contribution, but the importance of  p in promoting faulting processes within the TPE region is enhanced by its lubrication effects, as discussed in the following section. • If degassing is accompanied by a significant decrease of pressure from P m to P 0 < P m within the magma reservoir, we get additional purely deviatoric stress, positive for the rr com-

 components which θ and ϕϕ ponents and negative for the θ increase the magnitude of the maximum shear stress in the TPE region and decrease it in the embedding matrix. This additional contribution is shown in Fig. 7 assuming the same parameter values considered above and P m − P 0 = 0.2  p = 0.2 MPa. 6.2. Coulomb failure function for the shell model Solutions obtained in the previous sections show that moderate changes of temperature (e.g.  T = 10 K) and pore pressure (e.g.  p = 1 MPa) lead to significant stress changes within and outside the TPE shell; it may be mentioned that the difference between lithostatic and hydrostatic pressure at 3 km depth is typically 50 MPa, and this may be considered as the upper limit for the pore pressure changes. While the increasing shear stress shifts faults toward instability, the increasing confining pressure tends to lock them frictionally; moreover, the increasing pore pressure tends to unlock faults due to the decreasing effective normal compression. In order to assess the cumulative effect of these stress changes on fault stability, we consider the top portion of the shell (Fig. 5-b), where τzz = τrr acts vertically and τθθ = τxx horizontally (the same considerations apply if τϕϕ = τ y y is considered); we compute the change of the Coulomb failure function C F F due to thermo-poro-elastic deformation (see e.g. Nostro et al., 1997).

M.E. Belardinelli et al. / Earth and Planetary Science Letters 525 (2019) 115765

7

 T / p (when β → ∞) C F F is positive only for low friction coefficients. Fig. 9 shows C F F vs r /a for a few values of the friction coefficient f ; in panel (a) β = 0 (i.e.  T = 0), in panel (b) β = 1, in panel (c) β = 3 and in panel (d) β = ∞ (i.e.  p = 0). It is to be noted that the role of friction is relevant and it is completely different according to whether  T is low or if  p is low: in the former case (panels a and b) higher friction is associated with higher C F F while the opposite happens when  p → 0 (panel c and d). Furthermore, the role of f is negligible in the embedding matrix (r > b) w.r. to the TPE domain a < r < b. 7. Discussion

Fig. 8. Incremental Coulomb failure function C F F (normalized to μ1 ) vs. the friction coefficient f , computed on optimally oriented planes within the TPE shell, with T outer radius b = 1.5 a, for a few values of β =   p α H in r = a (solid lines) and r = b (dashed). The TPE medium is characterized by and H = 10 GPa.

μ = 6 GPa, ν = 0.2, α = 3 × 10−5 K−1

C F F = |τ | + f (σn +  p ) where τ is the change of the shear traction acting on the fault plane, f is the coefficient of friction, σn (which is negative if compressive) is the change of the normal traction and  p is the pore pressure change. In the matrix domain, where τkk = 0 (as shown in the previous section for both the spherical and the spherical shell TPE sources), no change of pore pressure is induced by the deformation of the TPE shell, since  p = − B3 τkk = 0 is the initial condition for  p outside V S , according to Skempton (1954) law, where B denotes the Skempton parameter. Of course, τ and σn depend on the fault dip δ and performing a rotation by an angle δ around the y-axis (Fig. 5-b) we get:

1

τ = τz x = (τθ θ − τrr ) sin 2δ 2 ⎧ 3 ⎨ −μ1 ar 3 sin 2δ   = ⎩ μ b3 −a3 a3 sin 2δ 1

a3

1

r3

in a ≤ r ≤ b in r > b

1

σn = τz z = (τrr + τθ θ ) + (τrr − τθ θ ) cos 2δ 2 2 ⎧ 3 ⎨ − 43 μ1 + μ1 ar 3 ( 13 + cos 2δ) in a ≤ r ≤ b  3 3 3 = ⎩ −μ b −a a ( 1 + cos 2δ) in r > b 1 a3 r3 3 where the incremental stresses within the TPE shell are computed from eqs. (19). The maximum value of C F F in the TPE shell (for assigned  T and  p) is found on “optimally oriented faults” for which tan 2δ = 1/ f ; the maximum of C F F in the embedding matrix is found when tan 2δ = −1/ f . These values are the same as provided by the Anderson theory for thrust and normal faults, respectively. The value of C F F /μ1 on optimally oriented planes is plotted in Fig. 8 as a function of the friction coefficient f , for a few values of β , the ratio  T / p (suitably normalized multiplying by α H ), when r = a (inner radius) and r = b (outer radius). C F F is larger if the  T contribution to 1 is smaller; this is due to the fact that both  T and  p contribute to 1 but only  p adds (algebraically) to the normal traction σn according to the modified Coulomb criterion. It is interesting to note that C F F is positive for low  T / p values (e.g. β ≤ 3) regardless of the friction coefficient. For high

In spite of its simplicity, the spherically symmetric ThermoPoro-Elastic (TPE) model provides important clues regarding the order of magnitude of deformation and stress fields induced by temperature and pore-pressure variations in geothermal areas. In particular, the spherical shell model shows that a significant deviatoric stress is generated by positive increments of pore pressure  p and temperature  T within a TPE source (Figs. 6 and 7): if the ground surface is envisaged to lay above the shell, the compressive stress axis is horizontal within the TPE region, so that thrust faulting mechanisms are promoted there. On the contrary, in the embedding matrix (in drained isothermal conditions) above the TPE region, the stress field is entirely deviatoric, with vertical compression axis, exciting normal faulting mechanisms. Displacements and stresses are identical (apart from a multiplicative factor), to those provided by a pressurized spherical cavity in an unbounded elastic medium. Of course, the models presented above are applicable only to deep TPE sources since the effects due to the proximity of the free surface are not included. For the spherical source described in section 6 (Fig. 4), the presence of the free surface may be simply introduced employing the technique described in Bonafede (1990) and Bonafede and Ferrari (2009), using the method of images combined with a Galerkin vector providing vanishing surface tractions (e.g. Fung, 1965). The solution obtained in this way for the TPE spherical source is exact. The same procedure may be applied to obtain approximate solutions for the TPE shell located in a halfspace (the approximation consists in neglecting the stress induced by the image source and the Galerkin vector within the inner fluidlike reservoir, and it is acceptable if the depth of the center of the sphere is much greater (e.g. 3 times) than its radius). Omitting explicit calculations, the displacement field obtained in this way at the free surface z = 0, if the TPE shell is confined below the surface and the overpressure in the magma chamber does not change while degassing (i.e. if P 0 = P m ), is obtained as

u z (ρ ) = u ρ (ρ ) =

4(1 − ν )

1 (b3 − a3 )

3 4(1 − ν ) 3

z0

( z02

1 (b3 − a3 )

+ ρ 2 )3/2

ρ

( z02 + ρ 2 )3/2

where u z is the vertical displacement component (positive up√ ward), u ρ the horizontal component, ρ = x2 + z2 is the horizontal distance from the z axis and − z0 is the vertical position of the shell center below the free surface. These equations show that the displacement field at the free surface due to inflation of the TPE shell is identical to the displacement provided by a Mogi source (Mogi, 1958) which approximates a spherical pressurized cavity in (b3 −a3 )

a half-space with overpressure  P = 43 μ1 a3 (the same relation already obtained for the unbounded medium). Accordingly, data collected outside the TPE domain and even geodetic data collected at the surface can be fitted assuming a pressurized cavity

8

M.E. Belardinelli et al. / Earth and Planetary Science Letters 525 (2019) 115765

Fig. 9. Incremental Coulomb failure function C F F (normalized to μ1 ) vs r /a, computed on optimally oriented planes for a few values of f : (a) β = 0 (i.e.  T = 0), (b) β = 1, (c) β = 3 and (d) β = ∞ (i.e.  p = 0). Same parameters as in Fig. 8.

and cannot provide unquestionable evidence of the presence of a TPE source at depth. However we do get such indication if faulting processes with inverse mechanisms take place above the magma chamber during an inflation episode. Shear failure processes with inverse mechanisms are indicative of the presence of a TPE region and related fractures would lessen compression and provide preferential paths for magma intrusion processes across the TPE region. Since τθθ = τϕϕ , inverse ring faults (striking along ϕˆ ) and inverse radial

faults (striking along θˆ ) are equally stimulated. It must be noted that the previous considerations neglect the presence of an initial stress field, with significant deviatoric components, onto which the TPE stresses are superposed: for instance, Corbi et al. (2015) describe the stress state in a caldera considering topographic unloading following caldera formation, magma chamber decompression, magma buoyancy forces and tectonic stresses. At shallow depths these stress components are probably long-lasting, but viscoelastic relaxation plausibly attenuates the deviatoric components at great depth, since the effective viscosity decreases sharply at higher temperatures. Assuming that a quiescence period, longer than the Maxwell relaxation time, preceded the degassing event, the TPE stress contributions dominate the deviatoric stress field and the maximum regional shear stress may be smaller than the one induced by the TPE within V S . In this case, thrust faults are promoted within the TPE shell with strike direction parallel to the horizontal tensional tectonic axis. The present TPE shell model features reverse faulting mechanisms within the TPE shell at depth, and normal faulting at shallower depth, providing a vertical separation between the compressional and tensional stress domains. Other explanations may be proposed for the appearance of heterogeneous fault mechanisms at shallow depth in inflating calderas. For instance, it may be easily shown that normal faulting is predicted by Mogi (1958) model vertically above a pressurized spherical cavity, while strike-slip faulting and reverse faulting are promoted laterally; Bonafede and Danesi (1997) show that strike-slip faulting must also be expected laterally on both sides of a vertically dipping pressurized dike and

thrust faulting is predicted by Okada (1992) solutions close to a horizontal tensile dislocation surface, as inferred by D’Auria et al. (2015) at Campi Flegrei from a joint inversion of geodetic and seismic data, during the 1982-84 unrest. Moreover, according to analogue models (e.g. Figures 7-8 in Acocella, 2007), inward dipping reverse ring faults accompany laterally the formation of a resurgent dome during magma chamber inflation experiments. According to the present model, the intrusion of a magmatic dike across the TPE region, must be accompanied by the release of shear tractions over crack walls, giving rise to a mixed tensile and shear dislocation with a reverse dip-slip component. The energetically favorite dike propagation in this case is along an obliquely dipping plane, according to the numerical model of Maccaferri et al. (2011), corroborated by analog laboratory experiments (Menand et al., 2010; Menand, 2011): the presence of a compressive stress along the horizontal deflects a dike, even if it started vertically dipping from the magma chamber; in extreme cases, the numerical simulations predict that a dike may become arrested as a horizontal sill. As said previously, the present TPE shell model is not specifically designed to reproduce the Campi Flegrei unrest episodes but a mixed-mode dislocation event most probably generated the 1982-84 inflation episode (Trasatti et al., 2011). Leveling and EDM data are best fitted employing a moment tensor source compatible with an inclined dike opening in mixed mode (tensile and reverseshear dislocation) at ∼ 5 km depth, striking EW and dipping ∼ 30◦ to the North. The inversion procedure was performed by Trasatti et al. (2011) taking into account realistic rigidity heterogeneities in the Campi Flegrei caldera, as inferred from seismic tomography (Zollo et al., 2008) and the moment tensor mechanism reproduces a plausible value for the density of the intruded magma, indicating that uplift cannot be explained only by degassing processes but it was accompanied by significant magma intrusion. It is noteworthy that even the more recent 2011-2013 inflation episode at Campi Flegrei (Trasatti et al., 2015) could be fitted by a nearly identical mixed dislocation mechanism. The ground deformation field, retrieved from dense SAR images along ascending and

M.E. Belardinelli et al. / Earth and Planetary Science Letters 525 (2019) 115765

descending orbits clearly suggests that the previously emplaced dike, possibly initiated more than 50 years ago, according to Del Gaudio et al. (2010), further propagated with a reverse dip-slip component accompanying dike opening. A similar deformation mechanism, considering the superposition of magma chamber inflation and aseismic inverse faulting over a system of ring faults was proposed by De Natale et al. (1997) and Troise et al. (2003) for the 1982-84 episode employing a homogeneous half-space model. The present model assumes that uniform temperature and pore pressure changes are present within the TPE shell. It is important to remark that in most cases  T and  p typically extend over different spatial domains: if temperature changes take place following migration of hot pressurized fluids across a colder permeable medium, then pressure changes typically migrate faster than temperature changes, unless permeability is extremely low (see e.g. figure 6-7 in Bonafede and Mazzanti, 1997 and Zencher et al., 2006). According to Fig. 9-a, we expect that larger C F F values would be attained within the pressurized shell external to the heated shell and, according to Fig. 9-d, a region may become seismically silent once it is pervaded by significant temperature increase. The model also predicts that a negative  T (with the same  p) has de-stabilizing effects. This mechanism has thus intriguing implications for earthquakes triggered by fluid injections in geothermal reservoirs (see e.g. Shapiro, 2015). Of course, the geometric shape of the TPE domain might be much more complex than assumed above: the spatial distribution of pore pressure is extremely sensitive to the permeability pattern and to the presence of fracture zones for which we usually have little - if any - “a priori” knowledge. In such cases, the method described by eqs. (11)-(13) cannot be applied since the Green’s function is no longer symmetric and resort must be made directly to equations (7)-(8). Furthermore, the previous models assume  p and  T values which are constant in time and space throughout the TPE region while they should more realistically propagate and vanish gradually when approaching the TPE boundary, where we should assign appropriate boundary conditions for pore pressure and temperature. Such generalizations are left for future developments, but the present simple models clearly show the importance of thermo-poro-elastic processes in producing rock deformation and triggering seismic activity with contrasting mechanisms in a restless caldera like Campi Flegrei. Finally, the TPE model may also apply (with appropriate changes of fluid parameter values) to describe transcrustal magmatic systems (TCMS) dominated by crystal mush (a system of crystals and melt in which the crystals form a continuous framework through which melt is distributed; e.g. Cashman et al., 2017). Acknowledgements Elisa Trasatti and an anonymous reviewer are acknowledged for constructive suggestions on the previous version of the present paper. Elisa Trasatti wishes to join the authors in the grateful memory of Prof. Enzo Boschi to whom this paper is dedicated. References Acocella, V., 2007. Understanding caldera structure and development: an overview of analogue models compared to natural calderas. Earth-Sci. Rev. 85, 125–160. Aki, K., Richards, P., 1980. Quantitative Seismology: Theory and Methods. A Series of Books in Geology, vol. 1. W.H. Freeman. Bonafede, M., 1990. Axi-symmetric deformation of a thermo-poro-elastic half-space: inflation of a magma chamber. Geophys. J. Int. 103 (2), 289–299.

9

Bonafede, M., 1991. Hot fluid migration: an efficient source of ground deformation: application to the 1982-1985 crisis at Campi Flegrei-Italy. J. Volcanol. Geotherm. Res. 48 (1), 187–198. Bonafede, M., Danesi, S., 1997. Near-field modifications of stress induced by dyke injection at shallow depth. Geophys. J. Int. 130, 435–448. Bonafede, M., Ferrari, C., 2009. Analytical models of deformation and residual gravity changes due to a Mogi source in a viscoelastic medium. Tectonophysics 471 (1), 4–13. Understanding stress and deformation in active volcanoes. Bonafede, M., Mazzanti, M., 1997. Hot fluid migration in compressible saturated porous media. Geophys. J. Int. 128, 383–398. Cashman, K.V., Sparks, R.S.J., Blundy, J.D., 2017. Vertically extensive and unstable magmatic systems: a unified view of igneous processes. Science 355, eaag3055. Chiodini, G., Todesco, M., Caliro, S., Del Gaudio, C., Macedonio, G., Russo, M., 2003. Magma degassing as a trigger of bradyseismic events: the case of Phlegrean Fields (Italy). Geophys. Res. Lett. 30 (8). Chiodini, G., Vandemeulebrouck, J., Caliro, S., D’Auria, L., Martino, P.D., Mangiacapra, A., Petrillo, Z., 2015. Evidence of thermal-driven processes triggering the 2005–2014 unrest at Campi Flegrei caldera. Earth Planet. Sci. Lett. 414, 58–67. Corbi, F., Rivalta, E., Pinel, V., Maccaferri, F., Bagnardic, M., Acocella, V., 2015. How caldera collapse shapes the shallow emplacement and transfer of magma in active volcanoes. Earth Planet. Sci. Lett. 2015 (431), 287–293. D’Auria, L., Giudicepietro, F., Aquino, I., Borriello, G., Del Gaudio, C., Lo Bascio, D., Martini, M., Ricciardi, G.P., Riciolino, P., Ricco, C., 2011. Repeated fluid-transfer episodes as a mechanism for the recent dynamics of Campi Flegrei caldera (1989-2010). J. Geophys. Res., Solid Earth 116, B04313. D’Auria, L., Massa, B., Cristiano, E., Del Gaudio, C., Giudicepietro, F., Ricciardi, G., Ricco, C., 2015. Retrieving the stress field within the Campi Flegrei caldera (southern Italy) through an integrated geodetical and seismological approach. Pure Appl. Geophys. 172 (11), 3247–3263. De Natale, G., Petrazzuoli, S.M., Pingue, F., 1997. The effect of collapse structures on ground deformations in calderas. Geophys. Res. Lett. 24 (13), 1555–1558. Del Gaudio, C., Aquino, I., Ricciardi, G.P., Ricco, C., Scandone, R., 2010. Unrest episodes at Campi Flegrei: a reconstruction of vertical ground movements during 1905-2009. J. Volcanol. Geotherm. Res. 195 (1), 48–56. Dragoni, M., Magnanensi, C., 1989. Displacement and stress produced by a pressurized, spherical magma chamber, surrounded by a viscoelastic shell. Phys. Earth Planet. Inter. 56 (3), 316–328. Eshelby, J.D., 1961. Elastic Inclusions and Inhomogeneities. Progress in Solid Mechanics, vol. 2. North-Holland Publishing Company, Amsterdam. Fung, Y., 1965. Foundations of Solid Mechanics. Hall, Englewood Cliffs, NJ. Landau, L., Lifschitz, E., 1967. Theorie de l’Elasticité. Editions Mir, Moscow. Lima, A., Vivo, B.D., Spera, F.J., Bodnar, R.J., Milia, A., Nunziata, C., Belkin, H.E., Cannatelli, C., 2009. Thermodynamic model for uplift and deflation episodes (bradyseism) associated with magmatic–hydrothermal activity at the Campi Flegrei (Italy). Earth-Sci. Rev. 97 (1), 44–58. Maccaferri, F., Bonafede, M., Rivalta, E., 2011. A quantitative study of the mechanisms governing dike propagation, dike arrest and sill formation. J. Volcanol. Geotherm. Res. 208 (1), 39–50. McTigue, D.F., 1986. Thermoelastic response of fluid-saturated porous rock. J. Geophys. Res., Solid Earth 91 (B9), 9533–9542. Menand, T., 2011. Physical controls and depth of emplacement of igneous bodies: a review. Tectonophysics 500 (1), 11–19. Emplacement of magma pulses and growth of magma bodies. Menand, T., Daniels, K.A., Benghiat, P., 2010. Dyke propagation and sill formation in a compressive tectonic environment. J. Geophys. Res., Solid Earth 115 (B8). Mogi, K., 1958. Relation between the eruptions of various volcanoes and deformations of the ground surfaces around them. Bull. Earthq. Res. Inst. 36, 99–134. Nostro, C., Cocco, M., Belardinelli, M.E., 1997. Static stress changes in extensional regimes: an application to southern Apennines (Italy). Bull. Seismol. Soc. Am. 87 (1), 234–248. Okada, Y., 1992. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 82 (2), 1018–1040. Rice, J.R., Cleary, M.P., 1976. Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents. Rev. Geophys. 14 (2), 227–241. Rinaldi, A., Todesco, M., Bonafede, M., 2010. Hydrothermal instability and ground displacement at the Campi Flegrei caldera. Phys. Earth Planet. Inter. 178 (3), 155–161. Segall, P., 2016. Repressurization following eruption from a magma chamber with a viscoelastic aureole. J. Geophys. Res., Solid Earth 121, 8501–8522. Shapiro, S.A., 2015. Fluid-Induced Seismicity. Cambridge University Press. Skempton, A.W., 1954. The pore-pressure coefficients A and B. Geotechnique 4 (4), 143–147. Todesco, M., Berrino, G., 2005. Modeling hydrothermal fluid circulation and gravity signals at the Phlegraean Fields caldera. Earth Planet. Sci. Lett. 240 (2), 328–338. Trasatti, E., Bonafede, M., Ferrari, C., Giunchi, C., Berrino, G., 2011. On deformation sources in volcanic areas: modeling the Campi Flegrei (Italy) 1982-84 unrest. Earth Planet. Sci. Lett. 306 (3), 175–185.

10

M.E. Belardinelli et al. / Earth and Planetary Science Letters 525 (2019) 115765

Trasatti, E., Polcari, M., Bonafede, M., Stramondo, S., 2015. Geodetic constraints to the source mechanism of the 2011-2013 unrest at Campi Flegrei (Italy) caldera. Geophys. Res. Lett. 42 (10), 3847–3854. Troise, C., Pingue, F., De Natale, G., 2003. Coulomb stress changes at calderas: modeling the seismicity of Campi Flegrei (southern Italy). J. Geophys. Res., Solid Earth 108 (B6).

Zencher, F., Bonafede, M., Stefansson, R., 2006. Near-lithostatic pore pressure at seismogenic depths: a thermoporoelastic model. Geophys. J. Int. 166 (3), 1318–1334. Zollo, A., Maercklin, N., Vassallo, M., Dello Iacono, D., Virieux, J., Gasparini, P., 2008. Seismic reflections reveal a massive melt layer feeding Campi Flegrei caldera. Geophys. Res. Lett. 35 (12).