Reconsideration of the energy balance in earthquake faulting

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 earth-quake 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 origin of 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 balance in dynamic shear faulting from the equation of motion for an elastic body subjected to tectonic-origin deviatoric stress. 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 introduction of 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.


Introduction
The occurrence of earthquakes is now understood as brittle shear fracture at a fault, which releases a part of the elastic potential energy (elastic strain 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, 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 in earthquake faulting from a current perspective.
At the end of the 1950s, Steketee (1958) introduced the elasticity theory of dislocation into geophysics and considered the static change in elastic strain energy by the formation of dislocation (tangential displacement discontinuity = fault slip) in a pre-stressed elastic body bounded by a traction-free surface.He concluded that the total elastic strain energy is always increased by the formation of dislocation irrespective of the initial stress state, which is 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 interpretation of the interaction energy between two stress systems (Eshelby 1956) and demonstrated that the formation of dislocation always decreases the total elastic strain energy 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 the total (elastic and gravitational) potential energy balances with the sum of the work done for shear faulting and the radiated seismic energy.In their formulations, the release of the gravitational potential energy is represented by the volume integral of the scalar product of gravity force and coseismic displacement over the entire earth.Savage and Walsh (1978) found that this volume integral could be transformed into the surface integral of the scalar product of gravity-origin shear stress and fault slip over the 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 the crustal deviatoric stress associated with the horizontal contraction of a homogeneous, isotropic, elastic earth model 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 the earth'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 affect the energetics in earthquake faulting substantially.To sum up, we may use the well-known formula of the static change in elastic strain energy U E due to earthquake faulting for non-gravitating elastic earth models (e.g., Aki and Richards 1980;Rivera and Kanamori 2005), in a good approximation.Here, σ initial ij represents the tectonic-origin initial stress, σ final ij the final stress (the sum of the initial stress and coseismic stress changes), and D final i 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 E R in an ingenious but a bit troublesome way to avoid the problem of gravitational potential energy change as where t initial and t final represent the rupture initiation and termination times, respectively.It should be noted here that 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 han- dle 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 leading edge of crack is not a singular point but a finite process ( 1) (2) zone, where the shear strength (stress) decreases from a peak 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 and Kanamori (2005) and Kanamori and Rivera (2006) derived a simplified expression, where τ is the average shear stress component paral- lel to the direction of slip D on the fault with its area S .The meaning of Eq. ( 3) is obvious; the excess of the released elastic strain energy over the work done for shear faulting gives the radiated seismic energy.Tracking back the history, this type of energy balance equation has been widely used for the evaluation of radiated seismic energy since Kanamori (1977), but something is wrong.Rudnicki and Freund (1981) have directly evaluated the seismic energy radiated from a point dislocation source as Here, M0 represents the second-order time derivative of the seismic moment defined by M 0 = µDS , and ρ , µ , V P , and V S are the density, rigidity, P-wave velocity, and S-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 the attempt 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 edge of 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.

Basic concepts
First of all, we briefly summarize the basic concepts of elastic potential energy (elastic strain energy) and moment tensor.Then, we discuss the representation of earthquake sources and the origin of background stress fields.By the way, throughout this article, Einstein summation convention is used to make mathematical expression concise. (3)

Decomposition of elastic potential energy
In a linearized theory, we can write the potential energy density U E of elastic forces in a quadratic form of the elastic components ε e ij of strain tensor ε ij ; that is, Here, the fourth-order tensor C ijkl denotes the elas- tic stiffness, which has the basic properties of C ijkl = C jikl = C ijlk .Then, by definition, the stress tensor components σ ij are given by the partial derivative of U E with respect to ε e ij as under the condition of C ijkl = C klij .Using this expression, we can rewrite Eq. ( 5) as When the elastic material is isotropic, as demonstrated by Walpole (1984), the constitutive Eq. ( 6) can be decomposed into two independent parts as.
where τ ij and s e ij denote the deviatoric stress and strain tensors defined by This means that the elastic deformation of isotropic material can be decomposed into two different types; that is, volumetric deformation and shearing deformation.As demonstrated by Matsu'ura and Terakawa (2021), substitution of Eq. ( 9) into Eq.( 7) yields with and Here, I 1 (ε e ij ) and I 1 (σ ij ) denote the first invariant of the strain tensor ε e ij and the stress tensor σ ij , respectively, and J 2 (s e ij ) and J 2 (τ ij ) denote the second invariant of the deviatoric strain tensor s e ij and the deviatoric stress tensor τ ij , respectively.From Eq. ( 10), we can see that the elastic potential energy density of isotropic material can be decomposed into two independent parts; that is, the volumetric part and the shearing part. (5) (6) 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 where � • � 2 denotes the Frobenius norm of a second- order tensor.One of the advantages of this representation is that we need not define the reference state to measure the strain ε ij .Another advantage is that the change in elastic potential energy density is also written as the sum of its volumetric part U v E and shearing part (Matsu'ura et al. 2019); that is, denoting the initial stress state and the static stress change from it by The important thing is that there is no interaction between the volumetric and shearing parts, and each of them consists of the interaction term of the stress changes 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 be applicable to the actual problems in the strict sense.However, as demonstrated by Moakher (2008), we can define the closest (in the sense of Euclidean distance) isotropic fourth-order tensor C IE ijkl to a given anisotropic stiffness tensor C A ijkl .The structure of C IE ijkl is the same as that in the isotropic case.The effects of anisotropy are reflected only in the values of two material constants, κ and µ .So, the logical consequences derived in this sec- tion are still applicable to the actual earth.

Physical meaning of moment tensor
Most theoretical problems in seismology are treated in 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 (13) �U earth and the linear elastic model, Backus and Mulcahy (1976) introduced the concept of moment tensor, which includes the seismic moment tensor (Kostrov 1974) as a special case.The following is a physical interpretation of the Backus-Mulcahy moment tensor based on the theoretical 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 S ex with its outer unit normal vector n = (n 1 , n 2 , n 3 ) in a Cartesian coor- dinate system (x 1 , x 2 , x 3 ) as shown in Fig. 1.In a line- arized theory, irrespective of the stress-strain relation of material, we can write the equation of motion of the solid body for t ≥ 0 as Here, üi represents the second-order time derivative of the i-component of a displacement vector u i , σ ij the ij -component of a stress tensor, and f i the i-component of a body force vector without time variation.The symbol ∂ j denotes partial differentiation with respect to the spa- tial coordinate x j .Let's suppose that the solid body was in an equilibrium state, 0 = ∂ j σ 0 ij (x) + f i (x) , for t < 0 .Then, in order to solve the equation of motion, we need to specify the relationship between the stress change �σ ij ( = σ ij − σ 0 ij ) and the displacement u i through the strain change �ε ij defined geometrically as 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 the region surrounding the defects.Thus, the strain �ε ij defined by Eq. ( 18) is not necessarily elastic anytime and (17) anywhere; that is, denoting the inelastic part of the strain �ε ij by �ε a ij , Then, from Eqs. ( 6) and ( 19), the stress-strain relation in the actual earth can be written as We now consider the occurrence of some inelastic event at t = 0 within the solid body in static equilib- rium.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 and the stress-strain relation (20) as with Here, �σ m ij is a model stress, calculated from the displacements u i through the geometrical strain under the assumption of linear elasticity, and m ij is called the moment tensor density, which is equivalent to the stress glut tensor, Ŵ ij ≡ �σ m ij − �σ ij , defined by Backus and Mulcahy (1976).Substituting Eq. ( 22) in Eq. ( 21), we obtain the equation of motion to be solved in the framework of linear elasticity as or with f m i (x, t) = −∂ j m ij (x, t) , which is called the body force equivalent (Burridge and Knopoff 1964).The inelastic strain referred to here is almost the same as the eigenstrain in Mura (1982), but slightly different from the 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 kl (x) itself is physical substance, but the moment (19) (20) tensor density m ij (x) defined by Eq. ( 24) is not physical substance.As can be seen from Eq. ( 25), the moment tensor density appears as a virtual source term in the theory of linear elasticity.If the distribution of inelastic strain changes is restricted within a region V s (Fig. 1), by definition, the moment tensor M ij is given by the volume integral of C ijkl �ε a kl (x) over V s as

Representation of earthquake sources
In the mid-1960s, it was convinced that the source of earthquakes 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 M ij defined by Eq. ( 27) can be written as In general, displacement discontinuity at the interface �(η) is represented as , and H(x − η) denotes the Heaviside step function.Then, denoting the Dirac delta function by δ(x − η) and a unit normal vector pointing from the backside of �(η) to the front side by n − (η) (Fig. 1), we can represent the distribution of inelas- tic strain �ε a ij (x) localized at �(η) as The isotropic part of the above inelastic strain is and so the deviatoric part can be written as 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): (27) Here, m � ij (η) represents the area density of moment ten- sor on �(η).
In the case of shear faulting, denoting the tangential component of displacement discontinuity 2), we can reduce Eq. ( 33) to 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 M kk of the seismic moment tensor is always zero.On the other hand, in the case of crack opening, denoting the normal component of displacement discontinuity i , we can reduce Eq. ( 33) to which has a nonzero isotropic part 1 3 M kk = κ � 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. (33)

Origin of background stress fields
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 part and the shearing part.This implies that we can treat the isotropic 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 question is the origin of the deviatoric part (tectonic stress).Our understanding is that the origin of tectonic stress is not the gravitational contraction of the elastic earth but the stable process of converting thermal energy into kinetic energy in the 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 background tectonic 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 S ex (Fig. 1).In the initial state, the elastic 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 ε a ij occurred within a region V s ( ⊂ V ).Then, the elastic poten- tial energy U E ≡ V U E dV produced by this process is evaluated as because Note that σ ij n j = 0 on S ex (the boundary condition of traction free) and When the source region V s consists of two independent parts, V s0 and V s1 (Fig. 3), applying the formula (37), we can write the produced elastic potential energy as ( 37) The proof of Eq. ( 40) is given as follows: and Following Andrews (1978), we regard ε a0 ij 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 read as the elastic potential energy after the occurrence of a sudden faulting with inelastic strain distribution ε a1 ij in the pre-stressed state of σ 0 ij .Therefore, changing the notations of the inelastic strain ε a1 ij , the stress change σ 1 ij , and the source region V s1 associated with the sud- den faulting to �ε a ij , �σ ij , and V s , respectively, we can rewrite Eq. (39) as with (40) (41) Here, U E 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 inelastic strain �ε a ij can be written as Substituting this expression in Eq. ( 44), we obtain which is nothing but Eq. ( 1).

Energetics in earthquake faulting
As a starting point of discussion, we directly derive a basic equation of mechanical energy balance in dynamic 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 seismic energy based on a simplified energy balance equation and the direct evaluation based on the analytical solution 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 a slip-weakening fault constitutive law from the viewpoint of energy balance.

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 and tectonic-origin deviatoric stress.So, the strain ε ij appear- ing in this section is elastic.The elastic body occupies the region V bounded by a traction-free external surface S ex and an internal surface S in enclosing a fault as shown in Fig. 1, where the direction of unit normal vectors n to the external 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 as ( 46) where ui is the i-component of a velocity vector, f g i is the i-component of a gravity force vector, and the symbol ∂ t represents 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 stress σ 0 ij (x) defined in Sect.2.4 and the gravity force f g i (x) .So, the stress appearing in Eq. ( 48) should be read as the sum of the background tectonic stress and the stress change measured from it; that is, σ ij (x, t) = σ 0 ij (x) + �σ ij (x, t).Now, following Sect.7.8 of Leigh (1968), we multiply both sides of Eq. ( 48) by ui : Here, the first term on the right-hand side can be rewritten as with εij = 1 2 (∂ j ui + ∂ i uj ) .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: Furthermore, applying the Gauss divergence theorem to the first term on the right-hand side and replacing the stress tensor σ ij in the second term with C ijkl ε kl , we can rewrite Eq. ( 51) as Here, we dropped the term S ex ui σ ij n j dS because the traction on the external surface S ex is always zero.Note that the left-hand side of Eq. ( 52) represents the increase rate of kinetic energy, the first term on the right-hand side the rate of work done by the traction acting on the internal surface S in , and the second and third terms the rates of change in elastic and gravitational potential energy, respectively.
In the case of shear faulting, shrinking the internal surface S in to a fault with its front side + and back side − and denoting the i-component of the tangential velocity jump ui ≡ u+ i − u− i across by Ḋi , we can rewrite the first term on the right-hand side of Eq. ( 52) as ( 49) As for the last gravity force term, applying the idea of interaction energy in Sect.2.4, Savage and Walsh (1978) have demonstrated that the volume integral over the whole elastic region V could be transformed into the sur- face integral over the fault as follows: where σ g ij represents the internal stress field that satisfies the equation of equilibrium 53) and (54) into Eq.( 52) yields When the internal stress field caused by the gravity force is isotropic ( σ ), the last term on the right-hand side of Eq. ( 55 We consider a shear faulting that initiated at t = t initial on a fault and terminated at t = t final .Hereafter, t initial and t final are abbreviated to t i and t f , (53) respectively, unless misleading.Let us suppose that the disturbances caused by the shear faulting died down by t = t stationary ( t s in abbreviation) except the seismic waves traveling outward from the source.Then, integrating each term of Eq. ( 57) with respect to t from t i to t s ( >> t f ), we obtain with and, from Eqs. ( 44) and ( 47), where represents the final rupture area �(t f ) .In this context, we can read K (t s ) as the radiated seismic energy E R , and so Eq. ( 61) is written as 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 dissipated energy (Fukuyama 2005).

Radiated seismic energy
Nowadays, using an inversion technique, we can estimate the dynamic process of earthquake faulting from seismological and/or geodetic data.A direct way to evaluate the radiated seismic energy of an earthquake is of course to perform the integration in Eq. ( 62).Another indirect way is to evaluate the right-hand side of Eq. ( 65), which relates the radiated seismic energy and the history of fault slip during the earthquake.Rivera and Kanamori (2005) simplified the energy balance Eq. ( 2) by exchanging the order of integration of the last frictionally dissipated energy term as ( 61)

Indirect evaluation of radiated seismic energy
Here, D denotes the vector magnitude of fault-slip, and τ = σ ij n − j ν i represents the scalar projection of the trac- tion vector T i = σ ij n − j in the slip direction ν i = D i D (Fig. 4).In this expression, the shear stress τ is assumed to be a single-valued function of the fault slip D ; that is τ = τ (D) .Then, we can replace the integration In the natural process of shear faulting, the shear stress τ defined above is always positive.Kanamori and Rivera (2006) further simplified the energy balance equation by neglecting the fracture surface energy term γ eff S and taking the spatial average of the shear stress τ (η) and the fault slip D(η) over the fault with its final dimension S as which is nothing but Eq. ( 3).This energy balance equation is often interpreted by using the τ −D diagram in (66 marked with + and -) gives the seismic wave energy E R radiated from the source.Such interpretation seems to be plausible; the most primitive one is found in Kanamori (1977), which gives a theoretical basis of the moment magnitude scale M w of seismic events.But something is wrong.For instance, when the growth of shear rupture is quasi-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.

Seismic waves and permanent deformation
The displacement field due to a moment tensor M pq located at the origin ( x = 0 ) of a Cartesian coordinate system is generally expressed as Here, G ip (x, t; 0, t ′ ) is the solution of the equation of motion to a unit impulsive force in the x p -direction, δ ip δ(x)δ(t − t ′ ) , applied to a point x = 0 at t = t ′ , and ∂ q G ip (x, t; 0, t ′ ) represents the partial derivative of G ip (x, t; 0, t ′ ) with respect to the spatial coordinate x q .
The explicit form of G ip for a homogeneous, isotropic, unbounded elastic medium, which is called Stokes' solution, is given in Sect.4.2 of Aki and Richards (1980) as (68) where r = x 2 1 + x 2 2 + x 2 3 and γ i = x i r ( = ∂ i r ) denote the 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, M pq = M 0 (t)(ν p n − q + ν q n − p ) , into Eq.( 68), we obtain the explicit expression of displacement field due to a seismic moment tensor as where a = γ • n − , b = γ • ν , and M 0 (t) is a cumulative seismic moment function.
From the difference in geometrical attenuation, the first term on the right-hand side of Eq. ( 70) is traditionally called the near-field term, the second and third terms the intermediate-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 distance r .Since the cumulative seismic moment function M 0 (t) monotonically increases with time t from 0 at t = t i up to M 0 at t = t f , its time derivative Ṁ0 (t − r/c) vanishes for t − r/c > t f .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 ( V P for P-wave and V S for S-wave).On the other hand, the near- and intermediate-field terms remain as permanent deformation after disturbance died down ( t − r/c > t f ), which causes the change in elastic potential energy.
During the rupture growth ( t i ≤ t ≤ t f ), the far-field terms 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 energy change �U E (t) for a sufficiently large time t ( ≥ t s ) as in Eqs. ( 62), (63), and (64).As discussed in Sect.3.3, the last (70)

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 x 1 -axis is taken in the direction of fault slip (then a ν = sin θ cos φ ) and the x 3 -axis in the direction of fault normal (then a n = cos θ ).From Eq. ( 70), we can see that the particle motion of the far-field P-wave (the fourth term) is perpendicular to the wave front, while that of the far-field S-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 mutually orthogonal particle velocities as (71a) where e r , e θ and e φ are the mutually orthogonal unit vec- tors that define a local coordinate system at a receiver point (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 t s ( >> t f ), the cumulative kinetic energy K (t s ) is given by the volume integral of the kinetic energy density 1 2 ρ| u| 2 over the sphere with the radius r = c(t s − t i ): with Substituting the explicit expressions of the mutually orthogonal particle given by Eq. ( 71) into Eq.( 72) and reducing the interval of integral for r to c(t s − t f ) ≤ r ≤ c(t s − t i ) , we can evaluate the contribu- tion to the total kinetic energy of each particle motion as follows: (71b) (74a) where T = t f − t i denotes the rupture duration.Finally, summing up the above expressions, we obtain Rudnicki and Freund (1981) has obtained the equivalent result in a slightly different way, though their radiated energy in Eq. ( 4) is just twice as large as the present evaluation, because they adopted the definition of radiated energy by Haskell (1964), which is essentially the same as K (t s ) 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 Ṁ0 (t) but the moment rate change M0 (t).
The seismic moment function M 0 (t) is defined as As pointed out by Ampuero et al. ( 2006), the fault slip distribution of actual earthquakes is heterogeneous in both time and space because of the spatiotemporal heterogeneity of stress drop.So, the moment rate function Ṁ0 (t) of actual earthquakes is not always unimodal.Nevertheless, in the following, we consider the case where the seismic moment monotonically increases with time for sake of simplicity.First, to compare the direct 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 M 0 (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) .
( 75) When the fault area concentrically expands from 0 at t = t i to S at t = t f , we may define the average rupture velocity as v r = √ S/2π �T with T = t f − t i .Figure 7 schematically shows the rupture velocity-dependence of the fault slip D , slip velocity Ḋ , and slip acceleration D , where the time t is normalized by the rupture duration �T = √ S/2π v r .From this diagram, we can see that the fault slip acceleration is proportional to the square of the average rupture velocity v r .On the contrary, the simpli- fied Eq. ( 67) appears to be independent of the rupture velocity.

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 of clack 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 strength of actual materials is finite.In the 1980s, the stick-slip experiments of rock samples revealed the existence of a 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, which is 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).

Griffith-type fracture criterion
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 G s (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 unstable rupture 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 total potential energy as the decrease in elastic strain energy.In the 1950s, on the basis of the analysis of stress and strain fields near crack tips, Irwin (1957) related the energy release rate G s to the stress intensity factor K s (the factor indicating the intensity of stress singularity at crack tips) as G s = cK 2 s 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 G s and K s 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 case of dynamically propagating cracks by considering the inflow rate of released elastic strain energy into the crack tip (e.g., Kostrov 1964Kostrov , 1966Kostrov , 1975;;Freund 1972Freund , 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 is given 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 v r .Figure 8 is a snapshot of stress and slip distributions along a crack surface ( x < c ) and its extension ( x > c ) at a moment t ( ≥ 0 ).In this case, the energy balance at the crack tip x = c(t) can be written as with (77) Here, the left-hand side of Eq. ( 77) represents the inflow rate of released elastic strain energy into the crack tip propagating at v r , and the right-hand side is the fracture surface energy at x = c(t) .As can be seen from Eq. ( 78), K inst 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 v r .When the crack extension is quasi-static ( v r = 0 ), the left-hand side of Eq. ( 77) corresponds to the static energy release rate G s , because the energy distribution factor becomes unity and the instantaneous stress intensity factor K inst is identical to the static stress intensity factor K s .As the rupture velocity v r increases from 0 to V S , the factor p(v r ) decreases monotonically from 1 to 0 (Fig. 9).This means that only a part ( p ) of the supplied elastic strain energy ( K 2 inst 2µ ) is spent for creating new crack area, and the remaining 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 K inst and γ eff , the factor p(v r ) must take an appropriate value at each moment.In other words, the rupture velocity v r 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 whole rupture process; a good numerical example is given in Das and Aki (1977).

J-integral and slip-weakening fault constitutive law
In the classical crack theory, as shown in Fig. 8, the shear stress τ just outside a crack has a 1 √ x − c(t)-type sin- gularity, but the fault slip D is zero everywhere.On the Fig. 8 A snapshot of stress and slip distributions along a shear crack extending in the x-direction at a rupture velocity v r .The red and blue lines represent the stress and fault slip distributions at a moment t , respectively Fig. 9 Dependence of the energy distribution factor p on the instantaneous rupture velocity v r 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 singular point, we obtain 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 resid- ual 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 shown in Fig. 11.Here, the area A surrounded by the external boundary ( Ŵ 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 per unit thickness is given by ( 79) where U E = 1 2 σ ij ε ij , T i = σ ij n j , and u i and n i denote the i -components of a displacement vector and a unit normal vector, respectively.Using a path independent integral J (Ŵ) defined by Rice (1968), we can evaluate the decrease of the total potential energy P by quasi-static unit crack extension as Note that the external boundary Ŵ ex is one of the coun- terclockwise integration paths Ŵ surrounding the pro- cess 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 G s by definition.By way of trial, we take a very thin rectangular integration path Ŵ c sur- rounding the process zone δc as shown Fig. 11.Then, from Eq. ( 81), we obtain with τ = σ yz and D = u + z − u − z , because n x = 0 and n y = ∓1 on ± .Assuming the shear stress τ to be a sin- gle-valued function of fault slip D , we can rewrite the above equation as where D c is the critical slip displacement for completing the process of slip weakening.From this equation, we can see 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; The above equation implies that if the relation τ = τ (D) is real, the Griffith-type fracture criterion G s = γ eff is automatically satisfied.In other words, we can use the shear stress-slip displacement relation τ = τ (D) as a fun- damental equation governing the whole rupture process instated of the Griffith-type fracture criterion.
The residual shear stress τ r on is usually not zero.In such a case, as pointed out by Palmer and Rice (1973), (83) we need to consider the effects of frictional resistance to fault slip over .In consequence, the energy balance Eq. ( 85) is modified as with Here, D final ( > D c ) denotes the final slip displacement on .In the energy balance Eq. ( 2) based on the classical crack theory, Kostrov (1974) and also Rivera and Kanamori (2005) used the divided expression on the right-hand side of Eq. ( 86) to represent the work done for shear faulting.However, from a physical point of view, it would be natural to use the unified expression D final 0 τ (D)dD in the midsection, because the fracture 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 last term 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 the equation 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 of fault constitutive law have been proposed; one is the rate-and state-dependent law (e.g., Dieterich 1979; (86) 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 behind the 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 , the shear 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 that rationally interprets the general features of the constitutive relations observed in stick-slip experiments by considering the abrasion of contacting irregular rock surfaces with the progress of fault slip.We show the constitutive relation curves computed from such a physical model in Fig. 13, where c = 2π k 1 repre- sents the upper limit of a wavelength range with same fractal dimension of irregular rock surfaces.From this figure, we can see that the critical slip displacement D c scales linearly with the upper wavelength limit c of fractal rock surfaces (Power et al. 1988), which gives a rationale of the scale dependence of the fracture surface energy γ eff ≃ 1 2 (τ p − τ r )D c suggested by Shibazaki and Matsu'ura (1998) and Ohnaka and Shen (1999).Thus, in the 1990s, the slip-weakening fault constitutive law was introduced as a fundamental equation governing the entire process of earthquake rupture (Matsu'ura et al. 1992;Shibazaki and Matsu'ura 1992).Here, c = 2π k 1 represents the upper limit of a wavelength range with same fractal dimension of irregular rock surfaces

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 a common feature of spontaneous rupture propagation, indicating the excess of the inflow rate of shear strain energy over the fracture surface energy.The dynamic rupture accelerated to a terminal velocity is eventually arrested by the existence of barriers with large fracture surface 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-field terms.As mentioned in Sect.3.2.2, the far-field term radiates from the source as traveling waves, while the near-and intermediate-field terms result in permanent deformation, which causes the change in elastic strain energy density there through the interaction with background stress.Figure 14 is a schematic diagram showing the spatial density of traveling seismic wave energy (green) and released elastic strain energy (blue) at a certain moment t after the termination of concentric 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 background stress field.As mentioned in Sect.3.3.1, the released elastic strain energy is partly spent for creating new fracture surfaces and partly spent for accelerating slip motion following the energy distribution factor p(v r ) .The accelera- tion of fault slip generates further dynamic disturbances, and so the above-mentioned cyclic process, accompanied with the radiation of seismic wave energy, continues as long as the dynamic rupture propagates.After the arrest of dynamic rupture, the disturbances except traveling seismic waves eventually die down.Therefore, in the case of dynamic rupture, not only the radiated seismic energy E R 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 inte- gral 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 energy balance equation in a stationary state ( t ≥ t s ), realized after all the dynamic disturbances except traveling seismic waves died down, as with 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 will continue 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 (determined from the equation of equilibrium) and the fault slip overshoots.Including such a transient adjustment process, which has been confirmed by Madariaga (1976) through a numerical simulation of a concentrically expanding circular shear crack, we represent the shear stress-fault slip relation by τ = τ dy (D) .The rupture growth path in the dynamic case is schematically shown in Fig. 15 by the thick grey curve, which is characterized by three (88 ui (x, t s ) ui (x, t s )dV .As a reference, we consider the case of quasi-static rupture growth, where both the rupture velocity v r and the slip velocity Ḋ are very slow in comparison with elastic wave velocities.In this case, as can be seen from Eq. ( 70), the disturbances generated by fault slip contain no far-field term (no traveling wave), and so the work done for shear faulting W balances with the released shear strain energy U E at any time.The shear stress- fault slip relation τ = τ qs (D) in the quasi-static case is characterized by the initial slip-weakening phase and the subsequent steady rupture growing phase, as shown by the thick white curve in Fig. 15.The energy balance equation corresponding to such a rupture growth path can be written as.The total release of elastic strain energy ( − U E ) 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 Here, the integral on the right-hand side corresponds to the area enclosed by the thick grey and white curves in Fig. 15.
From Fig. 15, we can see that the rupture growth path in the dynamic case is clearly different from that in the quasi-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 cor- responds 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 E R , 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 the quasi-static rupture growth path.

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,   τ qs (D)dD .In normal earthquakes, the dynamic rupture is gradually accelerated to a terminal velocity, indicating the excess of the inflow rate of shear strain energy over the fracture surface energy, and then eventually arrested by 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 the fault 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-

Fig. 1
Fig. 1 Geometry and notation of an elastic body V bounded by a traction-free external surface S ex and an internal closed surface S in .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

Fig. 2
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

Fig. 3
Fig. 3 Two independent source regions V s0 and V s1 in an elastic body V bounded by a traction-free external surface S ex .σ ij represents a common stress field ) vanishes: because D i n − i = 0 for shear faulting.Then, discarding the last term associated with gravity force, we obtain a basic equation of mechanical energy balance in dynamic shear faulting as with the increase rate of kinetic energy the rate of work done for shear faulting and the rate of change in elastic potential energy

Fig. 5
as follows: (1) the trapezoidal area under the thick broken line, which corresponds to the first term on the right-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, and 4) the remaining part (the difference between the areas

Fig. 4
Fig. 4Geometry 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

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

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

Fig. 7
Fig. 7 A schematic diagram showing the rupture-velocity dependence of the fault slip D , slip velocity Ḋ , and slip acceleration D .The vertical axis 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) − σ xz ∂ x u z )n x − (σ yz ∂ x u z )n y dŴ,

Fig. 10 A
Fig.10A 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 τ p to a residual value τ r with the increase of fault slip D

Fig. 12
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 D c represents the critical slip displacement

Fig. 14 A
Fig. 14 A schematic diagram showing the spatial density distributions of traveling seismic wave energy (green) and released elastic strain energy (blue) at a certain moment t after the termination of concentric expansion of a circular fault .In this diagram, the azimuthal variation of energy density is ignored initial (η) + τ final (η)]D final (η)dS − E R .initial (η) + τ final (η)]D final (η)dS.

Fig. 15 A
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 plotted in black.The symbols τ i , τ p , τ r , and τ f indicate the initial, peak, residual, and final values of the shear stress τ , respectively.The area enclosed by the thick gray and white curves gives the total seismic energy escaping from the mechanical system concerned

Fig. 16 A
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 seismic energy E R initial + τ final )D final = D final 0 τ dy (D)dD − E R S .Since the total release of elastic potential energy does not depend on the time history of rupture growth, from comparison of these two cases, we can finally evaluate the total seismic energy escaping from the mechanical system as −E R S = D final 0 [τ dy (D) − τ qs (D)]dD.