Matsu’ura ﻿Progress in Earth and Planetary Science(2024) 11:2
https://doi.org/10.1186/s40645-023-00602-x

Progress in Earth andPlanetary Science

Open Access

REVIEW

Reconsideration of the energy balancein earthquake faulting
Mitsuhiro Matsu’ura1*   

Abstract
The occurrence of earthquakes is now understood as brittle shear fracture releasing the elastic potential energy
stored in the earth. Since the 1950s, many studies on the energy balance in earthquake faulting have been done,
but there seems to be some incoherence among them. The essential reason is because various changes in conceptual framework happened during the last six decades, specifically the introduction of the new paradigm of plate
tectonics in the 1960s, the concept of moment tensor as source representation in the 1970s, and the fault constitutive
law governing rupture growth in the 1990s. Therefore, it will be worthwhile to reconsider the energetics of earthquake faulting from a current perspective. For this purpose, first of all, we summarize the basic concepts of elastic
potential energy and moment tensor and review the general representation of earthquake sources and the originof background crustal stress to confirm that the effect of earth’s self-gravitation is negligible in the energetics of shear
faulting. Next, as a starting point for discussion, we directly derive a basic equation of mechanical energy balancein dynamic shear faulting from the equation of motion for an elastic body subjected to tectonic-origin deviatoricstress. Then, we review the widely accepted formula for indirectly evaluating radiated seismic energy from a simplified energy balance equation and compare with the direct evaluation based on the analytical solution of displacement fields for a point dislocation source in order to call attention to inconsistency between them. The inconsistency
comes from the omission of the effects of rupture growth rate in the simplified energy balance equation. So, finally,
we review the energy balance at the tips of a propagating shear crack, which naturally leads to the introductionof the slip-weakening fault constitutive law as a fundamental equation governing earthquake rupture. Then, we discuss the whole process of earthquake rupture, consisting of initiation, acceleration, steady propagation, deceleration,and termination from the viewpoint of energy balance.
Keywords Theoretical seismology, Background stress field, Earthquake faulting, Moment tensor, Inelastic strain,
Elastic potential energy, Radiated seismic energy, Work done for shear faulting, Fault constitutive law, Rupture growth
rate*Correspondence:
Mitsuhiro Matsu’ura
m.matsuurajpgu@gmail.com
1
School of Science, University of Tokyo, 7‑3‑1 Hongo, Bunkyo‑ku,
Tokyo 113‑0033, Japan© The Author(s) 2024. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Matsu’ura ﻿Progress in Earth and Planetary Science(2024) 11:2

1 Introduction
The occurrence of earthquakes is now understood as
brittle shear fracture at a fault, which releases a part ofthe elastic potential energy (elastic strain energy) stored
in the earth. Since the 1950s, many studies on the energybalance in earthquake faulting have been done, but there
seems to be some incoherence among them. The essential
reason is because various changes in conceptual framework happened during the last six decades, for example
the introduction of the new paradigm of plate tectonics
in the 1960s, the concept of moment tensor as source
representation in the 1970s, and the fault constitutive
law governing shear rupture growth in the 1990s. Therefore, it will be worthwhile to reconsider the energetics inearthquake faulting from a current perspective.
At the end of the 1950s, Steketee (1958) introduced theelasticity theory of dislocation into geophysics and considered the static change in elastic strain energy by theformation of dislocation (tangential displacement discontinuity = fault slip) in a pre-stressed elastic body bounded
by a traction-free surface. He concluded that the totalelastic strain energy is always increased by the formationof dislocation irrespective of the initial stress state, whichis called the Steketee’s paradox. At the end of the 1960s(after the introduction of the new paradigm of plate tectonics), Savage (1969) gave a geophysical interpretationof the interaction energy between two stress systems(Eshelby 1956) and demonstrated that the formationof dislocation always decreases the total elastic strainenergy if it occurs so as to produce a stress drop. However, such an argument seems to be out of focus, because
the Steketee’s paradox results from a misinterpretation of
his correct formula, the first line of Eq. (6.10) in Steketee(1958). Actually, with this formula as a starting point,
we can easily derive essentially the same expression as
Eq. (4) in Savage (1969). After all, as Andrews (1978)
discussed later, the fundamental question arising here is
how the pre-stressed state can be formed in the isolated
elastic body, the earth.
By the way, the earth is a self-gravitating body, and so
we should consider the change in gravitational potential
energy as well as the change in elastic potential energy. In
the case of self-gravitating earth models, as pointed out
by Kostrov (1974) and Dahlen (1977), the release of thetotal (elastic and gravitational) potential energy balanceswith the sum of the work done for shear faulting and theradiated seismic energy. In their formulations, the releaseof the gravitational potential energy is represented by thevolume integral of the scalar product of gravity force andcoseismic displacement over the entire earth. Savage andWalsh (1978) found that this volume integral could betransformed into the surface integral of the scalar product of gravity-origin shear stress and fault slip over the

Page 2 of 19

area of earthquake faulting. The point is that the isotropic
part of gravity-origin stress is independent of the work
done for earthquake faulting. Anyway, they evaluated thecrustal deviatoric stress associated with the horizontal
contraction of a homogeneous, isotropic, elastic earthmodel by self-gravitation, and concluded that the change
in gravitational potential energy due to dip-slip faulting
is several orders of magnitude greater than the observed
seismic energy. If this is true, any energy balance equation for non-gravitating earth models will lose its meaning. In reality, because of the rheological property of theearth’s mantle and the steady seafloor spreading and oceanic plate subduction in long time scale, the stress field
caused by self-gravitation will disappear except for its
isotropic part. So, the gravity-origin stress does not affectthe energetics in earthquake faulting substantially. To
sum up, we may use the well-known formula of the staticchange in elastic strain energy UE due to earthquake
faulting for non-gravitating elastic earth models (e.g., Aki
and Richards 1980; Rivera and Kanamori 2005),

�UE = − 21

�

final
dS,
(σijinitial + σijfinal )n−
j Di

(1)

in a good approximation. Here, σijinitial represents the tectonic-origin initial stress, σijfinal the final stress (the sum of
the initial stress and coseismic stress changes), and Difinal
the final slip on a fault  with its unit normal n−j .
When the faulting is aseismic, we can regard Eq. (1) as
the energy balance equation because there is no seismic
wave radiation. However, the occurrence of earthquakes
is normally accompanied by the radiation of seismic
waves. Kostrov (1974) and also Rivera and Kanamori(2005) evaluated the radiated seismic energy ER in an
ingenious but a bit troublesome way to avoid the problem
of gravitational potential energy change as

final
1
ER = 2dS(σijinitial + σijfinal )n−j Di� tfinal (2)σij n−dS,
dt− γeff S −Ḋijtinitial

�(t)where tinitial and tfinal represent the rupture initiation andtermination times, respectively. It should be noted herethat we rewrote the Kostrov’s original expression (2.26)
by using the formula of integration by parts. The second
term γeff S on the right-hand side of the above equation
represents the total fracture surface energy. Historically,
the fracture surface energy γeff was introduced to handle the effects of artificial stress singularity at crack tips
in the classical theory of rupture propagation based on
the Griffith-type fracture criterion. In reality, the leadingedge of crack is not a singular point but a finite process

Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

zone, where the shear strength (stress) decreases from apeak value to a residual value with the increase of fault
slip. Then, as discussed later in detail, we can incorporate
the second term of fracture surface energy into the third
term of frictionally dissipated energy on the fault.
Anyway, with Eq. (2) as a starting point, Rivera andKanamori (2005) and Kanamori and Rivera (2006)
derived a simplified expression, Dfinalτ (D)dD,
ER S = 21 (τinitial + τfinal )Dfinal −0(3)where τ is the average shear stress component parallel to the direction of slip D on the fault  with its area
S . The meaning of Eq. (3) is obvious; the excess of thereleased elastic strain energy over the work done forshear faulting gives the radiated seismic energy. Tracking back the history, this type of energy balance equation has been widely used for the evaluation of radiatedseismic energy since Kanamori (1977), but something iswrong. Rudnicki and Freund (1981) have directly evaluated the seismic energy radiated from a point dislocation source as  tfinal

1
5
M̈02 (t)dt. (4)
ER =
2
+
3(V
V
)
P
S
30πρVP5
tinitial

Here, M̈0 represents the second-order time derivativeof the seismic moment defined by M0 = µDS , and ρ , µ,
VP, and VS are the density, rigidity, P-wave velocity, andS-wave velocity of the elastic medium, respectively. This
result clearly shows the dependence of the radiated seismic energy on the rapture growth rate, but the simplified
Eq. (3) appears to contain no information about rupture
growth rates.
Incidentally, Kostrov (1974), Rudnicki and Freund(1981), and Rivera and Kanamori (2005) have made theattempt to strictly evaluate the second and third terms of
Eq. (2), but their attempts were not completed because of
the singularity of stress and particle velocity at the edgeof a rupture zone. To resolve this problem, as discussed
later in detail, we need to consider the physical process of
rupture growth governed by a slip-weakening fault constitutive law.

2 Basic concepts
First of all, we briefly summarize the basic conceptsof elastic potential energy (elastic strain energy) andmoment tensor. Then, we discuss the representation ofearthquake sources and the origin of background stressfields. By the way, throughout this article, Einstein summation convention is used to make mathematical expression concise.

Page 3 of 19

2.1 Decomposition of elastic potential energy

In a linearized theory, we can write the potential energy
density UE of elastic forces in a quadratic form of theelastic components εije of strain tensor εij ; that is,e.
UE = 21 Cijkl εije εkl(5)Here, the fourth-order tensor Cijkl denotes the elastic stiffness, which has the basic properties ofCijkl = Cjikl = Cijlk . Then, by definition, the stress tensor
components σij are given by the partial derivative of UEwith respect to εije ase
σij = Cijkl εkl(6)under the condition of Cijkl = Cklij . Using this expression,
we can rewrite Eq. (5) as

UE = 12 σij εije .

(7)

When the elastic material is isotropic, as demonstrated by Walpole (1984), the constitutive Eq. (6) canbe decomposed into two independent parts as.
eand τij = 2µsije
σkk = 3κεkk(8)where τij and sije denote the deviatoric stress and strain
tensors defined byeδij .
τij = σij − 13 σkk δij and sije = εije − 13 εkk(9)This means that the elastic deformation of isotropic
material can be decomposed into two different types; thatis, volumetric deformation and shearing deformation. As
demonstrated by Matsu’ura and Terakawa (2021), substitution of Eq. (9) into Eq. (7) yields

UE = 12 σv εve + 12 σs εse

(10)

with

εve = I1 (εije ), εse = 2andσv = 13 I1 (σij ), σs =J2 (sije )(11)J2 (τij ).

(12)Here, I1 (εije ) and I1 (σij ) denote the first invariant of thestrain tensor εije and the stress tensor σij , respectively,
and J2 (sije ) and J2 (τij ) denote the second invariant of thedeviatoric strain tensor sije and the deviatoric stress tensor τij , respectively. From Eq. (10), we can see that theelastic potential energy density of isotropic material canbe decomposed into two independent parts; that is, thevolumetric part and the shearing part.

Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Page 4 of 19

Using the constitutive Eq. (8) of isotropic material,
we can represent the elastic potential energy density in
Eq. (10) with only stress tensor components as

UE (σij ) = (18κ)−1 I12 (σij ) + (2µ)−1 J2 (τij ) 22= (18κ)−1 σkk+ (4µ)−1 τij 2 ,(13)

where � · �2 denotes the Frobenius norm of a secondorder tensor. One of the advantages of this representationis that we need not define the reference state to measure
the strain εij . Another advantage is that the change inelastic potential energy density is also written as the sumof its volumetric part UEv and shearing part (Matsu’ura
et al. 2019); that is, denoting the initial stress state and0 δ + τ0
the static stress change from it by σij0 (= 13 σkk
ijij ) and1�σij (= 3 �σkk δij + �τij ), respectively,�UE ≡ UE (σij0 + �σij ) − UE (σij0 ) = �U vE + �U sE

(14)

with
2
0
�σkk + (18κ)−1 �σkk
�UEv = (9κ)−1 σkk0+ 21 �σkk �σkk ,= (9κ)−1 σkk2�UEs = (2µ)−1 τij0 �τij + (4µ)−1 �τij 2= (2µ)−1 τij0 + 12 �τij �τij .

(15)

(16)

The important thing is that there is no interaction
between the volumetric and shearing parts, and eachof them consists of the interaction term of the stresschanges and the initial stress field and the self-interaction
term of the stress changes.
The elastic property of the actual earth is not necessarily isotropic even in a macroscopic scale. So, the logical consequences derived in this section may not beapplicable to the actual problems in the strict sense.
However, as demonstrated by Moakher (2008), we candefine the closest (in the sense of Euclidean distance)
IE to a given anisotropic
isotropic fourth-order tensor Cijkl
IE is the same asA
stiffness tensor Cijkl . The structure of Cijkl
that in the isotropic case. The effects of anisotropy arereflected only in the values of two material constants, κand µ. So, the logical consequences derived in this section are still applicable to the actual earth.
2.2 Physical meaning of moment tensor

Most theoretical problems in seismology are treatedin the framework of linear elasticity for mathematical
simplicity. The actual earth is not a linear elastic body.
In the mid-1970s, to bridge the gap between the actual

Fig. 1 Geometry and notation of an elastic body V bounded
by a traction-free external surface Sex and an internal closed surface
Sin. The unit normal vectors n are taken to be outward to the surfaces.
n+ and n− indicate the unit normal vectors to the front side  +and back side  − of a fault , respectively

earth and the linear elastic model, Backus and Mulcahy(1976) introduced the concept of moment tensor, which
includes the seismic moment tensor (Kostrov 1974) asa special case. The following is a physical interpretationof the Backus–Mulcahy moment tensor based on thetheoretical consideration about the equation of motion
in continuum mechanics (Matsu’ura et al. 2019).
We consider a solid body occupying the region V
bounded by a traction-free surface Sex with its outer
unit normal vector n = (n1 , n2 , n3 ) in a Cartesian coordinate system (x1 , x2 , x3 ) as shown in Fig. 1. In a linearized theory, irrespective of the stress–strain relationof material, we can write the equation of motion of thesolid body for t ≥ 0 as

ρ üi (x, t) = ∂j σij (x, t) + fi (x).

(17)

Here, üi represents the second-order time derivative of
the i -component of a displacement vector ui , σij the ij-component of a stress tensor, and fi the i -component ofa body force vector without time variation. The symbol∂j denotes partial differentiation with respect to the spatial coordinate xj . Let’s suppose that the solid body wasin an equilibrium state, 0 = ∂j σij0 (x) + fi (x), for t < 0.
Then, in order to solve the equation of motion, we needto specify the relationship between the stress change �σij(= σij − σij0) and the displacement ui through the strainchange �εij defined geometrically as�εij = 12 (∂j ui + ∂i uj ).

(18)

The occurrence of inelastic deformation, such as brittle
fracture and/or plastic flow, at the defects distributed
within the earth brings about elastic deformation in theregion surrounding the defects. Thus, the strain �εij
defined by Eq. (18) is not necessarily elastic anytime andMatsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Page 5 of 19

anywhere; that is, denoting the inelastic part of the strain�εij by �εija,�εij = �εije + �εija .

(19)Then, from Eqs. (6) and (19), the stress–strain relation in
the actual earth can be written asa).
�σij = Cijkl (�εkl − �εkl

(20)

tensor density mij (x) defined by Eq. (24) is not physicalsubstance. As can be seen from Eq. (25), the moment
tensor density appears as a virtual source term in thetheory of linear elasticity. If the distribution of inelasticstrain changes is restricted within a region Vs (Fig. 1),
by definition, the moment tensor Mij is given by thea (x) over V asvolume integral of Cijkl �εklsa(x)dV .
Mij ≡mij (x)dV =Cijkl �εkl
(27)
Vs

We now consider the occurrence of some inelastic
event at t = 0 within the solid body in static equilibrium. Then, denoting the changes in stress and strain
due to the inelastic events by �σij (x, t) and �εij (x, t),respectively, we may rewrite the equation of motion(17) asρ üi (x, t) = ∂j �σij (x, t),(21)

and the stress–strain relation (20) as

�σij (x, t) = �σijm (x, t) − mij (x, t)(22)

�σijm (x, t) = Cijkl �εkl (x, t),(23)a(x, t).
mij (x, t) = Cijkl �εkl

(24)

with

Here, �σijm is a model stress, calculated from the displacements ui through the geometrical strain under
the assumption of linear elasticity, and mij is called themoment tensor density, which is equivalent to the stressglut tensor, Ŵij ≡ �σijm − �σij , defined by Backus andMulcahy (1976). Substituting Eq. (22) in Eq. (21), we
obtain the equation of motion to be solved in the framework of linear elasticity asρ üi (x, t) = ∂j �σijm (x, t) − ∂j mij (x, t)(25)

ρ üi (x, t) = ∂j �σijm (x, t) + fim (x, t)

(26)

or

with fim (x, t) = −∂j mij (x, t), which is called the bodyforce equivalent (Burridge and Knopoff 1964). The inelastic strain referred to here is almost the same as the
eigenstrain in Mura (1982), but slightly different fromthe stress-free transformation strain of an inclusion in
Eshelby (1957). The Mura’s eigenstrain is a virtual one,and so it can coexist with elastic strain, while the Eshelby’s inclusion and the elastic medium are mutually exclusive because both of them are substantial.
It should be noted here that the inelastic strain
a (x)
�εkl
itself is physical substance, but the momentVs

2.3 Representation of earthquake sources

In the mid-1960s, it was convinced that the source ofearthquakes is the development of tangential displacement discontinuity across a fault (Burridge and Knopoff 1964). We now consider displacement discontinuity
across an interface �(η) in an isotropic linear elastic body
V with the bulk modulus κ and the rigidity µ. In this case,the moment tensor Mij defined by Eq. (27) can be written
as
 

a
(x)δij dV .
Mij =2µ�sija (x) + κ�εkk(28)
Vs

In general, displacement discontinuity at the interface
�(η) is represented asui (x) = �ui (η)H(x − η),(29)−
where �ui (η) ≡ u+i (η) − ui (η) is the i-component of
displacement jump across �(η), and H(x − η) denotes
the Heaviside step function. Then, denoting the Diracdelta function by δ(x − η) and a unit normal vector
pointing from the backside of �(η) to the front side byn− (η) (Fig. 1), we can represent the distribution of inelastic strain �εija (x) localized at �(η) as−�εija (x) = 12 [�ui (η)n−j (η) + �uj (η)ni (η)]δ(x − η).
(30)
The isotropic part of the above inelastic strain isa(x)δij = �uk (η)n−�εkkk (η)δij δ(x − η),and so the deviatoric part can be written as−�sija (x) = 12 [�ui (η)n−j (η) + �uj (η)ni (η)]
− 31 �uk (η)n−k (η)δij δ(x − η).

(31)

(32)

Substitution of the expressions (31) and (32) into
Eq. (28) yields the moment-tensor representation of
displacement discontinuity across the fault (Matsu’ura
et al. 2019):

Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Page 6 of 19

2.4 Origin of background stress fields

Fig. 2 Geometry and notation of the displacement discontinuity
u across a fault . n− is the unit normal vector to the backside  −of the fault, and D represents the fault-parallel component of u

Mij =



�

m�
ij (η)dS

(33)

with
−
−
m�
ij (η) =µ[�ui (η)nj (η) + �uj (η)ni (η)]+ (κ − 23 µ)�uk (η)n−k (η)δij .

(34)

Here, m�ij (η) represents the area density of moment tensor on �(η).
In the case of shear faulting, denoting the tangential component of displacement discontinuity−
�ui − (�uk n−k )ni by Di or Dνi (Fig. 2), we can reduce
Eq. (33) to−Mij = µD(η)[νi (η)n−j (η) + νj (η)ni (η)]dS, (35)�

which is nothing but the Kostrov’s seismic moment tensor. Since the slip direction vector ν is perpendicular to
the normal vector n− (νk n−
k = 0), the trace Mkk of theseismic moment tensor is always zero. On the other
hand, in the case of crack opening, denoting the normal−component of displacement discontinuity (�uk n−k )ni by−Cni , we can reduce Eq. (33) to−−−Mij =µC(η)[n−i (η)nj (η) + nj (η)ni (η)]dS�+ (κ − 32 µ)C(η)δij dS,
�

(36)

which has a nonzero isotropic part 13 Mkk = κ � C(η)dS .
This difference between shear faulting and crack opening
is crucial when we consider the effects of gravity force on
the energy balance of the whole earth.

As demonstrated in Sect. 2.1, the elastic potential energy of
an isotropic body and also its change can be decomposed
into two independent parts, namely the volumetric partand the shearing part. This implies that we can treat theisotropic part and deviatoric part of stress fields separately.
The origin of the isotropic part of the earth’s stress (lithostatic pressure) is of course self-gravitation. The questionis the origin of the deviatoric part (tectonic stress). Our
understanding is that the origin of tectonic stress is not thegravitational contraction of the elastic earth but the stableprocess of converting thermal energy into kinetic energy inthe earth’s interior (mantle convection), which causes the
motion of tectonic plates. In this section, with the idea of
the interaction energy between two stress systems (Eshelby
1956), we conceptually explain the origin of backgroundtectonic stress. The derivation of interaction energy shown
below is based on Sect. 13 of Mura (1982).
We consider a non-gravitating elastic body V bounded
by a traction-free surface Sex (Fig. 1). In the initial state, theelastic body is assumed to be unstrained, and so its elastic
potential energy (elastic strain energy) is also zero. Let us
suppose some process with inelastic strain distribution εija
Vs (⊂ V ). Then, the elastic potenoccurred within a regiontial energy UE ≡ V UE dV produced by this process isevaluated as


e
1
1
UE = 2σij εij dV = 2
σij (εij − εija )dVVVσij ∂j ui dV − 21σij εija dV
= 12
(37)
V
Vs

= − 12σij εija dV ,Vsbecause σij εij = 21 σij (∂j ui + ∂i uj ) = σij ∂j ui and
σij ∂j ui dV =σij ui nj dV −∂j σij ui dV = 0.
V

Sex

V

(38)
Note that σij nj = 0 on Sex (the boundary conditionof traction free) and ∂j σij = 0 in V (the equation ofequilibrium).
When the source region Vs consists of two independentparts, Vs0 and Vs1 (Fig. 3), applying the formula (37), we canwrite the produced elastic potential energy asUE = − 12(σij0 + σij1 )εija0 dV − 21
(σij0 + σij1 )εija1 dV
Vs0
Vs1



0 a0
1 a1
1
1
σij εij dV − 2σij εij dV −σij0 εija1 dV ,
= −2
Vs0

Vs1

Vs1

(39)

Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Page 7 of 19

Here, UE represents the change in elastic potential
energy due to the sudden faulting. In the case of tangential displacement discontinuity (fault slip) at the interface
�(η), as shown in Sect. 2.3, the corresponding inelasticstrain �εija can be written as

Fig. 3 Two independent source regions Vs0 and Vs1in an elastic body V
bounded by a traction-free external surface Sex. σij represents a common
stress field

because

but Eq. (1).

σij1 εija0 dV =Vs1

σij0 εija1 dV .

(40)

The proof of Eq. (40) is given as follows:
σij1 εije0 dV ,σij1 (εij0 − εije0 )dV = −σij1 εija0 dV =
Vs0

Vs1

σij0 εija1 dV =

V

V





σij0 εije1 dV ,



e1
dV .
σkl0 εkl

σij0 (εij1 − εije1 )dV = −

V

V

(41)

(42)

and



V

�

with σijinitial = σij0 and σijfinal = σij0 + �σij , which is nothing

Vs0



−
�εija (x) = 12 [Di (η)n−j (η) + Dj (η)nj (η)]δ(x − η).
(46)
Substituting this expression in Eq. (44), we obtain
−�UE = −
(σij0 + 12 �σij ) 21 (Di n−
j + Dj nj )dS
�

(47)
1
(σijinitial + σijfinal )n−= −2
j Di dS

σij1 εije0 dV =



V

e1 e0
εij dV =Cijkl εklV

(43)

Following Andrews (1978), we regard εija0 as the inelastic strain distribution produced by the long-term
stable process of converting thermal energy into kinetic
energy in the earth’s interior. Then, Eq. (39) can be readas the elastic potential energy after the occurrence ofa sudden faulting with inelastic strain distribution εija1in the pre-stressed state of σij0 . Therefore, changing thenotations of the inelastic strain εija1, the stress change
σij1, and the source region Vs1 associated with the sudden faulting to �εija , �σij , and Vs , respectively, we canrewrite Eq. (39) as0�UE = UE − UE = −(σij0 + 12 �σij )�εija dV (44)
Vs

with

3 Energetics in earthquake faulting
As a starting point of discussion, we directly derivea basic equation of mechanical energy balance indynamic shear faulting from the equation of motion
for an elastic body subjected to gravity-origin isotropic
stress and tectonic-origin deviatoric stress. Then, we
review the indirect evaluation of radiated seismicenergy based on a simplified energy balance equation and the direct evaluation based on the analyticalsolution of displacement fields for a point dislocation
source, and call attention to inconsistency between
them. Finally, to resolve this inconsistency, we consider
the physical process of rupture growth governed by aslip-weakening fault constitutive law from the viewpoint of energy balance.
3.1 Mechanical energy balance in earthquake faulting

Based on the discussion in Sect. 2.4, we consider an elastic body subjected to gravity-origin isotropic stress andtectonic-origin deviatoric stress. So, the strain εij appearing in this section is elastic. The elastic body occupies theregion V bounded by a traction-free external surface Sex
and an internal surface Sin enclosing a fault  as shown in
Fig. 1, where the direction of unit normal vectors n to theexternal and internal surfaces is always taken outward.
As for the fault , in accordance with the convention, we
take a unit normal vector n− pointing from its backside
 − to front side  +. Although the equation of motion
in linear elasticity has been already given in Eq. (17), we
rewrite it for convenience of the later discussion asg

UE0 = − 12Vs0

σij0 εija0 dV .

(45)

ρ∂t u̇i (x, t) = ∂j σij (x, t) + fi (x) for t ≥ 0,(48)

Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Page 8 of 19gwhere u̇i is the i -component of a velocity vector, fi is thei -component of a gravity force vector, and the symbol ∂trepresents partial differentiation with respect to tine t . It
is supposed here that the elastic body was in an equilibrium state for t < 0 under the background tectonic stressg
σij0 (x) defined in Sect. 2.4 and the gravity force fi (x). So,
the stress appearing in Eq. (48) should be read as the sumof the background tectonic stress and the stress change
measured from it; that is, σij (x, t) = σij0 (x) + �σij (x, t).
Now, following Sect. 7.8 of Leigh (1968), we multiply
both sides of Eq. (48) by u̇i:g(49)

ρ u̇i ∂t u̇i = u̇i ∂j σij + u̇i fi .

Here, the first term on the right-hand side can be rewritten asu̇i ∂j σij = ∂j (u̇i σij ) − (∂j u̇i )σij = ∂j (u̇i σij ) − ε̇ij σij(50)
Then, substituting the above
expression into Eq. (49) and integrating each term over
the whole elastic region V , we obtain the following
energy balance equation:with ε̇ij = 21 (∂j u̇i + ∂i u̇j ).

V12 ρ∂t (u̇i u̇i )dV =V∂j (u̇i σij )dV −Vε̇ij σij dV +Vgu̇i fi dV .

(51)
Furthermore, applying the Gauss divergence theorem tothe first term on the right-hand side and replacing thestress tensor σij in the second term with Cijkl εkl , we canrewrite Eq. (51) as


1
ρ∂
(
u̇
u̇
)dV
=
u̇i σij nj dS
t i i
2
V
Sin

1
−
(52)
2 Cijkl ∂t (εij εkl )dVVgu̇i fi dV .
+V
Here, we dropped the term Sex u̇i σij nj dS because thetraction on the external surface Sex is always zero. Notethat the left-hand side of Eq. (52) represents the increase
rate of kinetic energy, the first term on the right-handside the rate of work done by the traction acting on theinternal surface Sin , and the second and third terms therates of change in elastic and gravitational potential
energy, respectively.
In the case of shear faulting, shrinking the internalsurface Sin to a fault  with its front side  + and backside  − and denoting the i-component of the tangen−tial velocity jump u̇i ≡ u̇+i − u̇i across  by Ḋi , we canrewrite the first term on the right-hand side of Eq. (52) asSinu̇i σij nj dS =

� + +� −u̇i σij nj dS = −�Ḋi σij n−
j dS.

(53)
As for the last gravity force term, applying the idea ofinteraction energy in Sect. 2.4, Savage and Walsh (1978)
have demonstrated that the volume integral over thewhole elastic region V could be transformed into the surface integral over the fault  as follows:



g
g
g
u̇i fi dV = −u̇i ∂j σij dV =Ḋi σij n−
j dS, (54)
V

V

�

g
where σij represents the internal stress field that satisfiesgg
the equation of equilibrium ∂j σij + fi = 0. Substitutionof Eqs. (53) and (54) into Eq. (52) yields
1∂t
σij n−
j Ḋi dS2 ρ u̇i u̇i dV = −
V
�

1
− ∂t
2 σij εij dVVg+
σij n−
j Ḋi dS.

(55)

�

When the internal stress field caused by the gravitygg
force is isotropic ( σij = 13 σkk δij ), the last term on theright-hand side of Eq. (55) vanishes:



g
g
g
−
1
1
σij n−
dS
=
Ḋ
σ
δ
n
dS
=
σkk n−
Ḋ
i
ij
i
j
j
i Ḋi dS = 0,
kk
3
3
�

�

�

(56)
=
0
for
shear
faulting.
Then,
discarding
thebecause Di n−ilast term associated with gravity force, we obtain a basic
equation of mechanical energy balance in dynamic shear
faulting as
K̇ (t) = −Ẇ (t) − U̇E (t)with the increase rate of kinetic energy1K̇ (t) = ∂t
2 ρ u̇i (x, t)u̇i (x, t)dV ,V

the rate of work done for shear faultingẆ (t) =σij (η, t)n−j (η)Ḋi (η, t)dS,�and the rate of change in elastic potential energy1U̇E (t) = ∂t
2 σij (x, t)εij (x, t)dV .
V

(57)

(58)

(59)

(60)

We consider a shear faulting that initiated at
t = tinitial on a fault  and terminated at t = tfinal .
Hereafter, tinitial and tfinal are abbreviated to ti and tf ,

Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Page 9 of 19

respectively, unless misleading. Let us suppose that the
disturbances caused by the shear faulting died down
by t = tstationary ( ts in abbreviation) except the seismic
waves traveling outward from the source. Then, integrating each term of Eq. (57) with respect to t from ti
to ts (>> tf ), we obtain(61)

K (ts ) = −W (ts ) − �UE (ts )with

K (ts ) ≡ ts

W (ts ) ≡

ti

 ts
ti



u̇i (x, ts )u̇i (x, ts )dV ,

σij (η, t)n−j (η)Ḋi (η, t)dS,

K̇ (t)dt = 12 ρẆ (t)dt =

 tf

dt

ti

V (ts )

�(t)

(62)

(63)

and, from Eqs. (44) and (47),
 ts
U̇E (t)dt
�UE (ts ) ≡ti
[σij (η, ti ) + σij (η, tf )]n−
= − 21j (η)Di (η, tf )dS,�

(64)
where  represents the final rupture area �(tf ). In thiscontext, we can read K (ts ) as the radiated seismic energy
ER , and so Eq. (61) is written as
ER = 12



�

final
dS −
(σijinitial + σijfinal )n−j Di

 tfinal
tinitial

dt



�(t)

σij n−
j Ḋi dS.

(65)
This expression is essentially the same as Eq. (2), if we
can incorporate the fracture surface energy term γeff S in
Eq. (2) into the last integral term of frictionally dissipatedenergy (Fukuyama 2005).

ER = 12



−

�



[τi (η) + τf (η)]Df (η)dS − γeff S�dS Df (η)(66)

τ (D|η )dD.

0Here, D denotes the vector magnitude of fault-slip, andτ = σij n−
j νi represents the scalar projection of the trac
tion vector Ti = σij n−
j in the slip direction νi = Di D
(Fig. 4). In this expression, the shear stress τ is assumed
to be a single-valued function of the fault slip D; thatis τ = τ (D). Then, we can replace the integration
 tf
 Df
ti τ (t)∂t D(t)dt with the integration 0 τ (D)dD . In
the natural process of shear faulting, the shear stress τ
defined above is always positive.
Kanamori and Rivera (2006) further simplified theenergy balance equation by neglecting the fracture surface energy term γeff S and taking the spatial average ofthe shear stress τ (η) and the fault slip D(η) over the fault
 with its final dimension S as
 Dfinalτ (D)dD,
ER S = 21 (τinitial + τfinal )Dfinal −0(67)which is nothing but Eq. (3). This energy balance equation is often interpreted by using the τ −D diagram in
Fig. 5 as follows: (1) the trapezoidal area under the thickbroken line, which corresponds to the first term on theright-hand side, represents the elastic potential energy
released by shear faulting, (2) the shaded area under the
thick solid curve (rupture growth path), which corresponds to the second integral term, represents the work
done for shear faulting, (3) the most part of the released
elastic potential energy is spent for the shear faulting, and4) the remaining part (the difference between the areas

3.2 Radiated seismic energy

Nowadays, using an inversion technique, we canestimate the dynamic process of earthquake faulting
from seismological and/or geodetic data. A direct
way to evaluate the radiated seismic energy of anearthquake is of course to perform the integrationin Eq. (62). Another indirect way is to evaluate theright-hand side of Eq. (65), which relates the radiatedseismic energy and the history of fault slip during theearthquake.
3.2.1 Indirect evaluation of radiated seismic energy

Rivera and Kanamori (2005) simplified the energy balance Eq. (2) by exchanging the order of integration of thelast frictionally dissipated energy term as

Fig. 4 Geometry and notation of a traction vector T on a fault . n−is the unit normal vector to the backside  − of the fault, and ν represents
the unit slip-direction vector

Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Page 10 of 19

where r =Fig. 5 A schematic diagram for interpreting the energy balancein earthquake faulting. The thick solid curve τ = τ (D) represents
a rupture growth path; τi, τp, and τf are the initial, peak, and final valuesof the shear stress τ , respectively. Further explanation is given in the text

marked with + and -) gives the seismic wave energy ER
radiated from the source. Such interpretation seems tobe plausible; the most primitive one is found in Kanamori(1977), which gives a theoretical basis of the moment
magnitude scale Mw of seismic events. But something iswrong. For instance, when the growth of shear rupture isquasi-static, the left-hand side of Eq. (67) must become
zero because of no seismic wave radiation, but the righthand side will not change anything because it appears to
be independent of rupture growth rates.
3.2.2 Seismic waves and permanent deformation

The displacement field due to a moment tensor Mpq
located at the origin ( x = 0) of a Cartesian coordinate
system is generally expressed as tui (x, t) =∂q Gip (x, t − t ′ ; 0, 0)Mpq (t ′ )dt ′ .
(68)0Here, Gip (x, t; 0, t ′ ) is the solution of the equation of
motion to a unit impulsive force in the xp-direction,
δip δ(x)δ(t − t ′ ), applied to a point x = 0 at t = t ′ , and∂q Gip (x, t; 0, t ′ ) represents the partial derivative ofGip (x, t; 0, t ′ ) with respect to the spatial coordinate xq .
The explicit form of Gip for a homogeneous, isotropic,
unbounded elastic medium, which is called Stokes’ solution, is given in Sect. 4.2 of Aki and Richards (1980) as
1 r/VS ′t δ(t − t ′ )dt ′4πρGip (x, t; 0, 0) = (3γi γp − δip ) 3
r r/VP1r+ γi γp 2 δ t −
VP
VP r


1
r
,
− (γi γp − δip ) 2 δ t −
VS
VS r
(69)



x12 + x22 + x32 and γi = xi r (= ∂i r ) denotethe distance and direction cosines from the source point
0 to the receiver point x.
In the case of shear faulting, substituting
the corresponding moment-tensor representation,−
Mpq = M0 (t)(νp n−
q + νq np ), into Eq. (68), we obtain the
explicit expression of displacement field due to a seismic
moment tensor as r/VS6(5abγi − aνi − bn−i )ui (x, t) =t ′ M0 (t − t ′ )dt ′4πρr 4r/VP
2(6abγi − aνi − bn−
r
i )
+
t
−
M
0
VP
4πρVP2 r 2
−
3(4abγi − aνi − bni )
r
t
−
M
−
0
VS
4πρVS2 r 2


2abγi
r
+
Ṁ0 t −
VP
4πρVP3 r(2abγi − aνi − bn−
r
i )
−
,
Ṁ0 t −
VS
4πρVS3 r(70)−
where a = γ · n , b = γ · ν, and M0 (t) is a cumulative
seismic moment function.
From the difference in geometrical attenuation, the firstterm on the right-hand side of Eq. (70) is traditionallycalled the near-field term, the second and third terms theintermediate-field terms, and the fourth and fifth terms
the far-field terms. Specifically, the first three terms decay
in amplitude with the square of the source–receiver distance r , while the last two terms decay with the distancer . Since the cumulative seismic moment function M0 (t)monotonically increases with time t from 0 at t = ti up
to M0 at t = tf , its time derivative Ṁ0 (t − r/c) vanishes
for t − r/c > tf . So, from Eq. (70), we can see that the
far-field terms become zero everywhere behind the wave
packets traveling outward at the phase velocity c (VP for
P-wave and VS for S-wave). On the other hand, the nearand intermediate-field terms remain as permanent deformation after disturbance died down (t − r/c > tf ), whichcauses the change in elastic potential energy.
During the rupture growth (ti ≤ t ≤ tf ), the far-fieldterms and the near- and intermediate-field terms interact
with each other, but such interaction will die down soon.
So, we may define the cumulative kinetic energy K (t), work
done for shear faulting W (t), and elastic potential energychange �UE (t) for a sufficiently large time t (≥ ts) as in
Eqs. (62), (63), and (64). As discussed in Sect. 3.3, the lastMatsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Page 11 of 19u̇Sθ (x, t) =1r) eθ ,
cos 2θ cos φ M̈0 (t −3VS4πρVS r(71b)

1
r
) eφ ,cos θ sin φ M̈0 (t −3V4πρVS rS
(71c)where er , eθ and eφ are the mutually orthogonal unit vectors that define a local coordinate system at a receiverpoint (Fig. 6).
The spherical fronts of the P- and S-waves expand outward with time at the phase velocity c (Fig. 6). So, for a
sufficiently large time ts (>> tf ), the cumulative kinetic
energy K (ts ) is given by the volume integral of the kinetic
energy density 12 ρ|u̇|2 over the sphere with the radiusr = c(ts − ti ): c(ts −ti )  2π
 π
2 2
1
K (ts ) =drdφ
2 ρ|u̇| r sin θdθu̇Sφ (x, t) = −0

Fig. 6 A Cartesian coordinate system fixed to a seismic source. The
x1-axis is taken in the direction of fault slip ν and the x3-axis is takenin the direction of fault normal n−. The mutually orthogonal unit
vectors, er , eθ and eφ , define the local coordinate system at a receiverpoint (r, θ , φ). The radial distance range of the wave packet travelingoutward at a phase velocity c is indicated by cT

two cumulative functions, W (t) and �UE (t), are mutuallydependent because the growth of fault rupture is controlled
by the influx of released potential energy.
3.2.3 Direct evaluation of radiated seismic energy

Now, we directly evaluate the total energy radiated as seismic waves from the explicit expressions of far-field terms
in Eq. (70). For this purpose, we take a new Cartesian coordinate system fixed to a seismic source as shown in Fig. 6,
where the x1-axis is taken in the direction of fault slip (thenaν = sin θ cos φ) and the x3-axis in the direction of fault
normal (then an = cos θ). From Eq. (70), we can see thatthe particle motion of the far-field P-wave (the fourth term)
is perpendicular to the wave front, while that of the far-fieldS-wave (the fifth term) is tangential. Furthermore, the tangential particle motion of S-wave can be decomposed into
two mutually orthogonal components. Following Sect. 4.3
of Aki and Richards (1980), we represent the mutuallyorthogonal particle velocities asu̇P (x, t) =1r
) er ,
sin 2θ cos φ M̈0 (t −3V
4πρVP r
P

(71a)

0

0

(72)

with

2
 2  2 

 
 

|u̇|2 = u̇P  + u̇Sθ  + u̇Sφ  .

(73)

Substituting the explicit expressions of the mutuallyorthogonal particle velocities given by Eq. (71) into
Eq. (72) and reducing the interval of integral for r toc(ts − tf ) ≤ r ≤ c(ts − ti ), we can evaluate the contribution to the total kinetic energy of each particle motion as
follows:
 VP (ts −ti )1rP2
M̈0 t −
K (ts ) =dr
VP
30πρVP6 VP (ts −tf )
(74a)
 �T
1
2 ′
′
=
M̈
(t
)dt
,
0
30πρVP5 0

K Sθ (ts ) ==K Sφ (ts ) =7
15 × 16πρVS6715 × 16πρVP51
3 × 16πρVS6 VS (ts −ti )rM̈02 t −drVS
VS (ts −tf ) �TM̈02 (t ′ )dt ′ ,0 VS (ts −ti )

(74b)

r
t−
dr
VS



M̈02
VS (ts −tf ) �T1
M̈02 (t ′ )dt ′ ,=3 × 16πρVP5 0(74c)

Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

where T = tf − ti denotes the rupture duration. Finally,summing up the above expressions, we obtain  �T
1
5
K (ts ) =M̈02 (t ′ )dt ′ .
2
+
3(V
V
)
P
S
60πρVP5
0(75)
Rudnicki and Freund (1981) has obtained the equivalent
result in a slightly different way, though their radiatedenergy in Eq. (4) is just twice as large as the present evaluation, because they adopted the definition of radiatedenergy by Haskell (1964), which is essentiallythe same asK (ts ) in Eq. (72) except for the factor 1 2. The important
thing is that the origin of the cumulative kinetic energy,
which can be read as the radiated seismic energy, is not
the moment rate Ṁ0 (t) but the moment rate changeM̈0 (t).
The seismic moment function M0 (t) is defined asM0 (t) ≡ µD(η, t)dS.
(76)
�(t)

As pointed out by Ampuero et al. (2006), the fault slip
distribution of actual earthquakes is heterogeneous inboth time and space because of the spatiotemporal heterogeneity of stress drop. So, the moment rate functionṀ0 (t) of actual earthquakes is not always unimodal.
Nevertheless, in the following, we consider the casewhere the seismic moment monotonically increases
with time for sake of simplicity. First, to compare thedirect evaluation of radiated seismic energy in Eq. (75)with the indirect evaluation in Eq. (67), we take the spatial average of the slip function D(η, t) over the fault with its final dimension S , and denote it by D(t); that is
M0 (t) = µSD(t). Then, to say it without worrying about
being misunderstood, the radiated seismic energy is proportional to the square of the fault slip acceleration D̈(t).

Page 12 of 19

When the fault area concentrically expands from 0 at
t = ti to S at t =√tf , wemay define the average rupture
velocity as vr = S/2π �T with T = tf − ti. Figure 7
schematically shows the rupture velocity-dependence ofthe fault slip D, slip velocity Ḋ, and slip acceleration D̈,
where √
the time
 t is normalized by the rupture duration
�T = S/2π vr. From this diagram, we can see that thefault slip acceleration is proportional to the square of theaverage rupture velocity vr. On the contrary, the simplified Eq. (67) appears to be independent of the rupture
velocity.
3.3 Energetics of rupture growth

To correctly evaluate the energy balance Eq. (2) or (65),
we need to make clear the physical mechanisms of rupture initiation, acceleration, deceleration, and termination. Until the 1970s, this problem has been treated with
the Griffith-type fracture criterion in the framework ofclack theory. One of the difficulties in the crack theory
is how to deal with the singularity of stress at crack tips,
though it does not exist in reality because the strengthof actual materials is finite. In the 1980s, the stick–slip
experiments of rock samples revealed the existence ofa breakdown zone just behind the propagating rupture
front, where the creation of new fracture surfaces is proceeding. We can use the shear stress–slip displacement
relation observed during the breakdown process, whichis called the fault constitutive relation, as a fundamental equation governing the rupture process instated of
the Griffith-type fracture criterion. Thus, in the 1990s, it
became possible to numerically simulate the entire process of earthquake generation by solving a coupled system of the equation of motion and the fault constitutive
relation under given initial and boundary conditions. In
this section, we discuss the physical process of rupture
growth governed by the slip-weakening fault constitutive
law and elucidate the true meaning of the energy balance
Eq. (2) or (65).
3.3.1 Griffith‑type fracture criterion

Fig. 7 A schematic diagram showing the rupture-velocity dependenceof the fault slip D, slip velocity Ḋ, and slip acceleration D̈. The verticalaxis t ′ represents
 the time normalized by the rupture duration√
�T = S/2π v r (S: the final fault area, v r: the average rupture velocity)

Early in the twentieth century, Griffith (1921) considered
the stability of a traction-free crack existing in the elastic
body subjected to an external surface traction and proposed an energetics-based global fracture criterion as follows. If the energy release rate Gs (the decrease in total
potential energy of the elastic body per unit crack extension) exceeds the fracture surface energy γeff (the energy
spent for newly creating unit crack area), the unstablerupture of crack begins. In the case of earthquake faulting, we need not consider the work done by the external surface traction, and so we can read the decrease inMatsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Page 13 of 19

Fig. 8 A snapshot of stress and slip distributions along a shear crack
extending in the x-direction at a rupture velocity vr . The red and bluelines represent the stress and fault slip distributions at a moment t ,respectively

total potential energy as the decrease in elastic strainenergy. In the 1950s, on the basis of the analysis of stressand strain fields near crack tips, Irwin (1957) related theenergy release rate Gs to the stress intensity factor Ks (thefactor indicating the intensity of stress singularity
at cracktips) as Gs = cKs2 with c = (3κ + 4µ) 4µ(3κ + µ) forMode I and II cracks and c = 1 2µ for Mode III crack.
Incidentally, this general relationship between Gs and Ks
can be obtained by substituting the asymptotic expressions of displacement and stress fields around a crack tip
into the formula evaluating potential energy changes due
to hypothetical crack growth (Rice and Drucker1967).
From the 1960s to the 1970s, the Griffith’s fracture
criterion for static cracks has been extended to the caseof dynamically propagating cracks by considering the
inflow rate of released elastic strain energy into the cracktip (e.g., Kostrov 1964, 1966, 1975; Freund 1972, 1979).
For the energetics of a propagating Mode III (anti-plane
shear) crack, first Kostrov (1966) and later Aki and Richards (1980) have obtained a complete solution. The derivation of the energy balance equation shown below isgiven in Sect. 15.2 of Aki and Richards (1980). We consider an anti-plane shear crack extending in the x-direction at an instantaneous rupture velocity vr. Figure 8 isa snapshot of stress and slip distributions along a crack
surface ( x < c ) and its extension ( x > c ) at a momentt (≥ 0). In this case, the energy balance at the crack tip
x = c(t) can be written as
1 − vr VS 21K
= γeff (vr < VS )(77)
2µ 1 + vr VS inst
withKinst =


�τ x′ , t − [c(t) − x′ ]/VS
2 c(t)√dx′ .
π c(t)−VS tc(t) − x′(78)

Fig. 9 Dependence of the energy distribution factor pon the instantaneous rupture velocity vr

Here, the left-hand side of Eq. (77) represents the inflow
rate of released elastic strain energy into the crack tip
propagating at vr, and the right-hand side is the fracturesurface energy at x = c(t). As can be seen from Eq. (78),
Kinst represents the cumulative effect of past stress drop�τ (x, t) on the present intensity of stress singularity at
x = c(t), which does not depend on the rupture velocity
vr. When the crack extension is quasi-static (vr = 0), theleft-hand side of Eq. (77) corresponds to the static energy
release rate
the energy distribution factor
 Gs, because
p(vr ) = (1 − vr VS ) (1 + vr VS ) becomes unity andthe instantaneous stress intensity factor Kinst is identicalto the static stress intensity factor Ks. As the rupture
velocity vr increases from 0 to VS, the factor p(vr )decreases monotonically from 1 to 0 (Fig. 9). This meansthat only a part ( p) of the supplied elastic strain energy22µ) is spent for creating new crack area, and the( Kinstremaining parts (1 − p) are spent for accelerating slip
motion, which causes seismic waves as shown in
Sect. 3.2.2. To satisfy the energy balance Eq. (77) for given
Kinst and γeff , the factor p(vr ) must take an appropriate
value at each moment. In other words, the rupture velocity vr at each moment is determined so as to satisfy
Eq. (77). Thus, we can say that the Griffith-type fracture
criterion is a fundamental equation governing the wholerupture process; a good numerical example is given in
Das and Aki (1977).
3.3.2 J‑integral and slip‑weakening fault constitutive law

In the classical crack theory, as shown √in Fig. 8, the shear
x − c(t)-type sinstress τ just outside a crack has a 1
gularity, but the fault slip D is zero everywhere. On theMatsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Page 14 of 19

P=

Fig. 10 A snapshot of stress and slip distributions
along a propagating shear crack with a finite process zone δc . In
the process zone, the shear stress τ degrees from a peak value τpto a residual value τr with the increase of fault slip DA

UE dS −

γeff = τx=c(t) × Dx=c(t) = ∞ × 0 (indeterminate).
(79)
This result is not irrational but physically meaningless.
In reality, the leading edge of crack would not be a singular point but a finite process zone, where the shear
stress (strength) decreases from a peak value τp to a residual value τr with the increase of fault slip D as shown in
Fig. 10.
For simplicity, we consider a two-dimensional antiplane shear crack  with a finite process zone δc as shownin Fig. 11. Here, the area A surrounded by the externalboundary (Ŵex ) and the upper and lower crack surfaces(+ and −) is elastic and subjected to some stress σij
and strain εij . When the residual shear stress on  is zero(τr = 0), the total potential energy of the elastic body perunit thickness is given by

(80)

Ŵ

(81)
we can evaluate the decrease of the total potential energy
P by quasi-static unit crack extension asδPδ�→0 δ�other hand, the fault slip just inside the crack takes a√c(t) − x-type distribution, but the shear stress is zero
everywhere. Then, regarding the fracture surface energy
γeff as the work done for the crack tip passing the singularpoint, we obtainTi ui dŴ,Ŵexwhere UE = 12 σij εij , Ti = σij nj , and ui and ni denote the i-components of a displacement vector and a unit normal
vector, respectively. Using a path independent integral
J (Ŵ) defined by Rice (1968),J (Ŵ) =(UE − σxz ∂x uz )nx − (σyz ∂x uz )ny dŴ,− limFig. 11 A two-dimensional anti-plane shear crack  with a finite
process zone δc embedded in an elastic body A. The closed curves
with arrowheads, Ŵex, Ŵ and Ŵc, represent three different integration
paths surrounding the process zone δc . The vectors T and u represent
the traction and displacement at a point on the integration path Ŵ,
respectively



(82)

= J (Ŵex ).

Note that the external boundary Ŵex is one of the counterclockwise integration paths Ŵ surrounding the process zone δc , which start from a point on the lower crack
surface − and end at the same point on the upper crack
surface +. From Eq. (82), we can see that J (Ŵex ) gives
the energy release rate Gs by definition. By way of trial,
we take a very thin rectangular integration path Ŵc surrounding the process zone δc as shown Fig. 11. Then,from Eq. (81), we obtain
J (Ŵc ) = −τ (x)∂x D(x)dx
(83)
�

−
with τ = σyz and D = u+z − uz , because nx = 0 andny = ∓1 on ±. Assuming the shear stress τ to be a single-valued function of fault slip D, we can rewrite theabove equation as Dc
J (Ŵc ) =
τ (D)dD,(84)
0where Dc is the critical slip displacement for completingthe process of slip weakening. From this equation, we cansee that J (Ŵc ) gives the fracture surface energy γeff . Since
the J-integral is path independent, from Eqs. (82) and(84), we obtain the following energy balance equation;
 Dc
Gs ≡ J (Ŵex ) = J (Ŵc ) =τ (D)dD ≡ γeff .
(85)
0

The above equation implies that if the relation τ = τ (D)is real, the Griffith-type fracture criterion Gs = γeff is
automatically satisfied. In other words, we can use theshear stress–slip displacement relation τ = τ (D) as a fundamental equation governing the whole rupture process
instated of the Griffith-type fracture criterion.
The residual shear stress τr on  is usually not zero. Insuch a case, as pointed out by Palmer and Rice (1973),

Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Page 15 of 19

Fig. 13 The constitutive relation curves computed from a physicalmodel considering the abrasion of contacting irregular rock surfaces
with the progress of fault slip (cited from Matsu’ura et al. 1992).
Here, c = 2π k1 represents the upper limit of a wavelength range
with same fractal dimension of irregular rock surfaces

Fig. 12 The constitutive relation between shear stress and slip
displacement observed in a stick–slip experiment of rock sample(cited from Ohnaka et al. 1987). The symbols τi, τp, and τr indicate
the initial, peak, and residual stress, respectively, and Dc represents
the critical slip displacement

we need to consider the effects of frictional resistanceto fault slip over . In consequence, the energy balance
Eq. (85) is modified as Dfinal Dfinal DfinalGs +τr dD =
τ (D)dD = γeff +
τr dD
0

0

0

(86)

with

γeff =

 Dc
0[τ (D) − τr ]dD.

(87)

Here, Dfinal (> Dc) denotes the final slip displacementon . In the energy balance Eq. (2) based on the classical crack theory, Kostrov (1974) and also Rivera andKanamori (2005) used the divided expression on theright-hand side of Eq. (86) to represent the work done
for shear faulting. However, from a physical point of
it would be natural to use the unified expression
view,
Dfinalτ (D)dD in the midsection, because the fracture0
surface energy γeff was introduced to handle the effects
of artificial stress singularity at crack tips. Thus, incorporating the fracture surface energy term γeff S into the lastterm of frictionally dissipated energy, we can reduce Eqs.
(2) to (65), which was derived from the general expression of mechanical energy balance in earthquake faulting.
Nowadays, most numerical simulations of earthquake
generation are done by solving a coupled system of theequation of motion and the fault constitutive law under
given initial and boundary conditions. In the 1980s, on
the basis of rock experiments, two different-types offault constitutive law have been proposed; one is therate- and state-dependent law (e.g., Dieterich 1979;

Ruina 1983) and another is the slip-dependent law (e.g.,
Ohnaka et al. 1987). From the viewpoint of energetics,
as discussed above, the slip dependence is more essential in earthquake rupture than the rate (slip velocity)dependence.
Through the stick–slip experiments of rock samples,
Ohnaka et al. (1987) and Ohnaka and Kuwahara (1990)
revealed the existence of a breakdown zone just behindthe propagating rupture front, where the creation of
new fracture surfaces is proceeding. In Fig. 12, we show
a typical example of the constitutive relations between
shear stress and fault slip observed in the breakdown
process; that is, with the progress of fault slip D , theshear stress τ first increases rapidly up to a peak valueτp , and then decreases gradually to a constant level τr .
Matsu’ura et al. (1992) developed a physical model thatrationally interprets the general features of the constitutive relations observed in stick–slip experimentsby considering the abrasion of contacting irregularrock surfaces with the progress of fault slip. We show
the constitutive relation curves computed from sucha physical model in Fig. 13, where c = 2π k1 represents the upper limit of a wavelength range with same
fractal dimension of irregular rock surfaces. From thisfigure, we can see that the critical slip displacement Dc
scales linearly with the upper wavelength limit c of
fractal rock surfaces (Power et al. 1988), which gives arationale of the scale dependence of the fracture surface
energy γeff ≃ 21 (τp − τr )Dc suggested by Shibazaki andMatsu’ura (1998) and Ohnaka and Shen (1999). Thus, inthe 1990s, the slip-weakening fault constitutive law was
introduced as a fundamental equation governing theentire process of earthquake rupture (Matsu’ura et al.
1992; Shibazaki and Matsu’ura 1992).

Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Fig. 14 A schematic diagram showing the spatial density distributionsof traveling seismic wave energy (green) and released elastic strainenergy (blue) at a certain moment t after the termination of concentric
expansion of a circular fault . In this diagram, the azimuthal variationof energy density is ignored

3.3.3 Effects of rupture velocity on energy balance

The dynamic rupture, suddenly started after a quasistatic nucleation process, is at first gradually and afterward rapidly accelerated to a terminal velocity. This is acommon feature of spontaneous rupture propagation,
indicating the excess of the inflow rate of shear strainenergy over the fracture surface energy. The dynamic
rupture accelerated to a terminal velocity is eventuallyarrested by the existence of barriers with large fracturesurface energy (e.g., Aki 1979) in most cases or the shortage of energy supply from the surrounding region (e.g.,
Day 1982) in some cases. To understand the energy flow
and transformation during such a sequential rupture process, the energy balance Eq. (57) would not be helpful,because it has been spatially integrated.
The slip motion with strength (stress) decrease in the
breakdown zone generates dynamic disturbances in its
surrounding region, which are soon separated into specific phases called the near-, intermediate-, and far-fieldterms. As mentioned in Sect. 3.2.2, the far-field term
radiates from the source as traveling waves, while thenear- and intermediate-field terms result in permanent
deformation, which causes the change in elastic strainenergy density there through the interaction with background stress. Figure 14 is a schematic diagram showing
the spatial density distributions of traveling seismic wave
energy (green) and released elastic strain energy (blue) ata certain moment t after the termination of concentric

Page 16 of 19

expansion of a circular fault . The traveling seismic
wave energy goes away from the central part at the Por S-wave velocity. The release of elastic strain energy
density occurs as a result of the interaction between the
stress changes caused by fault slip and a backgroundstress field. As mentioned in Sect. 3.3.1, the released elastic strain energy is partly spent for creating new fracturesurfaces and partly spent for accelerating slip motion following the energy distribution factor p(vr ). The acceleration of fault slip generates further dynamic disturbances,and so the above-mentioned cyclic process, accompaniedwith the radiation of seismic wave energy, continues aslong as the dynamic rupture propagates. After the arrestof dynamic rupture, the disturbances except travelingseismic waves eventually die down. Therefore, in the caseof dynamic rupture, not only the radiated seismic energy
ER but also the work done for shear faulting W must
depend on the time history of rupture growth with faultslip acceleration and deceleration.
As mentioned in Sect. 3.3.2, we can naturally incorporate the fracture surface energy term γeff S into the integral term of frictionally dissipated energy, and so Eq. (2)is reduced to Eq. (65). Then, exchanging the order of integration of the last term in Eq. (65), we obtain the energybalance equation in a stationary state (t ≥ ts), realized
after all the dynamic disturbances except traveling seismic waves died down, as1
ER = 2[τinitial (η) + τfinal (η)]Dfinal (η)dS
�

−



�

dS

 Dfinal (η)(88)τ (D|η )dD

0

with

ER ≡ K (ts ) = 12 ρV (ts )u̇i (x, ts )u̇i (x, ts )dV .

(89)

Here, the subscript “final” means the final stationary state
but not the state at the time of dynamic rupture termination. Even after the arrest of dynamic rupture, the shear
strain energy released in the past at a distant place willcontinue to flow into the source region for a short while.
As a result, the stress inside the rupture area decreases
gradually from a residual level τr (specified by the fault
constitutive law) to a stationary level τfinal (determinedfrom the equation of equilibrium) and the fault slip overshoots. Including such a transient adjustment process,which has been confirmed by Madariaga (1976) througha numerical simulation of a concentrically expanding
circular shear crack, we represent the shear stress–fault
slip relation by τ = τdy (D). The rupture growth pathin the dynamic case is schematically shown in Fig. 15
by the thick grey curve, which is characterized by threeMatsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Fig. 15 A schematic diagram showing the difference in rupture
growth path between the dynamic and quasi-static cases. The
thick gray curve τ = τdy (D) and white curve τ = τqs (D) represent
the dynamic and quasi-static rupture growth paths, respectively.
For reference, the slip-weakening curve τ = τfr (D) is also plottedin black. The symbols τi, τp, τr , and τf indicate the initial, peak, residual,
and final values of the shear stress τ , respectively. The area enclosedby the thick gray and white curves gives the total seismic energy
escaping from the mechanical system concerned

phases: the initial slip-weakening phase, the intermediatehigh-speed rupture propagation phase, and the final slip
deceleration and transient adjustment phase. The energybalance equation corresponding to such a rupture growth
path can be written as Dfinal (η)
dS
τdy (D|η )dD
�
 0
(90)
1
=
2 [τinitial (η) + τfinal (η)]Dfinal (η)dS − ER .
�The above equation states that the dissipated energy byfrictional faulting (W ) = the released elastic strain energy(−UE)—the radiated seismic energy ( ER).
As a reference, we consider the case of quasi-static
rupture growth, where both the rupture velocity vr andthe slip velocity Ḋ are very slow in comparison withelastic wave velocities. In this case, as can be seen from
Eq. (70), the disturbances generated by fault slip containno far-field term (no traveling wave), and so the work
done for shear faulting W balances with the released
shear strain energy UE at any time. The shear stress–
fault slip relation τ = τqs (D) in the quasi-static case ischaracterized by the initial slip-weakening phase andthe subsequent steady rupture growing phase, as shownby the thick white curve in Fig. 15. The energy balance
equation corresponding to such a rupture growth pathcan be written as.
 Dfinal (η)
dSτqs (D|η )dD
�
 0
(91)
1
=
2 [τinitial (η) + τfinal (η)]Dfinal (η)dS.
�Page 17 of 19

Fig. 16 A shear stress–fault slip diagram showing the dynamic rupture
growth path. The symbols τi, τp, τr, and τf indicate the initial, peak,residual, and final values of the shear stress τ , respectively. The difference
between the areas marked with + and − gives the radiated seismicenergy ER

The total release of elastic strain energy (−UE) does
not depend on the time history of rupture growth. Then,from comparison of these two cases, we can finally evaluate the total seismic energy escaping from the mechanical system concerned as Dfinal
−ER S =
[τdy (D) − τqs (D)]dD.
(92)0Here, the integral on the right-hand side corresponds tothe area enclosed by the thick grey and white curves inFig. 15.
From Fig. 15, we can see that the rupture growth pathin the dynamic case is clearly different from that in thequasi-static case. Now, to focus on the rupture growth
path in the dynamic case, we delete the slip-weakening
curve (black) and the quasi-static rupture growth path
(white) from the τ − D diagram in Fig. 15. Then, we
obtain the τ − D diagram as shown in Fig. 16, which corresponds to the τ − D diagram in Fig. 5. In this diagram,the trapezoidal area under the thick broken line represents the total release of elastic strain energy by shear
faulting likewise the diagram in Fig. 5. So, the difference
between the areas marked with + and—gives the radiated seismic energy ER , which is equivalent to the area
enclosed by the thick gray and white curves in Fig. 15.
However, their physical meanings are different, because
the thick broken line is only an auxiliary line but not thequasi-static rupture growth path.

4 Conclusions
We reconsidered the energetics of earthquake faulting
from the current perspective and obtained the following
conclusions. The origin of crustal deviatoric stress is not
the earth’s self-gravitation but the mechanical interaction
between adjacent plates, and so the effect of self-gravitation is negligible in the energetics of shear faulting. Then,Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

in quasi-static problems, the decrease of elastic potential energy (elastic strain energy) balances with the work
done for shear faulting.
The mechanical energy balance of a non-gravitating
earth model in dynamic shear faulting can be generally
written as K (kinetic energy) = –UE (change in elastic
potential energy)—W (work done for shear faulting). The
indirect evaluation of radiated seismic energy ER has
been often done by using a simplified energy balance
equation, written in the average shear stress component τ
parallel to the direction of slip D on a fault with area S asDER S = 21 (τinitial + τfinal )Dfinal − 0 final τ (D)dD. How-ever, this simplified equation appears to be misleading
because the direct evaluation of ER based on elastic dislocation theory shows its dependence on rupture growthrate.
The growth of shear rupture is controlled by the inflow
rate of released shear strain energy into the rupture front.
In the case of quasi-static rupture growth, the work
done for shear faulting W must always balance with thereleased potential energy −UE because of no seismic
wave radiation. Therefore, writing the τ −D relation in thequasi-static case as τ = τqs (D), we obtain the energy balDance equation 12 (τinitial + τfinal )Dfinal = 0 final τqs (D)dD.
In normal earthquakes, the dynamic rupture is gradually accelerated to a terminal velocity, indicating theexcess of the inflow rate of shear strain energy over thefracture surface energy, and then eventually arrestedby the existence of barriers with large fracture surface energy. Even after the arrest of dynamic rupture
propagation, the potential energy released in the past
at a distant place will continue to flow into the source
region for a short while. As a result, the stress level
inside the rupture area decreases gradually and thefault slip overshoots. Including this adjustment process, we write the τ −D relation in the dynamic case asτ = τdy (D). Then, we obtain the energy balance equa Dfinal1τdy (D)dD − ER S .
tion
2 (τinitial + τfinal )Dfinal = 0
Since the total release of elastic potential energy does notdepend on the time history of rupture growth, from comparison of these two cases, we can finally evaluate thetotal seismic energy escaping from the mechanical sysD
tem as −ER S = 0 final [τdy (D) − τqs (D)]dD.
Acknowledgements
I would like to thank Handling Editor and two reviewers for their useful comments and suggestions to improve the original manuscript.
Author contributions
This is a single-authored article.
Funding
No funding was received.

Page 18 of 19

Availability of data and material
Not applicable.

Declarations
Competing interests
I have no competing interest.
Received: 18 August 2023 Accepted: 8 December 2023

References
Aki K (1979) Characterization of barriers on an earthquake fault. J Geophys Res
84:6140–6148
Aki K, Richards PG (1980) Quantitative seismology: theory and methods. W. H.
Freeman and Company, San Francisco
Ampuero J-P, Ripperger J, Mai PM (2006) Properties of dynamic earthquake
ruptures with heterogeneous stress drop. In: Abercrombie R, McGarr
A, Toro GD, Kanamori H (eds) Earthquakes: radiated energy and thephysics of faulting. Geophys Mono Ser, vol 170. AGU, Washington DC, pp
255–261
Andrews DJ (1978) Coupling of energy between tectonic processes andearthquakes. J Geophys Res 83:2259–2264
Backus G, Mulcahy M (1976) Moment tensors and other phenomenologicaldescriptions of seismic sources—I: continuous displacements. Geophys J
R Astr Soc 46:341–361
Burridge R, Knopoff L (1964) Body force equivalents for seismic dislocations.
Bull Seism Soc Am 54:1875–1888
Dahlen FA (1977) The balance of energy in earthquake faulting. Geophys J R
Astr Soc 48:239–261
Das S, Aki K (1977) A numerical study of two-dimensional spontaneous rupture propagation. Geophys J R Astr Soc 50:643–668
Day SM (1982) Three-dimensional simulation of spontaneous rupture: theeffect of nonuniform prestress. Bull Seism Soc Am 72:1881–1902
Dieterich JH (1979) Modeling of rock friction 1: experimental results andconstitutive equations. J Geophys Res 84:2161–2168
Eshelby JD (1956) The continuum theory of lattice defects. Solid State Phys
3:79–144
Eshelby JD (1957) The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc R Soc Lond A241:376–396
Freund LB (1972) Energy flux into the tip of an extending crack in an elastic
solid. J Elast 2:341–349
Freund LB (1979) The mechanics of dynamic shear crack propagation. J Geophys Res 84:2199–2209
Fukuyama E (2005) Radiation energy measured at earthquake source. Geophys
Res Lett 32:L13308. https://​doi.​org/​10.​1029/​2005G​L0226​98
Griffith AA (1921) The phenomena of rupture and flow in solids. Philos Trans R
Soc Lond A221:163–198
Haskell NA (1964) Total energy and energy spectral density of elastic wave
radiation from propagating faults. Bull Seism Soc Am 54:1811–1841
Irwin GR (1957) Analysis of stresses and strains near the end of a crack traversing a plate. J Appl Mech 24:361–364
Kanamori H (1977) The energy release in great earthquakes. J Geophys Res
82:2981–2987
Kanamori H, Rivera L (2006) Energy partitioning during an earthquake. In:
Abercrombie R, McGarr A, Toro GD, Kanamori H (eds) Earthquakes: radiated energy and the physics of faulting. Geophys Mono Ser, vol 170. AGU,
Washington DC, pp 3–13
Kostrov BV (1964) Selfsimilar problems of propagation of shear cracks. J Appl
Math Mech 28:1077–1087
Kostrov BV (1966) Unsteady propagation of longitudinal shear cracks. J Appl
Math Mech 30:1241–1248
Kostrov VV (1974) Seismic moment and energy of earthquakes, and seismic
flow of rock. Izv Acad Sci USSR Phys Solid Earth 1:23–40 (translated toEnglish: Goodspeed F)

Matsu’ura P﻿ rogress in Earth and Planetary Science(2024) 11:2

Kostrov BV (1975) On the crack propagation with variable velocity. Int J Fract
11:47–56
Leigh DC (1968) Nonlinear continuum mechanics. McGraw-Hill, New York
Madariaga R (1976) Dynamics of an expanding circular fault. Bull Seism Soc
Am 66:639–666
Matsu’ura M, Kataoka H, Shibazaki B (1992) Slip-dependent friction law andnucleation processes in earthquake rupture. Tectonophys 211:135–148
Matsu’ura M, Noda A, Terakawa T (2019) Physical interpretation of moment
tensor and the energetics of shear faulting. Tectonophys 771:228228.
https://​doi.​org/​10.​1016/j.​tecto.​2019.​228228
Matsu’ura M, Terakawa T (2021) Decomposition of elastic potential energy anda rational metric for aftershock generation. Geophys J Int 227:162–168
Moakher M (2008) Fourth-order Cartesian tensors: old and new facts, notions
and applications. Q J Mech Appl Math 61:181–203
Mura T (1982) Micromechanics of defects in solids. Martinus Nijhoff Publishers,
Hague
Ohnaka M, Kuwahara Y (1990) Characteristic features of local breakdown neara crack-tip in the transition zone from nucleation to unstable rupture
during stick-slip shear failure. Tectonophys 175:197–220
Ohnaka M, Kuwahara Y, Yamamoto K (1987) Constitutive relations between
dynamic physical parameters near a tip of propagating slip zone during
stick-slip shear failure. Tectonophys 144:109–125
Ohnaka M, Shen L-F (1999) Scaling of the shear rupture process from nucleation to dynamic propagation: implications of geometric irregularity of the
rupturing surfaces. J Geophys Res 104:817–844
Palmer AC, Rice JR (1973) The growth of slip surfaces in the progressive failureof over-consolidated clay. Proc R Soc Lond Ser A 332:527–548
Power WL, Tullis TE, Weeks JD (1988) Roughness and wear during brittle faulting. J Geophys Res 93:15,268-15,278
Rice JR (1968) A path independent integral and the approximate analysis ofstrain concentration by notches and cracks. J Appl Mech 35:379–386
Rice JR, Drucker DC (1967) Energy changes in stressed bodies due to void andcrack growth. Int J Fract Mech 3:19–27
Rivera L, Kanamori H (2005) Representations of the radiated energy in earthquakes. Geophys J Int 162:148–155
Rudnicki JW, Freund LB (1981) On energy radiation form seismic sources. Bull
Seism Soc Am 71:583–595
Ruina AL (1983) Slip instability and state variable friction laws. J Geophys Res
88:10359–10370
Savage JC (1969) Steketee’s paradox. Bull Seism Soc Am 59:381–384
Savage JC, Walsh JB (1978) Gravitational energy and faulting. Bull Seism Soc
Am 68:1613–1622
Shibazaki B, Matsu’ura M (1992) Spontaneous process for nucleation,
dynamic propagation, and stop of earthquake rupture. Geophys Res Lett
19:1189–1192
Shibazaki B, Matsu’ura M (1998) Transition process from nucleation to highspeed rupture propagation: scaling from stick-slip experiments to natural
earthquakes. Geophys J Int 132:14–30
Steketee JA (1958) Some geophysical applications of the elasticity theory of
dislocations. Can J Phys 36:1168–1198
Walpole LJ (1984) Fourth rank tensors of the thirty-two crystal classes: multiplication tables. Proc R Soc Lond A391:149–179

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Page 19 of 19