A first principles approach to understand the physics of precursory accelerating seismicity

Observational studies from rock fractures to earthquakes indicate that fractures and many large earthquakes are preceded by accelerating seismic release rates (accelerated seismic deformation). This is characterized by cumulative Benioff strain that follows a power law time-to-failure relation of the form C(t) = K + A(Tf – t)m, where Tf is the failure time of the large event, and m is of the order of 0.2-0.4. More recent theoretical studies have been related to the behavior of seismicity prior to large earthquakes, to the excitation in proximity of a spinodal instability. These have show that the power-law activation associated with the spinodal instability is essentially identical to the power-law acceleration of Benioff strain observed prior to earthquakes with m = 0.25-0.3. In the present study, we provide an estimate of the generic local distribution of cracks, following the Wackentrapp-Hergarten-Neugebauer model for mode I propagation and concentration of microcracks in brittle solids due to remote stress. This is a coupled system that combines the equilibrium equation for the stress tensor with an evolution equation for the crack density integral. This inverse type result is obtained through the equilibrium equations for a solid body. We test models for the local distribution of cracks, with estimation of the stress tensor in terms of the crack density integral, through the NashMoser iterative method. Here, via the evolution equation, these estimates imply that the crack density integral grows according to a (Tf – t)0.3-law, in agreement with observations.


Introduction
The failure of a solid due to microcrack concentration and propagation has been studied in various disciplines, including material sciences [Scherbakov and Turcotte 2003], rock mechanics [Turcotte et al. 2003], solid earth geophysics [Main 1999], and seismology.
It is well accepted that there is a wide range of problems in geophysics where it is necessary to understand how rock deforms, from earthquake prediction to the driving forces of plate tectonics [Main 1999].The rock physics approach to understanding these geophysical processes is based on the premise that the macroscale behavior of rock is governed by microscale interactions.Rock deforms elastically and plastically under an applied stress, by fracturing, brittle flow, and frictional sliding on a fault.The magnitude and direction of the applied stress, the rate and duration of the loading, the ambient pressure and temperature, the presence of fluids, and the previous deformation history, all define the overall mechanical response [Main 1999, Scherbakov andTurcotte 2003].
Fracturing phenomena have attracted the interest of the scientific community, and particularly those concerning inhomogeneous materials.The main reason is that such phenomena are promising candidates as earthquake precursors [Vallianatos and Tzanis 2003a,b, Vallianatos et al. 2004, Varotsos et al. 1993, Vallianatos and Triantis 2008].
Over the past decade it has been credibly argued that rupture in heterogeneous media is a critical phenomenon [Cowie et al. 1995, Sammis andSornette 2002], while a mounting body of evidence indicates that the earthquake generation process can be viewed as a critical phenomenon that culminates in a large event that corresponds to some critical point (CP) [Sornette andSammis 1995, Rundle et al. 2000].According to the CP earthquake hypothesis, a failure in the crust can be thought of as a scaling up process, in which a failure on one scale of a fault network is part of the damage accumulation over a larger scale, which leads to long-range stress-stress correlation and an increase (acceleration) in the seismic release rates prior to large earthquakes [Vallianatos and Tzanis 2003b].The culminating large event will only occur when the network is in a critical state, balancing on the verge of disorder, and characterized by both extreme susceptibility to external factors and strong correlation between its different parts.
Such a process can be described by a power-law time-tofailure relation of any quantity estimated from the earthquake magnitude.This scaling law has been justified in Special Issue: EARTHQUAKE PRECURSORS

Subject classification:
Earthquake faults: properties and evolution, Seismicity fracture Benioff.
the laboratory in terms of run-away crack propagation and empirical expressions for accelerating (tertiary) creep preceding failure [Bufe andVarnes 1993, Vallianatos andTriantis 2008].However, it can also result naturally from the many-body interactions between small cracks that form before an impending rupture [Sornette andSammis 1995, Bowman et al. 1998].
Previous studies have determined that the cumulative Benioff strain [Vallianatos and Tzanis 2003b] is particularly useful when smaller events are also of interest and magnitude scaling is desirable, while the cumulative moment is dominated by larger earthquakes and the event count does not allow for magnitude scaling.Earlier studies empirically determined that typically the cumulative Benioff strain critical exponent is m ~0.3.There are also two theoretical predictions for the value of the critical exponent.Using a model seismogenic crust governed by a particular damaged rheology, Ben-Zion and Lyachovsky [2002] derived a nonsingular power-law relation for cumulative Benioff strain that is proportional to (t ct) 0.3 , i.e. m = 0.3.Rundle et al. [2000] adopted a very different approach, by relating the behavior of seismicity prior to a large earthquake to the excitation in the proximity of a spinodal instability.They showed that the power-law activation associated with the spinodal instability is essentially identical to the power-law acceleration of Benioff strain observed prior to an earthquake, and that m = 0.25.
Previous studies have developed and reviewed techniques to identify regions of accelerating seismicity and to attempt the prediction of the time and magnitude of the next large earthquake, sometimes with remarkable success [Bufe and Varnes 1993, Bufe 2004, Sornette and Sammis 1995, Brehm and Brail 1998, 1999, Bowman et al. 1998, Papazachos and Papazachos 2001].We note, however, that until recently, the CP model has been largely conceptual and based on an analogy with phase transitions, and it has drawn support from theoretical simulations that involve simple models, such as cellular automata.There have been no effective physical models to describe the observations and the evolution of the earthquake cycle.
To study the physics of seismicity in a similar way to the fracture of a solid, numerous approaches have appeared that have dealt with the propagation of cracks, and furthermore, with static [Taguchi 1989] and dynamic [Wackertapp et al. 2000, Mignan et al. 2007] distributions of microcrack concentration.
In the dynamic distributions of microcrack concentration, it is assumed that the space-and-timedependent microcrack concentration is governed by the applied microscopic external stress field.These have been studied according to the deviation of microcrack distributions from random ones during the deformation of a solid, and particularly how the microcrack concentration increases, and then the cracks coalesce and localize, and finally the solid fails [Scherbakov and Turcotte 2003].
In the present study, based on the approach presented by Wackertapp et al. [2000], we formulate the time dependence of accelerated deformation based on first principles.In Wackertapp et al. [2000] they propose a model for the propagation and concentration of microcracks in brittle solids.First, they derive the equations for the time evolution of an integral quantity related to the crack density.Apart from a factor that characterizes the material, the rate of growth of the crack density integral is determined by the stress tensor.Therefore, we apply partial differential equation techniques in an elastic equilibrium system.In this way, we estimate locally the stress tensor in terms of the elastic tensor or the crack density integral, provided the latter shows a local variation that we can determine to recover the observed laws.

First principles and methodology
To set up the notation, we have to review some standard concepts for the formulation of an elastic equilibrium system.From the mathematical physics point of view, we will focus our attention on an open domain W 0 in space that contains the origin and is equipped with the Euclidean metric denoted by (• , •).The latter description is equivalent with an earth that contains the earthquake zone.The elastic constant tensor field is a rank 4 tensor field, where D defines at each point in W 0 a positive definite quadratic form on the space of rank-2, symmetric tensors [Landau and Lifshitz 1936].Furthermore, if u _ is the deformation vector field, then we have the strain tensor field S : where i, j = 1, 2, 3. Then through the elastic constant tensor D and Hook's law, we introduce the stress tensor field v: where we use the summation convention.The equilibrium equations for the elastic medium are then expressed for k = 1,2,3 as: with the forces exerted on the boundary of the region W 0 giving rise to Dirichlet boundary conditions in the form of the prescribed deformation of its boundary.We extend in three dimensional space the studies of Wackertapp et al. [2000] that derived a system of an ordinary differential equation coupled to the elliptic system of elasticity, introducing the crack density integral: (1) where the c-factor has the dimension of inverse force, and t is the crack density as a function of the crack length l, where the crack orientation is described by the pair of angles j _ = =(j 1 , j 2 ).We refer to this as the Wackertapp, Hergarten and Neugebauer (WHN) model [Wackertapp et al. 2000].Furthermore, in this WHN model, the strain tensor is split into two components, one of which accounts for the contribution of the cracks through the tensor field C, which depends on the fourth-order moment of the crack density integral, as: where dX is the elementary solid angle, and where we integrate over the space of all directions in 3-space, denoted as D. We conclude that the elastic tensor is modified as: (2) (3) According to the WHN model, using I(i _), the strain in the domain W is given by: In addition, the cumulative strain stored in the domain M is expressed by: where we have introduced the tension field as a function of the direction field n _: Applying the triangle inequality, we conclude that: We proceed now to introduce the following assumption on the local distribution of cracks in a spatial region M for exponents `, e and positive constants c 01 , c 02 > 0: (4) where vol(M) is the volume of the domain M.This introduces the bounds on the local variation of the tensor field C, implying smooth limits in the physical variation of the crack population properties.For the index i = i 1 i 2 i 3 we set: and thus taking into account the radial distribution of the cracks.
We proceed now to determine the exponents `, e, to obtain the condition under which the observed law holds.This is accomplished through examination of the elasticity equation and the crack or strain propagation equation.Specifically, we divide our domain into regions where inequalities of type of Equation ( 4) hold, and then we estimate the variation of the strain tensor through the stored energy inside this domain.Moreover, this estimate completes the elasticity system as a 'div-curl' system that allows for the estimation of the local variation of the strain in terms of the local concentration of cracks.The standard Nash-Moser technique was modified by Pliakis [2011]: it is combined with the local harmonic approximation method, which provides local weights for the derivation of the estimates.Therefore, the local variation of the stress tensor is given by this elliptic system.In the Appendix we recall the standard identities and the standard arguments that allow us to derive the necessary estimates.Additionally, standard arguments provide the existence of a solution under a smoothly varying crack density integral.Wackertapp et al. [2000] provided numerical evidence for the behavior of the tension field x(n _).We proceed now to formulate the bounds of the physical quantities involved in our analysis.For this, we use the typical Nash-Moser iteration used by Pliakis and Minardi [2009] and Pliakis [2011] for domains defined by polynomials.
Specifically, we have the following result from following the arguments in Pliakis [2011]: Let M 0 be a domain in space and M an arbitrary subdomain with a smooth boundary.Then we denote by M t the points in M at distance at least t > 0 from the boundary.Then the stress field inside M at every moment shows the following estimates of its maximum and minimum values: The latter implies that in a stressed region of the earth we can define a subdomain in which the stress field can be bounded with upper and lower bounds scaled to the volume and the tensor field C that depends on the crack density integral.To define the evolution of the system, we proceed now to study the dynamics of crack concentration.Wackertapp et al. [2000] suggested that the crack density integral varies according to the cumulative strain stored in the region due to the presence of cracks: # and therefore the growth of the cumulative strain is scaled as in the expression: (5) Furthermore using the above estimates we see that: To derive a similar expression with the observed accelerated deformation law, after substituting C and integrating, Equation ( 5) should then be written as: This latter equation is in perfect agreement with that proposed by Main [1999] to scale crack growth.In the case when (1+4e)/(1+3e) ~ -2 i.e. e ~ -3/10 we get that:

Conclusions and perspectives
The present approach leads to an estimate of the temporal growth of the cumulative strain that is governed by an integrodifferential equation, provided the crack density has a specific local variation expressed also in the form: This condition suggests that as the crack density increases, then their variation averaged in all possible directions vanishes: these tend to align in a specific direction (which is reasonable, because as cracks condense and align, they merge to produce a fracture).The cumulative strain tensor follows a (T ft) 1/3 law.The WHN model [Wackertapp et al. 2000] comprises the development of cracks in the elastic constants of the material (as brittle solids), and accordingly, the evolution of cracks is determined by a simple equation where the rate of increase of cracks is proportional to the existing cracks.However, the constant of proportionality depends on the elastic constants of the material.This leads to a nonlinear system of equations.This latter approach is consistent with the evolution process in an earthquake zone, where microfracturing increases during the preparation of a main event, which leads to an increase in moderate earthquake events [Papazachos andPapazachos 2001, Vallianatos andTzanis 2003b].This is also in agreement with the stress accumulation model proposed by King and Bowman [2003], where the accelerating energy release is due to a decrease in the size of a stress shadow that was left by one or more previous events.
To obtain this precise rate, we determine the dependence of the stress tensor on the crack density integral.This is derived from the elastic equilibrium equations of the medium, provided the crack density integral satisfies a particular form of local estimates, the precise form of which is set by the requirement that they reproduce the observed law.We apply the Moser iteration scheme for the elliptic system of elasticity, incorporating the local estimates of the crack density integral through weighted inequalities that was originated by Pliakis and Minardi [2009] and Pliakis [2011].Finally, the differential equation describing the microcrack density evolution leads to the 1/3-law.Furthermore, we expand the WHN approach to a new three-dimensional formulation.Introducing Hook's law: and using the harmonic approximation method for the norm of C, D: This brief description of the mathematical method allows us to produce the local weights that are necessary for the estimates.Then using localization in conjunction with Hardy inequality for the terms curl(C); curl(v) as well as Equation ( 4), after iteration, we derive the basic estimates of the stress tensor.The full treatment is given in Pliakis [2011] and Papakostas and Pliakis [2011].