Elastoplastic source model for microseismicity and acoustic emission

isotropic rocks and ‘wet’ isotropic rocks. The model highlights theoretical links between stress state, geomechanical parameters and conventional representations of the moment tensor Hudson source type


Background and definitions
Microseismicity and acoustic emission research is a part of earthquake science which focuses on weak seismic signals that often lack a clear main shock and form so-called microseismic clouds and earthquake swarms. In comparison to normal earthquakes, microseismic events are characterized by smaller magnitudes, higher frequencies, shorter duration and a more complex source mechanism (Foulger et al. 2004;Kamei et al. 2015;Eaton 2018). Microseismicity is observed both in natural conditions and during geotechnical operations (Vavryčuk 2002;Fischer & Guest 2011). The induced activity can be associated with hydraulic fracturing (Baig & Urbancic 2010;Ellsworth 2013;Holland 2013;Schultz et al. 2020), stimulation of geothermal reservoirs (Majer et al. 2007;Deichmann & Giardini 2009;Dorbath et al. 2009;Zang et al. 2014), waste water injection (Keranen et al. 2014), mining (Manthei & Eisenblätter 2008;Sen et al. 2013), tunnel excavations and borehole breakouts (Paul Young & Martin 1993;Hazzard & Young 2002;Read 2004).
The precise term definition is not well agreed and depends on application (Hardy 2003;Grosse & Ohtsu 2008;Eaton 2018). In this work, we reserve the word 'microseismicity' for small events (M < 0) that occur in geological environment (either naturally or induced by human operations) with a frequency of 10 1 −10 3 Hz, whereas mainly ultrasound 'acoustic emissions' (10 4 −10 6 Hz) is referred to laboratory deformation experiments on rock specimens.
The acoustic emissions are elastic waves generated by a sudden release of mechanical energy by means of micro/grain-scale deformation processes, for example microfracturing inside and between mineral grains, grain crushing and boundary sliding, deformation twinning and phase transitions (Zang et al. 1996;Fortin et al. 2009). The acoustic emissions are considered as a precursor of large-scale faulting, and can be used to obtain continuous data at various stages of the deformation process: starting from the distributed plastic failure towards a localized macroscopic shear (Lockner 1993;Amitrano 2006;Manthei & Eisenblätter 2008;Dresen et al. 2010;Renard et al. 2017;Lei 2017;Wei et al. 2019). The spatial distribution of microseismic/acoustic emission events indicates the location of cracks and can be used to delineate the failure region. The source mechanism provides information on the mode of failure: tensile crack, shear crack or their combination (Enoki & Kishi 1988;Grosse & Ohtsu 2008;Lei 2017;Eaton 2018).
In this paper, we suggest a simple mathematical model that helps to compute the seismic moment tensor in numerical models of rock failure based on continuum mechanics approach. This model can be used for interpretation of laboratory acoustic emission data as well as to link microseismic signals to the progressive rock deformation associated with hydraulic fracturing, mine excavation and other geotechnical operations.

Observations
Microseismic monitoring data often indicate distributed clouds of seismic events around a borehole together with a more aligned seismicity. In applications related to fluid injection and hydraulic fracturing operations, microseismic events around a wellbore are often explained by activation of pre-existing fractures facilitated by pore pressure diffusion into the surrounding rock (Zoback 2010;Fischer & Guest 2011;Shapiro 2015). The shear events are present alongside with tensile events expected for fluid-saturated rocks (Šílený et al. 2009;Fischer & Guest 2011;Eaton et al. 2014b;Zhang et al. 2019).
Significant isotropic components are common for microseismic sources (Foulger et al. 2004;Baig & Urbancic 2010;Fischer & Guest 2011;Sen et al. 2013;Zhang et al. 2016). The observation of non-double couple (non-DC) components in the source mechanism of both natural micro-earthquakes and operational microseismicity motivated to revise models of the seismic source mechanism. The suggested models include heterogeneity and elastic anisotropy (Vavryčuk 2005), tensile failure facilitated by fluid injection (Šílený et al. 2009;Fischer & Guest 2011), co-seismic brittle damage, complex faulting and dilatant jogs (Julian et al. 1998;Zhang et al. 2016) and other complexities of the source compared to idealized slip on a fault plane in homogeneous medium (Frohlich 1994;Julian et al. 1998).
Laboratory experiments are important to better understand how plastic failure mechanisms are linked to the percentage of isotropic and deviatoric components of the seismic moment tensor. Aker et al. (2014) investigated the acoustic emissions associated with developing faults originating from a horizontal borehole in sandstone sample under triaxial compression. Vera-Rodriguez et al. (2017) performed a hydraulic fracturing experiment in sandstone under triaxial compression and recorded the acoustic emission. The full moment tensor solutions were obtained in both studies. The majority of events of Aker et al. (2014) and the events corresponding to the dry fracture period of Vera-Rodriguez et al. (2017) are situated in the second quadrant (upper-left part) of the Hudson plot (Hudson et al. 1989) which corresponds to a combination of DC, compensated linear vector dipole (CLVD) and volumetric extension (+V) components. A small number of events occupy the forth quadrant (right-lower part) corresponding to DC, CLVD and volumetric compression (-V). Aker et al. (2014) made several other interesting observations regarding the isotropic component of moment tensor: (i) events with large moment magnitude have a small volumetric component, (ii) the percentage of DC component increases with progressive loading and (iii) events close to the borehole have a larger isotropic component, compared to the events farther away from the borehole. In this paper, we will try to explain these and some other microseismic observations based on classical continuum mechanics.

Approaches to mechanical modeling of seismic events
There is no standard methodology for modeling the mechanical processes causing seismic events. In computational earthquake seismology, the dynamic rupture simulation (DRS) is typically used to predict the propagation of rupture on a fault plane and associated seismic motions in the entire model (Andrews 1999;Dalguer & Day 2006;Aagaard et al. 2013;Galvez et al. 2014;Lyakhovsky et al. 2016;Harris et al. 2018). The earthquake cycle can be modelled using time-and rate-dependent friction laws on the fault (Kanamori & Brodsky 2004;Scholz 2019). Effects of fluid and off-fault plasticity can also be included in DRS Viesca et al. 2008;Roten et al. 2017;Jin & Zoback 2018;Okubo et al. 2019). The DRS method has mainly been applied in simulations of earthquakes with a DC source mechanism.
Another class of models represents earth materials as an assembly of discrete particles that might be connected by beams or bonds of other nature (see Lisjak et al. 2013;Wangen 2018). The particle flow code (PFC, Itasca 1999;Potyondy & Cundall 2004) or hybrid methods including PFC (Itasca 2000;Cai et al. 2007) seem to be particularly popular in simulation of fracture and acoustic emission in geotechnical applications. The constitutive relations are typically reproduced by calibrating microstructural parameters. The PFC modeling results have been applied in the operational design of various structures such as mines, boreholes and tunnels (Hazzard & Young 2002;Cai et al. 2007). A microseismic moment tensor response due to a dry tunnel failure was modelled by Hazzard & Young (2002). This method was later extended to fluid-saturated rocks by Zhao & Paul Young (2011). Raziperchikolaee et al. (2014) developed a similar microstructural modeling approach to obtain source mechanisms of acoustic emission during deformation at rock joints including effects of fluid flow. The PFC method was shown to reproduce various types of inelastic behaviour of rocks and associated complex source mechanisms of acoustic emission (Hazzard & Young 2002;Potyondy & Cundall 2004;Raziperchikolaee et al. 2014). However, the numerical micromechanical parameters (such as particle size distribution, bond strength and contact friction) are difficult to relate to physical properties. Downloaded from https://academic.oup.com/gji/article/227/1/33/6289921 by guest on 08 July 2021 In contrast, our model is built on classical continuum mechanics approach which uses standard rock mechanical properties. The model allows to derive analytical functional relations between geomechanical and seismic parameters. We follow the non-associative plasticity theory that helps i) to explain some of the acoustic emission observations made in rock testing and ii) to apply the experimental constitutive laws of rock failure to study microseismicty in heterogeneous geological environments. The full moment tensor source, associated with dilatant plastic shear bands and zones of localized tensile failure, can be explicitly calculated for each grid element in our numerical model. We do not include slip on a fault since we specifically focus on microseismicity and acoustic emission that can be considered as a precursor to macroscale shear failure.

Outline of this study
In the following, we modify the standard 'stress glut' formalism for the seismic moment tensor sources (Aki & Richards 1980;Ben-Menahem & Singh 1981;Vavryčuk 2005;Ampuero & Dahlen 2005;Vavryčuk 2011;Madariaga 2015). Specifically, the moment density tensor describing a general (non-DC) source mechanism can be represented using an inelastic strain tensor, for example (Madariaga 2015). Here, we extend the existing theory by including a stress-and stress-path-dependent rock behaviour during inelastic deformation. We show how to calculate the plastic strain tensor and provide analytic expressions of the moment density tensor for plastic failure in 2-D plane strain conditions. This simple physical model allows for better understanding the first-order relations between seismic observations and failure process associated with microseismic events. The model describes the shear and volumetric source mechanisms as a result of developing plastic strain primarily in the form of dilatant shear bands and tensile failure. We predict quantitatively the radiation pattern of seismic waves due to local changes of the stress tensor induced by spontaneous strain localization. To illustrate our concept, we perform numerical simulations assuming the 'dry' and 'wet' rock scenaria in 2-D plane strain loading conditions. As example, we calculate the earthquake source mechanisms induced by a failure around a cylindrical opening in the rock mass (e.g. boreholes or tunnels), the process that was extensively studied both in nature and in the lab. The accompanying datasets and MATLAB scripts to illustrate the model are available in zenodo.org, under doi: 10.5281/zenodo.4911325.

THEORY
In Section 2.1, we present equations of motion describing seismic wave propagation. Following previous works, we introduce inelastic strain tensor to describe the seismic source. After that (Section 2.2), we suggest an interpretation of inelastic deformation in term of instantaneous plastic strain tensor. In Section 2.3, we show that the plastic strain tensor can be obtained from derivatives of the plastic potential in the stress space. The plasticity constitutive relations are highly nonlinear, and, except for simple cases, require a numerical solution. In Section 2.4, we discuss the standard approaches for the mathematical description of the plastic failure criteria and the plastic potential. In Section 2.5, we extend the model to the case of poro-elastoplastic deformation including the effects of fluid pressure diffusion and seepage forces.

Equations of motion
Our approach is based on classical continuum mechanics and generally follows (Aki & Richards 1980;Ben-Zion & Ampuero 2009;Madariaga 2015). The standard elastic wave equations can be written for small stress increments σ . The source term should account for inelastic deformation responsible for initiation of seismic waves. The linearized equations of motion can be written as: (1) where the incremental stress tensor σ and displacements u vary with time t and position x. The equivalent body force f is related to the seismic moment density tensor: We use continuum mechanics sign convention and consider tensile stresses and strains as positive, such as a pure implosion source corresponds to a negative unit moment density tensor (Dahlen & Tromp 1998). The seismic moment tensor can be defined as a 'stress glut' (Backus & Mulcahy 1976) or stress produced by inelastic deformation of a body that is elastic everywhere outside the source region (Madariaga 2015): where D is the elastic stiffness tensor and ε T is inelastic or transformational strain which develops in the earthquakes source region (Eshelby 1957;Kostrov 1974;Backus & Mulcahy 1976;Ben-Zion & Ampuero 2009). The moment tensor representation of seismic sources is useful for efficient solution of forward and inverse problems of seismic wave propagation (Moczo et al. 2014) and allows for better understanding of failure mechanisms associated with hydraulic fracturing and borehole breakout processes (Baig & Urbancic 2010;Vera-Rodriguez et al. 2017). Yet, physical interpretation of the seismic moment tensor remains somewhat ambiguous. We aim to clarify this issue for the particular case of incipient faulting induced by stress perturbation.

Moment tensor source due to material failure
The inelastic deformation in the crust is a strongly nonlinear function of stress tensor components. A number of inelastic processes may contribute to the acoustic emission including fracturing, formation of shear bands, distributed plastic flow and creep. The fracturing and plastic flow are associated with instantaneous development of permanent rate-independent deformation, when stresses reach a threshold value, and can be modeled using plasticity theory (Kachanov 1971;Rice & Rudnicki 1980). The viscous flow due to creep is a slow, time-dependent accumulation of permanent deformation that does not necessarily require a threshold stress.
The seismic waveform analysis of natural and induced earthquakes suggests that the dynamic moment tensor source (3) can often be represented as a product of static moment tensor and a source time function (Aki & Richards 1980;Vavryčuk & Kühn 2012): The static partM(x) controls the seismic energy radiation pattern, whereas the time-dependent behaviour M(t) accounts for the history of stress release. Here, we focus on the static seismic moment tensor and do not directly address the source dynamics and wave propagation.
We propose the instantaneous rate-independent plastic failure as a process responsible forM(x). The inelastic strains (ε T ) in eq (3) can be approximated by the plastic strains at time t 0 : With this assumption, the concept of 'strain glut' (Ampuero & Dahlen 2005) acquires a meaning of the plastic strain tensor since ε p controls the pattern of seismic waves radiation. We do not consider other processes, such as creep and material damage that may also affect the radiation pattern.
The moment tensor is the volume integral of the moment density and can be expressed using eq. (3) For isotropic rocks, the moment tensor density becomeŝ where I is identity matrix, tr indicates the sum of the diagonal components and λ, μ are the Lamé parameters.

Quasi-static elastoplastic problem
In order to obtain the static moment tensor (6), we have to know the plastic strain tensor ε p . In the quasi-static approximation, implied by eq. (5), the elastoplastic problem can be solved based on the static form of the equations of motion (1), the incremental form of the generalized Hooke's law: the decomposition of the total strain into elastic and plastic components (Kachanov 1971): and the elastoplastic rheology described by following equations (Davis & Selvadurai 2005;Yarushina et al. 2010): where is a positive scalar plastic multiplier, F is the yield function and Q is the plastic flow potential. Closed expressions for the yield function F = F(σ ) and the plastic flow potential Q = Q(σ ) are known, so that the partial derivatives ∂ σ F and ∂ σ Q can be calculated analytically. Some examples of such functions will be given in the next sections of the paper. Eq. (11) is the condition for plastic loading which ensures stresses to remain at the yield surface. The solution of the elastoplastic problem requires boundary conditions that can be supplied as a stress at the external boundary of the model. The total plastic strain in the moment tensor equation (6) can be obtained as a sum of N incremental plastic strains (defined by eq. 10) over time or loading path: The incremental form of the Hooke's law (8) can be rearranged using eqs (9), (10) and (11) to find the elastoplastic rheology relation between the increments of the total strain and the stress increments: where Downloaded from https://academic.oup.com/gji/article/227/1/33/6289921 by guest on 08 July 2021 Note that there is no general relation between stresses and deformation in the plastic regime. The stress increments and strain increments are connected through the stress-dependent elastoplastic material matrix D ep . The strain increments are related to displacements as (Davis & Selvadurai 2005;Yu 2007

Yield surface and flow potential
The condition for plastic yielding is generally defined as a function of stress: In principal stress coordinates, this equation defines a surface such as shown in Fig. 1. The outward normal to this surface is ∇ F = ∂ F ∂σ . In the elastic regime, the stresses are inside this surface [F(σ ) < 0]. During progressive plastic deformation the stresses stay on the yield surface [F(σ ) = 0 and ∇F = 0].
The plastic strain rate tensor can be related to stresses via experimentally proven plastic flow rules. The plastic flow is 'associated' if both the principal stresses and strain rates (written as vectors) are normal to the yield surface (Davis & Selvadurai 2005). In this case, the following rule is applied to define the plastic strain tensor: In plasticity theory (Kachanov 1971;Yu 2007), a potential Q was introduced to account for some observations that the associated plasticity theory could not predict: (i) smaller plastic dilation; (ii) deviation of kinematic and stress characteristics and (iii) phenomenon of strain localization in shear bands (Rice & Rudnicki 1980). In the case of 'non-associated' flow, described by eq. (10), the vector of principal values of strain tensor is orthogonal to the plastic flow surface Q(σ ) = 0 which can deviate from the yield surface F(σ ) = 0. It is essential to note that plastic yielding occurs only when stresses reach the yield surface and remain on it.
Different types of failure are described using different failure criteria: tensile (Griffith), shear (Mohr-Coulomb) and compaction (elliptical cap) models. In the absence of hardening (also called 'perfect plasticity'), the yield function depends only on stress components. Here, the yield function follows the standard Mohr-Coulomb relation for shear failure (Jaeger et al. 2009): where φ is the friction angle, c is the cohesion or the shear strength at zero normal stress. The normal stress is defined as negative in compression This sign convention also implies that the dilatation is positive. The plastic flow potential takes the form (Davis & Selvadurai 2005) where σ 1 and σ 3 are maximum and minimum principal stresses (σ 1 > σ 2 > σ 3 ), and ψ is dilation angle. In general non-associative flow, the friction and dilation angle are different, and ψ < φ. The tensile failure criterion and plastic flow potential can be taken in the following simplest form (Rozhko et al. 2007) The failure mode in the plastic state is controlled by the failure criterion with the largest value of the yield function [F(σ ) ≤ 0] In the principal stress axes, this combined yield surface represents a truncated hexagonal cone (Fig. 1a). The plastic flow defined either by eq. (10) or eq. (18) is always dilational whenever the shape of the yield surface expands with increasing mean stress. This conclusion can be drawn from the energy considerations (see Davis & Selvadurai 2005).

Effect of pore fluid pressure
Previous sections highlight the relation between plastic strain tensor and seismic moment tensor. Based on laboratory and field observations, it has long been recognized that fluid flow is important factor controlling nucleation and growth of fractures and associated seismicity (Sibson 1981;Masuda et al. 1990;Miller et al. 1996;Li et al. 2016). Therefore, the effects of pore fluid pressure should be considered for the interpretation of microseismicity data. For a general description of seismic sources in fluid-saturated porous rocks, one has to solve a coupled fluid flow and solid deformation problem (Biot 1941;Rice & Cleary 1976;Rudnicki 1985;Yarushina & Podladchikov 2015). A basic demonstration of fluid pressure effects on the failure process is also possible following a less rigorous 'decoupled' modelling approach, for example (Rozhko et al. 2007). Despite the simplifying assumptions, the decoupled model captures two important characteristics of porous rock deformation: (1) the effective stress controlling deformation according to Terzaghi's principle and (2) the seepage force arising due to the gradient of fluid pressure. In the following, we will modify our constitutive and momentum conservation equations to account for these two effects.
The macroscopic incremental stresses in the isotropic porous elastic medium with fluid pressure perturbation p f can be related to macroscopic incremental strains as (Biot 1941;Coussy 2005;Yarushina & Podladchikov 2015): where the Biot-Willis coefficient α (0 < α < 1) describes the efficiency of pore pressure counteracting the effect of confining pressure in changing bulk volume. Force balance eq. (1) can now be rewritten using effective stresses and the poroelastic stress-strain relation. In the static limit, it reduces to Assuming that α is constant, this equation is similar to the standard static equilibrium equation that we used for dry rocks but also includes a force term in the form of the fluid pressure gradient (so-called 'seepage force'; Rozhko (2010)). In practice, the fluid flow can often be modelled separately from solid deformation (e.g. Rozhko 2010;Shapiro 2015). In this case, the pore pressure can be obtained by solving diffusion equatioṅ with specified initial and boundary conditions. The hydraulic diffusivity D d controls the position of the fluid diffusion front l d = √ 4π D d t. A number of analytical solutions are available for simple geometries. For example, the pore pressure around a pressurized spherical cavity, as an analog to fluid injection through an open short borehole section, has the form (Shapiro 2015): where p 0 is the fluid pressure at the surface of the cavity (R 0 ) and r is the radial distance. The failure of poro-elastoplastic media is controlled by Terzaghi's effective stress σ + p f I (Rice 1975;Coussy 2005). Therefore, the failure criteria and plastic flow potentials in eqs (21), (23) and (19) are modified to Downloaded from https://academic.oup.com/gji/article/227/1/33/6289921 by guest on 08 July 2021 for shear failure, and for tensile failure (Rozhko et al. 2007;Jaeger et al. 2009).

Moment tensor
The link of moment tensor components and the failure process can be illustrated based on a 2-D example in plane strain deformation around a cylindrical opening. This setup describes failure around tunnels and boreholes which was a focus of a number of laboratory and field studies. Eq. (10) written in terms of incremental principal strains ( ε 1 > ε 2 > ε 3 ) and principal stresses (σ 1 > σ 2 > σ 3 ) in plane strain are The principal plastic strain increments can also be expressed using the Mohr-Coulomb plastic flow potential in eq. (21) The volumetric incremental plastic strain becomes This indicates that the bulk plastic strain is always dilational since the plastic multiplier is positive ( > 0). The moment density tensorM described by eqs (6) and (10) depends on both plastic strains and elastic stiffness tensor. For isotropic solid, the principal moments can be written using the elastic Láme parameters λ and μ, the plastic dilation angle ψ, and the plastic multiplier : where ν = λ 2(λ+μ) is the Poisson ratio. In the practically interesting case of a transversely isotropic elastic material, the moment tensor density becomeŝ where D ij are the components of the elastic stiffness matrix. In both cases, the sum of diagonal elements of the moment tensor can be non-zero due to plastic dilation that would give rise to isotropic source component.

Radiation pattern of seismic waves
For a general failure caused by a combination of shear and volumetric deformation in x−y plane, the corresponding plastic strain tensor is If rocks are isotropic, we can use equation (7) to obtain the moment densitŷ The far-field radiation pattern of the outgoing seismic waves on the focal sphere surrounding the source is proportional to the far-field amplitude of the outgoing waves (Aki & Richards 1980). Following Dahlen & Tromp (1998), the radiation pattern can be expressed using normalized wave vector p and polarization vector η where symbol ':' denotes the double contraction operation. For the particular case of P waves, the wave propagation direction coincides with the polarization vector . Moment tensor can be obtained from the moment density by the surface integration: where A is the surface element of the plastic region, dl is an out-of-plane (vertical) length scale required to preserve the dimension of the moment (Nm).

Hudson source type parameters
A useful representation of the moment tensor is given by Hudson parameters κ and τ (Hudson et al. 1989). These parameters characterize the isotropic (κ) and shear components (τ ) of the source. For a general moment tensor, the principal moments can be derived using eigendecomposition. The Hudson parameters written for principal moments (eigenvalues) M 1 > M 2 > M 3 of the moment tensor M are where Volume expansion at the source region corresponds to κ > 0 whereas κ < 0 indicates implosion. A CLVD-source with a tensile largest principal moment corresponds to τ < 0. A DC-source is obtained if both parameters are zero (κ = 0, τ = 0).

Model setup
We illustrate the concepts developed in chapter 3 based on 2D numerical examples of plastic failure induced by external loading of a circular opening. This model geometry is often chosen in the analysis of borehole stability or tunnel excavation (e.g. Jaeger et al. 2009). The standard approach to estimate the minimum critical stress before rock failure is based on the elastic solution by Kirsch (1898). Further evolution of failure region and predictions of associated acoustic emission or miscroseismicity must be based on elastoplastic numerical solutions, that take into account development of the failure region.
We consider a vertical cylindrical opening in heterogeneous volume assuming a random distribution of material parameters (such as cohesion or Poisson ratio). The simulation of a random field was performed using the spectral randomization method (Räss et al. 2019). We use an in-house elastoplastic finite-element code for numerical simulations in which both shear and tensile failure are implemented (Yarushina et al. 2010;Minakov et al. 2018). The efficiency of the method and benchmarks of the code versus analytical solutions were presented in Yarushina et al. (2010). The numerical grid is structured but non-uniform, and it has a finer resolution close to the borehole. The time step is adaptive, depending on maximum incremental strain values. The domain is circular (to avoid stress concentrations at the model edges) with the outer radius of L and the inner radius of the circular opening R 0 (Fig. 2a). We perform numerical experiments for four special Downloaded from https://academic.oup.com/gji/article/227/1/33/6289921 by guest on 08 July 2021 (a) (b) Figure 2. Finite-element model geometry (a) and the random cohesion field (b) used for model initialization. Stresses (σ ∞ x , σ ∞ y ) are applied at the external radius L. The inner boundary (R 0 ) is free (dry case) or subjected to fluid pressure p 0 (wet case). Model parameters are specified in Table 1. cases assuming different rheology and boundary conditions: (1) a combination of hydrostatic compression and pure shear ( σ ∞ x > σ ∞ y ) in an isotropic dry rock, (2) a combination of hydrostatic compression and pure shear (σ ∞ x > σ ∞ y ) in an isotropic elastically heterogeneous dry rock, (3) hydrostatic compression in an anisotropic dry rock, where external stresses are σ ∞ x = σ ∞ y and (4) non-hydrostatic compression of an isotropic fluid-saturated porous rock. In all simulations, we use a random cohesion field with an exponential correlation function (Fig. 2b). The solution strategy follows an incremental first-order forward Euler scheme with a correction of the drift from the yield surface. The details of numerical implementation of the finite-element model can be found in Yarushina et al. (2010) and Minakov et al. (2018). The seismic source is defined according to eq. (45), such as acoustic waves are assumed to be emitted from within every cell of the model with non-zero plastic strain. Material and other numerical model parameters are listed in Table 1.

Isotropic dry rocks
To reproduce typical pattern of borehole failure under tectonic stress conditions, we apply a combined pressure and shear loading at the external boundary. The simulation results in Fig. 3 correspond to the pressure (mean stress) load of ≈60 MPa and shear load of ≈20 MPa (see also Table 1). The components of the elastic stiffness matrix in plane strain (Zienkiewicz & Taylor 2005) depend only on the Poisson ratio ν. The Young modulus E = 2μ(1 + ν) scales the magnitude of the components, and the stiffness matrix becomes Elastic and plastic strains (see eq. 9), caused by the applied load, are shown in Figs 3(a), (c) and (b), (d), respectively. The region of maximum elastic shear strain coincides with the maximum of plastic shear strain. The bulk elastic deformation is compressional whereas the plastic deformation exhibits dilatancy across shear bands. The loading of isotropic rocks beyond the failure limit leads to formation of two failure regions on the opposite sides of the borehole. The borehole breakouts occur in the direction orthogonal to the maximum compression, as expected. The plastic failure is due to development of localized shear bands arising spontaneously at small heterogeneities around the borehole.
Development of plastic deformation in the form of shear bands in Fig. 3 is linked to radiation of seismic waves in our model. We use eq. (6) with the stiffness matrix eq. (52) and the plastic strain tensor, obtained using finite element simulations, to compute the seismic moment tensor density on each element of the numerical grid. The scalar moment M 0 characterizes the intensity of seismic events: where the vertical moment M z = ν(M x + M y ) was adopted using the plane strain approximation. The seismic moment magnitude is derived from the scalar moment using standard relation (e.g. Stein & Wysession 2003) expressed in SI units: The scalar moment in Fig. 4(a) correlates with the maximum plastic shear strain at the borehole . The seismic source region is localized around the borehole and elongated in the direction orthogonal to the maximum compression. The numerical values of seismic magnitude M w in our model depend on the area of plastic zone and scale with the model size as R 2 0 (Fig. 5). The range of event magnitudes is obtained based Downloaded from https://academic.oup.com/gji/article/227/1/33/6289921 by guest on 08 July 2021 (a) (b) Figure 4. Normalized magnitude of the moment tensor for non-hydrostatic compression in isotropic medium (a) and transversely isotropic (n = E 1 /E 2 = 2, ν 1 = ν 2 = 0.3, μ 2 = E 2 /(2(1 + ν))) medium (b). The seismic event with largest moment magnitude is shown by the target symbol. The displacement field is shown by arrows. The spatial variability of the seismic source mechanism associated with the borehole failure can be analysed using Hudson source type parameters τ and κ (introduced in Section 3.3). The non-zero values of these parameters correspond to the regions of plastic deformation around the borehole (Fig. 6). The parameters τ and κ are inversely correlated, and the magnitude of isotropic extension κ > 0 is proportional to the CLVD component. The region of τ and κ with smallest magnitudes corresponds to the maximum moment magnitudes M 0 (cf. Fig. 4a). The Hudson parameters attain their maximum absolute value at the inner radius, oriented by 45 • angle with respect to the maximum compression directions (Fig. 6).

Transversely isotropic dry rocks
The anisotropy of elastic properties is commonly observed in rocks during laboratory experiments and geotechnical operations. The anisotropy can affect both the elastic rock response to loading and the failure pattern (Choens et al. 2019). The elastic anisotropy can also contribute to non-DC source mechanisms (Vavryčuk 2005). In Figs 7 and 4(b), we present results for rocks with preferred orientation of elastic properties. We choose a transversely isotropic elastic medium where the directional dependence of elastic properties is caused by stratification. Our y-axis in Fig. 2 is the rotational symmetry axis and perpendicular to the strata. The model stratification plane is parallel to x−z plane where z is the out-of-plane axis. The expression for transversely isotropic elastic stiffness matrix D in plane strain given by Zienkiewicz & Taylor (2005) is where n = E x /E y , m = μ y /E y ; E x , ν x are the elastic moduli in plane of the strata, and E y , ν y and μ y moduli are normal to the plane. In order to study the effect of anisotropy, we apply incrementally only external isotropic pressure (maximum mean stress is 60 MPa); so that, the deviation from the radial symmetry of the failure pattern around the borehole should only depend on the elastic properties.
The anisotropy is implemented using the elastic stiffness matrix (55). The Young modulus in the strong (x) direction is set twice as in the weak (y) direction (n = E x /E y = 2) whereas the Poisson ratio is taken constant (ν y = ν x = 0.3).
The elastic anisotropy results in elongation of the deformation region in the weak y direction (Figs 4b and 7). The failure pattern is generally the same as in the case of applied far-field shear (Fig. 4a). However, in the anisotropic case, some strong events are also obtained on the western (left) side of the borehole wall. Both τ and κ parameters change the sign (Fig. 7), and the predicted source mechanism implies a significant implosion component along the weak y-direction.

Isotropic fluid-saturated rocks
The third modeling setup corresponds to a combined stress and fluid pressure loading of isotropic poroelastic rocks. We extended the model for dry isotropic rocks described in Section 4.2 by implementing equations in Section 2.5. The solution of decoupled poro-elastoplastic problem requires two main modifications compared to the dry case: (i) adding the body force term in the force balance equations caused by the fluid pressure gradients (eq. 26); (ii) adding the fluid pressure in the failure criteria eqs (29) and (31) to counteract the effect of confining pressure. The far-field shear stress boundary condition was the same while the isotropic stress was three times smaller compared to the dry isotropic case. The expression for fluid pore pressure (eq. 28) allows to derive an analytical form of the body force. Our elastoplastic model is independent of time. Therefore, we choose the position of the fluid diffusion front to be constant, and equal to the twice of the model's inner radius (l d = 2R 0 ). Following our incremental loading approach, the fluid pressure at the borehole wall (p 0 ) increases with a small step during the simulation towards its maximum value (see Table 1). In our setup, the added fluid pressure inhibits the shear failure resulting in overall smaller plastic strain at the borehole (Fig. 8a). The Mode-II failure (shear events) pattern is similar to the dry case. As expected, Mode-I failure (tensile events) appears in the direction of maximum compression stress which is in the x-direction in Fig. 8. The events with largest magnitude are located at the borehole wall. In addition, the shear failure at random heterogeneities induced by fluid pressure results in a number of small-magnitude events away from the borehole. These microevents follow the shear deformation bands along the elevated shear stress directions (45 o angle). With progressive loading, the deformation bands may eventually develop into fault planes with the orientation controlled by effective pressure and friction angle. Geological observations indicate that faults may form in or along deformation bands in porous rocks (Fossen 2020).

Hudson plots and radiation patterns
We project the obtained source mechanisms on a standard Hudson τ -κ diagram for the events with M 0 larger than 0.1 per cent of the maximum value (Fig. 9). The domain of the calculated moment tensors is known as shear-tensile failure (Vavryčuk 2011;Fischer & Guest 2011;Jia et al. 2018). We use eq. (43) together with eqs (46)-(51) to calculate a theoretical dilatant shear crack curve in Fig. 9. The location of a source mechanism along this curve depends on the relative magnitude of the shear and bulk strain components in eq. (42) and the dilation angle ψ. Most events project onto the straight line from '+Crack' to '-Crack'. The slope of this line is determined by the Poisson ratio ν. The '+Crack' and 'DC' points correspond to two fundamental cracking mode in rock mechanics. The opening-dominated fracture would be located at the '+Crack' point whereas shear fracture would be at the centre of the diagram (double couple, 'DC').
The events in the isotropic numerical model shown in Figs 6 and 9(a) occupy a compact region along the theoretical curve in quadrant II of the plot. The reulst shown in Fig. 9(b) correspond to the model where we superimposed the random cohesion field ( Fig. 2b) and a random Poisson ratio field assuming the same statistical properties of the small-scale heterogeneity (see Table 1). From eqs (6), (46), (47) and (52), the variation of Poisson ratio leads to dispersion along the axis of the isotropic source component κ. This result is in agreement with our numerical simulations in Fig. 9(b). The results in Fig. 9(c) correspond to the anisotropic numerical model, also shown in Fig. 7. In this case, the source parameters occupy a wider region along the theoretical curve in quadrants II and IV. The obtained source mechanisms span from explosive crack across the DC region in Fig. 9(c) and towards the events with an implosive crack component. The model including the fluid pressure predicts both tensile events concentrated at the '+Crack' point (mode-I crack) and shear events distributed similarly to the dry isotropic case (Fig. 9d). The amplitude of seismic waveforms, measured at a large distance, is proportional to the intensity of wave radiation at the source. We calculate the seismic energy radiation pattern due to plastic failure using eq. (44). The radiation pattern in Fig. 10 corresponds to the grid element with the largest magnitude M 0 . Both P and SH waves are radiated. In isotropic model (Fig. 10a), the P-wave radiation is dominated by the tensile source component with smaller compression lobes. In the anisotropic case (Fig. 10b), the compressional event has an opposite sign of the P-wave radiation whereas the SH-radiation (Fig. 10c) is not affected. This observation can potentially be used as a diagnostic tool in real data applications.
The stereographic projection of the P-wave radiation on the focal plane is shown in Fig. 11. The beachballs for the isotropic and anisotropic case are very similar, and most likely would be very difficult to distinguish in practice. The orientation of the stress directions (pressure 'P' and tension 'T' axes in Fig. 11a) based on isotropic model coincides with the orientation of the shear stresses applied at the external boundary of the model indicating that the focal mechanisms should provide a good estimate of tectonic stress directions. In transversely isotropic medium (Fig. 11b), 'P'-and 'T'-axes do not reflect the boundary conditions but result from the anisotropic rock properties: the strong and the weak elastic directions, respectively.

Statistics of seismic events
The relation between the frequency of occurrence and the seismic moment of tectonic earthquakes generally follows the Gutenberg-Richter (GR) relation (Gutenberg & Richter 1944;Stein & Wysession 2003). This relation corresponds to a cumulative power-law distribution with exponent equal −1. In earthquake seismology, the b-value is defined as the negative of the power-law exponent. We analysed the statistics of seismic events with the progress of deformation in our model (Fig. 12, left-hand column). The cumulative distribution of normalized moment magnitude (log 10 M 0 /R 2 0 − 5) is compared to the GR relation on four model snapshots in pseudo-time corresponding to the incremental loading (Fig. 12, middle column). The maximum likelihood estimate (MLE) of the b-value can be obtained for a subset of logarithmic moment magnitude data (M) above the threshold value (M c ) (Aki 1965 whereM is the mean magnitude. M c has been selected as the magnitude with the maximum probability density in the dataset (Wiemer & Wyss 2000;Eaton 2018). The estimated b-value (b MLE ) shown on a log-log plot in Fig. 12 (middle column) varies from 0.96 to 1.8. The power-law distribution fits the simulation data only in the medium-magnitude range. The observed cumulative distribution saturates at low magnitudes and falls-off more rapidly than GR for large-magnitude events.
The microseismicity data, such as obtained during hydraulic fracturing monitoring (Eaton et al. 2014a;Eaton 2018) and enhanced geothermal reservoir operations (Dorbath et al. 2009;Leary et al. 2020), often show similar deviations from the GR relation towards a distribution characterized by a smaller probability for both small-and large-magnitude events. This type of behaviour can be approximated by a probability density combining an exponential and power-law functions (Kagan 2002). In particular, the Weibull distribution (Weibull et al. 1951) belongs to this class of functions, and has been used to approximate the temporal statistics of seismic events (Abaimov et al. 2008).
In Fig. 12, we compare the probability density distribution for simulated seismic moment data with different parametric statistical distributions fitted to the observations. The statistics of the seismic events evolves with the progress of deformation (from top down in Fig. 12; see also Table 2). At early stages, different distributions can approximate the observations including the power law, Gaussian, exponential and Weibull models. When the deformation is well developed, the power-law and Gaussian models fail to approximate data, whereas the exponential and Weibull distribution provide best fit to the simulations (r-squared ≥0.98). The role of material parameters, stress heterogeneity   and other factors controlling the statistics of seismic events is a topic of ongoing research, and can be addressed in the framework of statistical mechanics (Heimpel 1996;Rundle et al. 1997;Ampuero et al. 2006;Vallianatos et al. 2016;Leary et al. 2020).

Non-DC source mechanisms
The elastoplastic formulation for the seismic source mechanism (eq. 6) naturally includes both a shear and isotropic component, and can be applied for interpretation of microseismic events with a non-DC source mechanism. The moment density tensor is a product of the tensor of elastic constants and the plastic strain tensor. Thus, the source parameters must be controlled by both of these quantities. Our model predictions in Fig. 9 generally agree with previous publications in that a variety of the synthetic source mechanisms are all located along the shear-tensile failure domain on the Hudson diagram [+Crack to -Crack line, Julian et al. (1998)]. The tensile-shear failure model in application to fluid injection operations have previously been discussed by Fischer &Guest (2011) andJia et al. (2018). Our interpretation regarding the presence of both tensile and shear events during hydraulic injection is somewhat different from previous works (e.g. Šílený et al. 2009;Fischer & Guest 2011;Vavryčuk 2011) who attributed the presence of isotropic components to opening of pre-existed fractures with preferred orientation along σ 1 -direction. In our continuum mechanics model, both the mode and geometry of the rock failure is obtained in numerical simulations based on plasticity theory. The isotropic moment tensor component in shear failure regime arises naturally due to plastic dilatation but also can be related to elastic anisotopy. In addition, both mode-I and mode-II fractures can nucleate simultaneously in certain regions of the model depending on the distribution of stresses and fluid pressure. The source parameters corresponding to the two cracking modes are distinct and clearly identified on the Hudson plot in Fig. 9d.
In the isotropic dry rock regime (Fig. 9a), the event location along the +Crack to -Crack line is determined by the amount of isotropic plastic strain. In our numerical experiments of wellbore failure, the larger dilation angle ψ implies κ and τ which are closer to the explosive shear crack while the lower ψ imply a larger DC source component. The case of zero plastic dilation corresponds to a DC source mechanism where τ and κ are both zero for all events around the borehole. The moment tensor parameters associated with the zero dilation angle plot close to the centre of the Hudson diagram.
The borehole failure and hydraulic fracturing in the lab are often associated with some implosive events (Baig & Urbancic 2010;Aker et al. 2014;Vera-Rodriguez et al. 2017). This observation disagrees with our result of non-negative plastic dilation implied by eq. (35). Implosive earthquakes imply a negative trace of the deviatoric moment tensor: Since the Poisson ratio is positive, this simply means that M 1 < −M 3 . The moment tensor is a product of the plastic strain and the elastic stiffness tensor. We can relax the assumption of elastic isotropy and find the difference of elastic moduli in strong and weak directions (transverse isotropy) required to produce a micro-earthquake with an implosive source mechanism.
Consider a diagonal form of the moment tensor which we find by multiplying anisotropic stiffness matrix (55) with the diagonal form of the plastic strain tensor. We assume here that D 11 > D 22 , and M 1 is positive for extension. Thus, the term in the first parentheses in eq. (57) can be written Using this relation, we predict three classes of the seismic source mechanisms, depending on the combination of elastic moduli D and plastic dilation angle ψ: η > sin ψ, Implosion, tr(M) < 0 (59) η = sin ψ, Double couple, tr(M) = 0 (60) where we defined η as The borehole breakouts in dry transversely isotropic rocks may result in pure DC and implosive events, as well as tensile events (Fig. 9), despite the plastic dilation is always positive. Note, that this result applies to a plastic shear failure around wellbore. The plastic compaction and corresponding implosion source mechanism may occur if stresses reach elliptic cap surface. This plastic regime was suggested for high pressures and is generally referred as the shear-enhanced compaction regime (Fortin et al. 2009).
The presented plane strain model provides useful analytic relations but it can be applied strictly to problems with cylindrical geometry. The general theory linking the seismic moment tensor source with the plastic strain tensor (described in Section 2) can still be used for complex 3-D geometries. The theoretical analysis in 3-D becomes less transparent, and purely numerical solutions will probably be preferred.

Interpretation of laboratory acoustic emission data
Our model can explain some recent observations made in the lab (Aker et al. 2014;Vera-Rodriguez et al. 2017) and hydraulic fracture stimulations (Baig & Urbancic 2010). The majority of microseismic events are located in quadrants II-IV of the Hudson plot. This directly follows from the representation (43) of the seismic moment tensor combined with the shear and pressure loading. For constant elastic parameters, the Hudson parameters plot along the main diagonal. Since the plastic flow defined by the yield criterion eq. (24) is dilational, the obtained source mechanism is a combination of explosive crack and DC components.
Another interesting observation made in the borehole breakout experiments (Aker et al. 2014) is that the events with a large moment magnitude have a small volumetric component. This can be explained in terms of distribution of principal stresses and principal moments due to stress concentration around the borehole. The results of our dry rock simulations (Figs 4a and b) indicate that the minimum and maximum principal moments have the opposite sign in the failure regions (M 1 > 0,M 3 < 0). This implies that M 0 (in eq. 53) is larger for a smaller isotropic source component (M 1 +M 3 , Figs 4a, b and 7).
The laboratory acoustic emission data suggest that percentage of DC component increases with progressive loading (Aker et al. 2014). This can be explained by a decrease of plastic dilation ψ with time and loading observed in deformation experiments (Davis & Selvadurai 2005). Note that both the bulk and shear plastic strain are largest close to the borehole (Fig. 3). Therefore, the events close to the borehole should have a larger isotropic component compared to more distal events.

Implication for crustal stress measurements
Diagnostic observations on in situ stress state, such as boreholes breakouts and earthquake focal mechanisms, are important for understanding geological evolution of the lithosphere and for practical geotechnical applications (Zoback 1992;Amadei & Stephansson 1997). A robust determination of crustal stress from deep borehole data requires integration of various methods (Schmitt et al. 2012). In isotropic rocks, the regions of shear failure at the borehole wall are located along the minimum tectonic stress direction. Our numerical models indicate that the presence of elastic anisotropy results in a similar failure pattern under hydrostatic load (Figs 4a and b). Can we really distinguish the effect of tectonic stresses from the intrinsic elastic anisotropy? In our model, the microseismic source mechanisms alone cannot provide a definite answer. The numerical model results shown in Figs 6 and 7 indicate that the effect of anisotropy on source parameters is much more pronounced compare to the failure pattern alone. The Hudson parameters τ and κ in the anisotropic case span a large range of values and even can change the sign (Fig. 9c). Thus, a joint analysis of microseismic moment tensor data and failure pattern around borehole can be used for better constraining the tectonic stress.
Another question is to what extent we can distinguish the effect of heterogeneity and anisotropy using microseismic observations? Indeed, most upper crustal rocks are both heterogeneous and anisotropic. The Poisson ratio variation results in spreading of seismic events along the κ-axis (Fig. 9) whereas the anisotropy would spread the events along the diagonal (+Crack-DC) direction. This is another diagnostic observation that can be potentially used by future high-resolution moment tensor inversion studies.

C O N C L U S I O N S
Plasticity theory is useful to extrapolate the laboratory acoustic emission data to microseismic observations of crustal deformation. We proposed a new representation of seismic sources derived from the non-associative plastic flow law. This model does not require a pre-defined fracture geometry. It helps predicting the localized failure pattern and corresponding seismic response. Using a quasi-static incremental formulation, plastic strains can be incorporated as a part of the effective body force in the wave propagation equation. The traditional definition of the moment tensor via the stress or strain 'glut', reformulated in terms of plastic strain tensor, can provide more information on the physical mechanisms inside a non-planar source region. The elastoplastic representation for moment tensor sources includes both shear and isotropic source components through the standard parametrization of the plastic flow potential, and can be used to explain the source mechanisms observed in hydraulic fracturing and acoustic emission experiments.
Downloaded from https://academic.oup.com/gji/article/227/1/33/6289921 by guest on 08 July 2021 The 2-D numerical examples apply to the borehole/tunnel stability problem or seismicity around a fluid-saturated vertical channel in porous rocks, and facilitate joint analysis of rock deformations and seismic moment tensor data. On the Hudson diagram, the source parameters predicted by the numerical model, spread along the '+Crack' to '-Crack' line known as the tensile-shear failure domain. The spread across the theoretical line is controlled by elastic parameters. The position of source mechanism along this line is determined by the plastic dilation. The fluid pressure in our model does not change this relation but can lead to a separate tensile failure region located at '+Crack' point.
The developed model can also be useful for better characterization of material parameters and for a more robust estimate of tectonic stresses using moment tensor solutions. The orientation of borehole breakouts in dry rocks depends on both tectonic stress and anisotropic elastic parameters. In dry transversely isotropic rocks, dilatant failure of a wellbore may lead to a tensile, DC or implosive source mechanism depending on the relation between the plastic dilation (sin ψ) and the parameter η defined in eq. (62). Spatial analysis of moment tensor solutions can reveal the relative contribution of material anisotropy versus tectonic stress.
The pore fluid pressure controls the preferred failure mode around a pressurized circular opening, and can lead to two distinct regions corresponding to the tensile and shear failure. In our poro-elastoplastic model setup, random heterogeneities of crustal strength can be seismically activated along the maximum shear stress directions due to the mechanical effects of fluid diffusion and seepage forces. Our approach is compatible with existing geomechanical and seismic wave propagation solvers and enables their closer integration for microseismic monitoring.

A C K N O W L E D G E M E N T S
We thank David Eaton, Eyvind Aker and two anonymous reviewers for constructive criticism and useful suggestions. For preparation of colour figures in this article, we used scientific colourmaps designed by Fabio Crameri (https://www.fabiocrameri.ch/colourmaps/). VY acknowledges funding by Research Council of Norway (Project 280953) and Ministry of Science and Higher Education of the Russian Federation(Project 075-15-2019-1890). AM acknowledges funding by Research Council of Norway (Project 223272, The Centre for Earth Evolution and Dynamics).

DATA AVA I L A B I L I T Y
The data and MATLAB scripts to produce figures underlying this paper are available in Zenodo.org, under doi: 10.5281/zenodo.4911325