Layout 6


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

Dimitrios Pliakis1, Taxiarchis Papakostas1, Filippos Vallianatos1, 2, *

1 Technological Educational Institute of  Crete, Laboratory of  Geophysics and Seismology, Crete, Greece
2 University College London, Department of  Earth Sciences, London, United Kingdom

ANNALS OF GEOPHYSICS, 55, 1, 2012; doi: 10.4401/ag-5363

ABSTRACT

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

1. 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 and
Turcotte 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 and Sornette 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 and Sammis 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-to-
failure relation of any quantity estimated from the
earthquake magnitude. This scaling law has been justified in

Special Issue: EARTHQUAKE PRECURSORS

165

Article history
Received August 10, 2011; accepted February 1, 2012.
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 and Varnes 1993, Vallianatos and
Triantis 2008]. However, it can also result naturally from the
many-body interactions between small cracks that form
before an impending rupture [Sornette and Sammis 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 non-
singular power-law relation for cumulative Benioff strain that
is proportional to (tc – t)

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

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

S x
u

x
u

2
1

kl
l

k

k

l

2

2

2

2
= +c m

SD klij ijkl=v

div ( ) x 0j
jk

2

2
v

v
= =

,I l ldl
0

j ct j=
3

^ ^h h#

166

PLIAKIS ET AL.



167

ON THE PHYSICS OF ACCELERATING SEISMICITY

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_ =
=(j1, j2). 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 c01, c02 > 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 = i1i2i3 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 M0 be a domain in space and M an arbitrary
subdomain with a smooth boundary. Then we denote by Mt
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:

C I n n n n d
D

ijkl i j k li X= ^ h#

( )D I Sv = u

( ) ( ( ))D I DC I D1 1= + -u

( ) ( , 0)L I x d dM
D

x
M

v y= ##

( ) ( ( ))C C I n2
1

M
D

ijrs ij rs
2

M M

v v x= =# ##

( ) ( , ( ))n n nx v=

( )MLC In n n n d d Id d
D D

ijkl i j k l x x
M M M

#= =v y v y# ## ##

vol ( ) curl ( ) vol ( )c C CC cM M
e

` `
e

01 02# #

curl ( ) ( )C x
C

x
C

n A I di
i i

D

jk
k

j

i i i jk
j

k

1 2 32

2

2

2
X= - = #

( )A I n x
I

n x
I

jk j
k

k
j2

2

2

2
= -

(vol ( ))max c CM
e

10
M

M

3`
2

2
3

#v
t

a k#

(vol ( ))max c CM
e

11
M

M

`

2
3 2

3

$v
t

a k#



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:

3. 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 (Tf – t)

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 and Papazachos

2001, Vallianatos and Tzanis 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.

Acknowledgements. This study was supported in part by the
ARCHIMEDES III Programme of the Ministry of Education of Greece,
and the European Union in the framework of the project entitled
‘Interdisciplinary Multi-Scale Research of Earthquake Physics and
Seismotectonics at the Front of the Hellenic Arc (IMPACT ARC)’.

References
Ben-Zion, V. and V. Lyachovsky (2002). Accelerated seismic

release and related aspects of seismicity patterns on ear-
thquake faults, Pure Appl. Geophys., 14, 2385-2412.

Bowman, D., G. Ouillon, C.G. Sammis and D. Sornette
(1998). An observational test of the critical earthquake
concept, J. Geophys. Res., 24, 372.

Brehm, D.J. and L.W. Brail (1999). Intermediate term ear-
thquake using the modified time to failure method in
Southern California, B. Seismol. Soc. Am, 89, 275-293.

Brehm, D.J. and L.W. Brail (1998). Intermediate-term ear-
thquake using precursory events in the new Madrid sei-
smic zone, B. Seismol. Soc. Am., 88, 564-580.

Bufe, C.G. (2004). Comparing the November 2002 Denali
and November 2001 Kunlun earthquakes, B. Seismol.
Soc. Am., 94, 1159-1165.

Bufe, C.G. and D.J. Varnes (1993). Predicitive modeling of
the seismic cycle of the greater San Fransisco Bay re-
gion, J. Geophys. Res., 98, 9871-9883.

Cowie, D., P.A. Vanneste and T. Sornette (1995). Multi-
fractal scaling properties of a growing fault population,
Geophys. J. Int., 122, 457-469.

Gilbarg, D. and N. Trudinger (2003). Elliptic Equations of
the Second Order, Springer Verlag.

King, G.C.P. and D. Bowman (2003). The evolution of re-
gional seismicity between large earthquakes, J. Geo-
phys. Res., 108, 2096-2106.

( )
( ) ( )MC

L

dt
d c

n I
c

2 2
M

D

2

M

= =
r

x
r

6 @##

C

dt
d

C
e

M

2
3

+ fa k#

vol ( ) vol ( )C
C

C CM
M

3 1 3 1
3 1

`
e

`

e

3 1

1

&+ +
+

+

+
+

c m

C
Cdt

d
1
1 4

e

e

3+ +
+

( )C C c T tf f2 3
1

= - -

.n I d I
D

3

1

#4 +X -#

PLIAKIS ET AL.

168



169

ON THE PHYSICS OF ACCELERATING SEISMICITY

Landau, L.D. and E.G. Lifshitz (1936). Theory of Elasticity,
Pergamon Press.

Main, I.G. (1999). Applicability of time-to-failure analysis
to accelerated strain before earthquakes and volcanic
eruptions, Geophys. J. Int., 139. F1-F6.

Mignan, A., G. King and D. Bowman (2007). A mathema-
tical formulation of accelerating moment release based
on the stress accumulation model, J. Geophys. Res.,
112, 307-308.

Papakostas, T. and D. Pliakis (2011). Local existence of
Euler-Einstein system, Growth estimates of this type
(in preparation).

Papazachos, B.C. and C.B. Papazachos (2001). Precursory
accelerated Benioff strain in the Aegean era, Annals of
Geophysics, 44, 461-471.

Pliakis, D. (2002). Generalized Hardy’s Inequality, arXiv:
math/0203092v2 [math.AP].

Pliakis, D. (2011). On the volume of nodal sets. Harnack
inequalities in general domains (submitted)

Pliakis, D. and S. Minardi (2009). Phase front retrieval by
means of an iterative shadow graphic method, J. Opt.
Soc. Am. A, 26, 99-107.

Rundle, T., W. Klein, D.L. Turcotte and D.B. Malamud
(2000). Precursory local existence of einstein-euler sy-
stems, Pure Appl. Geophys., 157, 2165-2182.

Sammis, C.G. and D. Sornette (2002). Damage and self si-
milarity in fracture, Proc. Nat. Assoc. Sci. USA, 99,
2501-2508.

Scherbakov, R. and D.L. Turcotte (2003). Damage and self
similarity in fracture, Theor. Appl. Mechan., 39, 245-
258.

Sornette, D. and C.G. Sammis (1995). Complex critical ex-
ponents from renormalization group theory of ear-
thquakes: implications for earthquake predictions, J.
Phys. I., 5, 607-619.

Taguchi, Y. (1989). Fracture propagation governed by the
laplace equation, Physica A, 156, 741-755.

Turcotte, D., W. Newman and R. Scherbakov (2003). Micro
and macroscopic models of rock fracture, Geophys. J.
Int., 152, 718-728.

Vallianatos, F. and P. Triantis (2008). Scaling in pressure sti-
mulated currents related with rock fracture, Physica A,
152, 718-728.

Vallianatos, F. and A. Tzanis (2003a). Distributed power law
seismicity changes and crustal deformation in the sw
Hellenic arc, Nat. Hazard Earth Sys., 3, 179-195.

Vallianatos, F. and A. Tzanis (2003b). On the nature, sca-
ling and spectral properties of pre-seismic ulf signals,
Nat. Hazard Earth Sys., 3, 718-728.

Vallianatos, F., D. Triantis, A. Tzanis, C. Anastasiadis and I.
Stavrakas (2004). Electric earthquake precursors: from
laboratory results to field observations, Phys. Chem.
Earth, 29 (4-9), 339-351.

Varotsos, K., P. Alexopoulos and M. Lazaridou (1993). La-
test aspects of earthquake prediction in Greece based
on seismic electric signals II, Tectonophysics, 224, 1-37.

Wackertapp, S., J. Hergarten and H. J. Neugebauer (2000).
A model for propagation and concentration of micro-
cracks, Geophys. Res. Lett., 27, 1771-1774.

* Corresponding author: Filippos Vallianatos, 
Technological Educational Institute of Crete, Laboratory of Geophysics
and Seismology, Crete, Greece;
email: f.vallian@chania.teicrete.gr

© 2012 by the Istituto Nazionale di Geofisica e Vulcanologia. All rights
reserved.



Appendix
The estimates indicated in the main text will be achieved

here with the Moser iteration technique [Gilbarg and
Trudinger 1978], with the ‘tricks’ presented by Pliakis and
Minardi [2009], Pliakis [2011], and Papakostas and Pliakis
[2011]. We summarize them here. The building blocks are:

1. The Sobolev inequality [Gilbarg and Trudinger 1978]
in R3; q = 3p/(3 – p), p < 3 reads as for f ! C0

∞ (R3):

2. Generalized Hardy inequality [Pliakis 2002] for a

polynomial function in P: Rn $ R – here n = 3 we have that
for f ! C0

∞ (Rn\{P = 0}):

Applying Gauss theorem, we see that for an arbitrary
smooth cut-off | and a smooth (p + 1)-symmetric tensor
field U: 

where i = i1… ip,:

and the magnitude of a (p + 1)-tensor:

Choosing a smooth cut-off {: | = |U|k – 1�, we are led
for U = v to the following inequality for p < 2, q = 3p/(3 – p)

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

f C f
q p

4#

( )P
P

f C P f
2

2 2

M M

4
4## #

) ( )div ( curlU C U U U2
2 2 2 22

M M

4 4#| |+ +# #

div ( ) , ( )curlU
U

U x
U

x
U

xi
i

i

i i
j

j

j

jk

j

k j

k

2

2

2

2

2

2
= = -

.U U ,...
,...

i i
i i

2 2

p

p

1 1

1 1

=
+

+

/

( )

( )

, curl ( )

C

C N

/ /

2( 1)

qk q p p

k

k

p
p

1 1

2
2

M M

M

4# #

#

{v {v

f { v v
-

-

p

a a

a

k k

k

# #

#

C S S Dijrs rs ij ijrs rs,v v= =

=0,C C C
2

3 =
2X

t t

=0, DD D
2

3 =
2X

t t

PLIAKIS ET AL.

170