Post-seismic relaxation theory on a laterally heterogeneous viscoelastic model

Summary

1 Introduction

Time-dependent deformation of the Earth after an imposed stress step includes numerous processes spanning many time scales. In modelling of crustal deformation phenomena, the principal processes are the static elastic offset and the subsequent relaxation of Earth's ductile regions over periods ranging from days to years. In the case of relaxation after an earthquake (e.g. Elsasser 1969; Nur & Mavko 1974; Anderson 1975), an earthquake occurring in an idealized elastic region in the upper crust is followed by time-dependent readjustment as underlying viscoelastic material relaxes in response to the imposed static stresses. Modelling of post-seismic relaxation in such a geometry has gained in importance in recent years due to the constraints provided by space-based geodetic data, which has allowed estimation of depth-dependent viscosity structure in numerous regions. Such models have served well the purpose of representing the first order viscoelastic structure but have been limited because they are usually laterally homogeneous, i.e. viscosity (and all other viscoelastic parameters) varies only as a function of depth in a flat- or spherically-layered structure. For the post-seismic relaxation problem and many other purposes, it is desirable to account for laterally varying viscoelastic structure. Especially for continental deformation (where most measurements are available), lateral variations in elastic plate thickness and temperature at depth, which is likely to be reflected in the mechanical behaviour of a ductile lower crust and upper mantle, are exhibited at horizontal scales of 15 km and greater (Lowry et al. 2000). Similar lateral heterogeneity in viscoelastic structure is implied by large-scale continental seismic tomography. The 3-D S-wave velocity structure around eastern Asia down to a depth of 720 km imaged by Friederich (2003) reveals strong lateral variations in velocity structure at a horizontal scale of 100–3000 km which, if due to thermal perturbations, would apply lateral variations in viscosity of several orders of magnitude.

This paper seeks to model the response of continental or oceanic lithosphere to earthquakes or to applied loads (e.g. the ongoing India–Eurasia collision), with account for vertically-varying and/or laterally varying material properties. A survey of current methods to describe this response can be described with three categories:

• Fully analytic. An example is the response of an elastic layer over a homogeneous viscoelastic half-space in the case of strike-slip faulting along an infinitely long, straight fault (Nur & Mavko 1974; Savage & Prescott 1978).

• Semi-analytic. This is defined to include spectral methods in either a half-space (e.g. Rundle 1980; Iwasaki 1985) or spherical (Pollitz 1992; Piersanti et al. 1997; Pollitz 1997) geometry, designed for a vertically-layered earth model. The response may be calculated exactly for any single spectral component, but integration over wavenumber or summation over spherical harmonic degrees is necessary in order to obtain the full response.

• Fully numerical. The finite element method has found several applications to post-seismic relaxation problems (e.g. Wahr & Wyss 1980; Melosh & Raefsky 1980; Cohen 1982; Lyzenga et al. 1991; Saucier & Humphrieys 1993; Hearn 1998; Freed & Lin 2001; Hearn et al. 2002).

The methods employed in each category carry respective advantages and disadvantages. Calculations with fully and semi-analytic solutions are obviously much faster than those with fully numerical solutions. The semi-analytic solutions allow for the inclusion of finite faults within a vertically-layered medium, and the same computational effort is required for both small and large spatial domains. However, they can not account for lateral variations in viscoelastic properties, and they are limited to a linear constitutive relationship. Fully numerical solutions overcome these limitations but are computationally more expensive, especially as the size of the spatial domain increases.

A natural extension of the semi-analytic approach is to use seismic scattering theory (e.g. Hudson & Heritage 1981) to allow inclusion of lateral variations in viscoelastic properties. In the context of the normal mode theory of wave propagation, this method has found application to both Born (i.e. single-scattering) scattering (e.g. Gaherty et al. 1996) and multiple scattering (e.g. Friederich 1999) problems. The method takes advantage of the fact that the normal mode solution for a 1-D earth model can be rapidly calculated; to construct the solution on the 3-D model requires essentially numerous evaluations of wavefields on the 1-D model, and the numerical effort scales roughly with the size of the laterally heterogeneous region. In the case of post-seismic relaxation of an elastic-viscoelastic stratified medium, the normal mode solution of Pollitz (1992, 1997) for 1-D stratification is amenable to the same extension to the 3-D case.

A viscoelastic normal mode approach for relaxation on a 3-D viscoelastic model has been previously developed by Tromp & Mitrovica (2000) for postglacial rebound applications. That work differs from the present one in several respects: it is designed for surface loads rather than buried seismic sources; it is a global treatment employing splitting matrices of very long wavelength structure rather than a local treatment employing scattering integrals over local structure.

A local approach based on perturbation theory was developed by Du et al. (1994) and Cervelli et al. (1999) for the problem of the static response to seismic sources on a laterally heterogeneous half-space. Their work bears many similarities to the present one in the use of scattering theory and the use of a representation theorem as an expedient to integrate over the relevant region of equivalent sources. Besides the restrictions to the static case and a ‘zeroth order’ solution for a 2-D homogeneous elastic half-space, it differs from the present treatment in the use of direct Green functions, rather than normal modes, to represent the deformation field.

In this paper a scattering-based solution is developed to the problem of post-seismic relaxation on a 3-D laterally heterogeneous viscoelastic model, starting from the normal mode solution of the 1-D case. The present treatment considers only lateral variations in isotropic elastic parameters and viscosity, and account is made for the first order effect of Earth's gravitational acceleration. Following the development of the solution, examples are presented to illustrate its effectiveness in problems of post-seismic crustal deformation where lateral variations in material properties may be encountered.

2 Normal Mode Representation Of Post-Seismic Relaxation

2.1 Response On Laterally Homogeneous Earth Model

Let elastic parameters {μ₀(r), κ₀(r)}, density ρ₀(r), and viscosity η₀(r) define a spherically layered earth model in a spherical r−θ−φ geometry, where μ₀, κ₀, ρ₀, and η₀ are the depth-dependent shear modulus, bulk modulus, density, and viscosity, respectively. Let c₀ be the associated elastic tensor (Aki & Richards 1980). It is most convenient to work in the Laplace transform domain. The Laplace transform of a function with independent variable s is

(1)

It is assumed that constitutive relation between tensor strain rate and stress is linear, as is the case for the well-known Maxwell viscoelastic solid, standard linear solid, and Burgers viscoelastic solid (e.g. Yuen & Peltier 1982). Viscoelasticity is realized by application of the correspondence principle, through which the elastic parameters are replaced with s-dependent elastic parameters, e.g.

(2)

for a Maxwell viscoelastic solid.

The combined elastic and post-seismic response of a spherically layered viscoelastic earth model to a seismic source is expressed as a sum of normal modes. Assuming the seismic source is a point source in space acting at point r₀, the displacement field on the laterally homogeneous viscoelastic model may be represented as a superposition of vector spherical harmonics associated with individual spheroidal and toroidal modes (Gilbert & Dziewonski 1975; Pollitz 1997). In order to describe both the instantaneous elastic response and post-seismic deformation following an earthquake, both ‘static’ and ‘viscoelastic’ modes must be included in the formulation. The displacement field may be written as

(3)

where Φ⁰_j is a potential function which depends upon the source properties as well as the source-observation point geometry, which determines the sense of excitation of mode j at the observation point. The sum in (3) is over all spheroidal (j=S) and toroidal (j=T) static and viscoelastic modes in the basis set. The operator O^j is defined by

(4)

(5)

where U and V are spheroidal mode displacement eigenfunctions and W is the toroidal mode displacement eigenfunction, each of which depends on the viscoelastic stratification, the characteristic inverse decay time of the considered mode s_S or s_T, and the corresponding spherical harmonic degrees l_S or l_T, as well as the radius of observation. is the surface gradient operator.

In general, the static mode basis set consists of a countably infinite number of dispersion branches, i.e. those associated with Earth's free oscillations. For a spherically stratified model with a finite number of homogeneous layers, the viscoelastic mode basis set consists of a finite number of dispersion branches (Pollitz 1997). The sum over j then represents a double sum over all spherical harmonic degrees and the dispersion branches. In order to represent deformation fields up to a wavelength of, say, 16 km (0.0251 radians on an earth model of radius 6371 km), requires a spherical harmonic expansion up to a maximum total degree of 2π/0.0251 ∼ 2500.

The potential function Φ⁰_j is given by

(6)

In (6), M is the moment tensor, F is a single force. Under the geometry of Fig. 1, E_j is the excitation function tensor () for moment tensor excitation involving products of the strain tensor associated with mode j evaluated at r₀ and corresponding geometrical functions, E′_j is the excitation function vector () for force excitation. The function ψ_j(s) is given by

(7)

For a static mode j, ε_j is the kinetic energy integral:

where ρ(r) is density and ω_j is the angular frequency of the jth degenerate free oscillation on the spherical earth model. For a viscoelastic mode j, ε_j is a group ‘diffusion rate’ for mode j defined by eq. (19) of Pollitz (1997) and eq. (30) of Pollitz (1992):

(8)

The integration is carried out from radius 0 to Earth's radius R. It is assumed in (8) that relaxation occurs only in shear, not in the bulk modulus. We shall assume henceforth that the mode eigenfunctions U, V, and W are normalized such that ε_j= 1 for all modes j.

1

Scattering geometry for scatterer —observation point interactions, with angular distance θ and azimuth γ. The scattering region on the unit sphere is denoted by Γ. For source —scatterer interactions, the angular distance and azimuth are Δ and φ, respectively.

Open in new tabDownload slide

The foregoing results are valid for any time history of earthquake faulting embodied in M(r₀; s) (similarly for a history of forcing F(r₀; s)). In the case of an impulsive seismic source, M(r₀; s) ∼ 1/s, in which case ψ_j(s) in eq. (7) behaves as ∼1/s for a static mode and 1/[s (s+s_j)] for a viscoelastic mode. From eqs (3) and (6), this yields a displacement field with time-dependent behaviour like ∼H (t) (a step function acting at time t= 0) for any static mode and ∼1 − exp[−s_jt] for a single viscoelastic mode. It is evident that the ‘static’ modes on a laterally homogeneous earth model coincide with the purely elastic displacement field (i.e. the static offset immediately after a seismic event), and the ‘viscoelastic’ modes on the laterally homogeneous earth model coincide with the post-seismic displacement field (i.e. the evolving deformation field at time t > 0). It will become clear in the next section that the viscoelastic modes on the laterally heterogeneous Earth contribute to only the post-seismic displacement field, but the static modes on the laterally heterogeneous Earth contribute to both the elastic and post-seismic displacement fields because of the effects of mode coupling.

2.2 Response On Laterally Heterogeneous Earth Model

Let a laterally heterogeneous earth model be defined in terms of a density distribution and an elastic tensor c decomposed into spherically symmetric and aspherical parts:

(9)

The response of a laterally heterogeneous earth model to the same seismic excitation as in the previous section can be expressed as

(10)

(11)

Eqs (3), (10), and (11) define a decomposition of the total displacement into the incident displacement u⁰ and a scattered displacement u¹:

(12)

(13)

Seismic scattering theory specifies the scattered displacement u¹ to be the convolution of responses to equivalent seismic sources distributed over the volume V of the Earth. In terms of the potential functions this takes the form ()

(14)

where δm is a moment tensor density given by

(15)

and δf is a single force density given by

(16)

where g₀ is the gravitational acceleration and the operator H is given by eq. (B5). Substitution of eqs (13) through (16) into (12) shows that eq. (12) is an integral equation for u because the total displacement field u also appears in the right-hand-side of (15)–(16).

We account for lateral variations in the isotropic elastic parameters {μ, κ}, density ρ, and viscosity η. By substituting (10) into (15)–(16), (15)–(16) into (14), and employing the expressions (A1)–(A4) for the excitation functions in (14), we arrive at the following expressions for :

(17)

(18)

The interaction kernels δω_j′j are given in. In (17)–(18), the integration is carried out over all points on the unit sphere; ∂_ν and ∂_α are partial derivatives acting on the primed coordinates with respect to variables ν and α which measure angular distance towards due East and South, respectively, from the argument of the considered function; γ is the azimuth of the minor arc from scattering point to observation point measured positive counter-clockwise from due South; θ is the angular distance along the minor arc from to. The functions a (l), b (l), and c (l) are defined by eqs (A5)–(A7).

The application of seismic scattering theory to post-seismic relaxation on a laterally heterogeneous viscoelastic structure problem bears many similarities to its application to seismic wave propagation on a laterally heterogeneous earth model. It can be verified that the form of perturbed scattered displacement field in the viscoelastic problem, prescribed by eqs (13) and (17)–(18), is identical to the form of perturbed scattered displacement in the seismic wave propagation problem prescribed by eq. (34) of Friederich (1999). The only difference in the two sets of results are in the various constant factors that appear ahead of the scattering integrals—factors which are particular to the different problems being treated. This exact correspondence of the structure of the results in the two treatments reflects the fact that the post-seismic relaxation problem being treated in the present paper is based on the same equations that govern the seismic wave propagation problem, with the understanding that the post-seismic problem is evaluated in the zero frequency limit because of the neglect of inertial terms.

Taking necessary first and second order spatial derivatives of the potential functions in eqs (11) and (17), (18), we may define an integral equation for the resulting state vector of the form

(19)

where A is a matrix with scalar elements and

(20)

In (20), the index j spans the entire spheroidal and toroidal mode basis set. The aspherical potential state vectors given by the ensemble of and provide a complete description of the 3-D deformation field. Furthermore these state vectors contain only the minimum amount of information necessary to describe the deformation field since the indicated spatial derivatives are necessary in order to evaluate the right-hand sides of eqs (17) and (18).

3 Numerical Implementation

3.1 Iterative Solution

It is necessary to build a solution to (19) which involves known quantities as well as discretize the problem. Let a laterally heterogeneous volume be described by an area Γ on the unit sphere extending to a certain depth. We assume that Γ may be represented by N discrete cells. Let be the locations of the cells on the unit sphere, and Δ²_k the corresponding area of the kth cell. Similarly, let {r_d, d= 1, …, D} be the locations of the depth points with corresponding depth increment Δ_d. Since both Φ⁰ and A are known at these discretized locations, eq. (19) could, in principle, be solved by inversion for and the corresponding displacement field then obtained by eq. (10). Such an inversion is impractical since Φ may contain millions of unknowns. A better procedure which we shall adopt involves iterative evaluation of eqs (17) and (18). There are two possibilities for this evaluation. First we define the state vectors

(21)

The discretized form of eq. (19) is

(22)

An iterative solution to (22) applicable for sufficiently weak heterogeneity (small δκ, δμ(s), δρ and/or small heterogeneous volume) is

(23)

This can be solved iteratively by defining B⁽⁰⁾(s) =E⁰(s) and

(24)

In the general case where aspherical structure may be strong one may employ an iterative solution. The solution to (22) is constructed by defining B⁽⁰⁾(s) =E⁰(s) and for β > 0

(25)

where I is the identity matrix. Then the solution of (22) is

(26)

Calculation of the product KB⁽ⁿ⁾ in (24) involves essentially evaluation of the right-hand sides of (17) and (18) with elements of B⁽ⁿ⁾ replacing the various potential functions. Since B⁽⁰⁾ and K are known, eqs (25) and (26) define an iterative solution to (22) in terms of known quantities. Which solution method—eq. (23) or (26)—is appropriate depends upon the strength of the heterogeneity, the length scales of the relaxation, and size of the heterogeneous region. In the examples to be considered in, we use the solution (23) for weak heterogeneity.

3.2 Evaluation Of Scattering Integrals

The next issue concerns the calculation of the matrix product K B⁽ⁿ⁾ in (23). As remarked above this involves evaluation of Φ¹_j using the scattering integrals on the right-hand sides of (17) and (18). For example, at first sight an evaluation of for a single and single mode j requires ∼N×D×M computations, where M is the total number of modes in the basis set. To evaluate for all and all modes would apparently require N²×D×M² computations, which could be extremely large.

Efficiency is gained in two ways. First, we note that each interaction kernel δω⁽ⁿ⁾_j′j is a sum of separable functions of the mode indices j′ and j, respectively. For the effects of δμ—heterogeneity we have

(27)

for known functions f and g, and additional terms with similar structure are added for the effects of δκ and δρ—heterogeneities. Using azimuthally symmetric deformation as an example, define

(28)

Then the azimuthally symmetric part of the aspherical spheroidal potential can be written in the form

(29)

where θ is the angular distance between and. Similar expressions can be written for the azimuthally asymmetric parts of Φ¹_S and Φ¹_T. Use of the function defined in eq. (28) results in a reduction in the numerical effort to evaluate the aspherical potentials by a factor of M compared with the direct evaluation of (17) or (18).

Second, since each Legendre function X^m_l satisfies a Helmholtz equation, direct evaluation of the Laplacian of the aspherical potential functions in eqs (17) and (18) shows that each aspherical potential satisfies a Helmholtz equation:

(30)

Suppose that we knew the values of Φ¹_S and Φ¹_T and their lateral gradients on a closed contour C enclosing Γ on the unit sphere (Fig. 2). Then it can be shown (e.g. Dahlen 1980) that the values of Φ¹_S and Φ¹_T elsewhere are given by Green's theorem on the unit sphere:

(31)

where satisfies

(32)

where is the Dirac δ-function. Note that points outward or inward according to whether is inside or outside C (Fig. 2). It should be noted that Dahlen's (1980)eq. (33) corresponding to our eq. (31) is in error by a factor of −1.

2

Geometry for the representation theorem eq. (31). The contour C is taken to enclose the scattering region Γ. Values of aspherical potentials may be evaluated either within or outside of C, with local normal to pointing outward or inward, respectively.

Open in new tabDownload slide

3.3 Procedure For Evaluating Aspherical Potentials

The considerations of the previous two sections suggest the following strategy for evaluating the aspherical potentials.

• Define a suitable contour C enclosing Γ on the unit sphere, and represent this contour with points with corresponding segment lengths Δ_c.

• Evaluate and for and their corresponding first and second spatial derivatives (i.e. ∂_ν, ∂_α, ∇²₁, ∂²_να, and ∂²_αα) using the discretized sum over Γ, e.g. eq. (22).

• 33
  Evaluate and its first and second spatial derivatives for using discretized forms of eqs (31) and (32), i.e.

  

  (33)

  Note that this procedure applies more generally to the evaluation of the matrix product KB⁽ⁿ⁾ when iterating (23) to calculate B⁽ⁿ⁺¹⁾.

The total numerical effort to calculate Φ¹_S and Φ¹_T is a sum of the efforts to evaluate (28) and to accomplish steps 2 and 3 above. The evaluation of (28) for all and all r_d requires ∼N×D×M computations. The evaluation of step 2 for all modes also requires ∼N×D×M computations. The evaluation of step 3 for all modes and all requires ∼N× C ×M computations. Thus the total computation time scales like

(34)

If for a given minimum wavelength scale we assume D∼N^1/2 and C ∼N^1/2, then

(35)

That is, for a given minimum wavelength scale the numerical effort scales with the volume of the heterogeneous region. The coupled mode method is able to achieve this level of efficiency essentially because of the representation theorem (31) that applies to each component of the viscoelastic mode sum.

3.4 Time Domain Evaluation

Until this point our treatment has been inclusive of both ‘purely elastic’ deformation and purely viscoelastic deformation. Strictly speaking, for an impulsive seismic source these are represented by aspherical potentials Φ_j which vary according to ∼1/s or a more complicated s-dependence, respectively. In the time domain these yield displacements that behave with time as the Heaviside step function at time t= 0 or a time-varying signal, respectively. The purely elastic deformation is a by-product of the static modes. As remarked in, on the laterally homogeneous Earth the ‘static’ modes are associated with only purely elastic deformation, but on the laterally heterogeneous Earth they are associated with both elastic and viscoelastic deformation because of the effect of viscoelastic mode → static mode coupling as well as static mode → static mode coupling. The aspherical static mode potentials also contain the purely elastic displacement field contributed by both the laterally homogeneous static potentials and the effects of static mode → static mode coupling. In this paper we are interested primarily in the viscoelastic deformation field, i.e. the evolution of the displacement field at times t > 0. We must therefore reckon with both a poorly approximated purely elastic deformation component and a well-approximated viscoelastic component when obtaining time domain results. Fortunately the estimation of these two separate components are practically independent using the numerical inverse Laplace transform described below. This is very convenient as only a limited number of static modes will be employed, and they cannot adequately represent the short-wavelength static deformation field, but they will serve the purpose of representing the long wavelength effects of viscoelastic mode → static mode and static mode → viscoelastic mode coupling.

In order to obtain deformation in the time domain, it is necessary to evaluate the inverse Laplace transform of u (r ; s) or, equivalently, the state vector E(s). Eq. (23) again is

(36)

Let the zeros of det(I−K (s)) define a set E₁. These are then the poles of [I−K(s)]⁻¹. Except in degenerate cases these are simple poles. It can be shown for a spherically symmetric model that the poles must lie on the negative real s-axis (e.g. Pollitz 1992). I assume without proof that this is the case for a general 3-D model. From eq. (6)E⁰(s) has simple poles at E₂={−s_S}∪{−s_T} for all spheroidal and toroidal modes on the laterally homogeneous model, also on the negative real s-axis. For an impulsive source, M (s) has a pole at s= 0, thus E⁰(s) has an additional pole at s= 0. These results allow us to represent E(s) as

(37)

The s_n are bounded by the inverse Maxwell relaxation time of the weakest volume occupying the earth model; let the largest pole of E(s) be denoted −s_max. (Note that any choice of s_max greater than this value is also satisfactory from the point of view of the following technique, which is designed to account for any temporal signal with inverse decay time up to the largest pole of E(s).)

A collocation technique is used to represent displacement or strain components of E(s) in terms of a small number of poles. To represent a function F(s) which possesses simple poles between s= 0 and s=−s_max we employ an expansion using 7 inverse relaxation times

(38)

with τ_j taking on the values 2 s_max, s_max, s_max/2, s_max/5, s_max/10, s_max/100, and s_max/500. In the time domain

(39)

where H (t) is the Heaviside function. In order to estimate the weighting coefficients B_j I utilize 12 evaluations of F(s) in the complex s-plane. I have found that accurate reconstruction of F(s) can be performed by taking samples of the function at the sample points in the complex s-plane shown in Fig. 3(a). The B_j are estimated from these 12 evaluations by least squares inversion. The stability of this inversion is very sensitive to the choice of sample points. In experiments with different configurations of sample points, I found that the condition number of the inversion for the B_j is generally high, but it is lowest for the configuration shown in Fig. 3(a). A small amount of damping is introduced to scale up the smallest eigenvalues when the system is solved with singular value decomposition.

To demonstrate the accuracy of the representations (38) and (39), consider a single component of relaxation, i.e.

(40)

ɛ₁(t) is the exact inverse Laplace transform of E₁(s). In Fig. 3(b)ɛ₁(t) is compared with the numerical approximation f(t) for a range of s₁ within −2s_max≤s₁ < 0. The agreement between the exact and numerical approximation is excellent. There is a negligible elastic offset B₀H(t) found in the numerical approximation even though such a component is included in the parametrization of the numerical inverse transform (eq. 38). Needless to say, any deformation component equal to a constant multiple of H(t) will be estimated exactly by the inverse transform. Thus, the estimation of time-dependent deformation is practically decoupled from the estimation of the elastic offset. The reconstructive power of the inverse Laplace transform breaks down for times t > ∼ 150/s_max, i.e, about 150 Maxwell relaxation times. A greater number of sample points and collocation points would be needed to achieve accurate results at longer times.

3

(a) Set of 12 collocation points on the negative real s-axis and the set of sample points in the first quadrant of the complex s-plane used to estimate the expansion coefficients of F(s) (eq. 38). (b) Comparison of exact (eq. 40) and numerically estimated (eq. 39) inverse Laplace transform of a single relaxation component with inverse relaxation time s₁.

Open in new tabDownload slide

4 Relaxation On A 3-D Structure

4.1 Model Description

A synthetic viscoelastic model is constructed in order to demonstrate the effect of a shallow heterogeneous structure upon the relaxation pattern. The 1-D background model is shown in Fig. 4(a). It is characterized by constant elastic parameters λ and μ throughout the volume of the Earth, and a three-layer viscoelastic structure consisting of an elastic upper crust (0–16 km depth), viscoelastic lower crust (16–30 km) of viscosity 10¹⁹ Pa s, and viscoelastic mantle (>30 km) of viscosity 3 × 10¹⁸ Pa s. The 3-D structure consists of the background model plus a 100 × 50 × 65 km³ heterogeneous volume of uniform viscosity η_h (Figs 4b and c) occupying the depth range 1–80 km. We shall explore the relaxation patterns of the 1-D and 3-D models driven by seismic sources confined to the upper crust. In this section we shall consider only the non-gravitational case (g-effects are neglected).

4

(a) Viscoelastic structure on the 1-D background earth model, consisting of 3 layers with uniform elastic parameters but varying viscosity. (b), (c) Viscoelastic structure on the laterally heterogeneous earth model. It is the 1-D structure modified to have a 100 × 50 × 65 km³ heterogeneous volume of viscosity η_h. (The spherical boundaries of the heterogeneous region appear rectangular because of the very small sphericity at this scale.)

Open in new tabDownload slide

The viscoelastic mode decay times for the 1-D earth model are shown in Fig. 5. Each point represents the wavelength (=2πR/(l+ 1/2)) and relaxation time (=1/s_j) of the jth viscoelastic mode. To represent post-seismic deformation down to 16 km wavelength there are 17 500 separate modes which collapse onto 7 spheroidal and 2 toroidal dispersion branches. To the maximum l employed, this set of viscoelastic modes is complete for the given 1-D earth model. In general, the fundamental mode with the longest relaxation times tends to dominate post-seismic deformation near the Earth's surface on the 1-D model, but all dispersion branches become important when coupling of modes at depth on the 3-D model is factored in. The corresponding normalized eigenfunctions of the modes for a horizontal wavelength of 30 km (spherical harmonic degree 1333) are displayed in Fig. 6(a).

5

(a) Spheroidal and (b) toroidal dispersion of viscoelastic normal modes on the 1-D earth model. Wavelength refers to 2πR/(1 + 1/2), where R is Earth's radius and l is spherical harmonic degree number.

Open in new tabDownload slide

6

Displacement eigenfunctions of (a) all 7 spheroidal and 2 toroidal viscoelastic modes and (b) the first 5 static modes (free oscillations) on a spherically stratified viscoelastic medium (described in ).

Open in new tabDownload slide

In order to concentrate the static mode energy near the volume of interest with as few modes as possible, they are calculated with a zero-displacement boundary condition at a depth of 5 wavelengths. The eigenfunctions of the first 5 static modes (i.e. free oscillations) at a horizontal wavelength of 30 km are shown in Fig. 6(b). These modes capture a significant part of the total static displacement field in the depth range of faulting as well as the depth range of the laterally heterogeneous structure. As remarked previously, this (and any finite) set of static modes being employed is incomplete because, in principle, an infinite number of Earth's free oscillation dispersion branches are needed for a complete basis set.

The heterogeneous volume has been divided into 32 768 cells of dimension 32 × 32 × 32, i.e. N= 1024 and D= 32. For one seismic source, the computations to be described required 18 h on an Ultra 10 single processor, taking the ascending series solution (eq. 23) up to third order in the heterogeneity. Larger problems are certainly feasible: computer memory is the limiting factor. The mode basis set and model discretization chosen here are the largest that could be achieved with 512 Mb read-access memory.

4.2 Examples

Post-seismic deformation is calculated using two different fault geometries. The sources are: a right-lateral strike-slip East–West trending fault of length 50 km, width 10 km, extending from 5 to 15 km depth; a thrust fault on a 50 km long, East-West trending plane dipping 45° towards the South. In both cases slip U= 1 meter is assigned. Comparison of results on the 3-D and 1-D models enables us to isolate the effect of the 3-D heterogeneity. No attempt will be made here to validate the calculated solution on the 3-D model. This awaits future comparison with independent calculations (e.g. finite element model).

All calculations of the response to finite faulting are realized by discretizing the fault plane into a sufficient number of small patches. The response to an effective point source contributed by each patch is evaluated, then summed to obtain the potential functions Φ⁰_j on the laterally homogeneous earth model. The algorithm presented in is then used to build the solution on the laterally heterogeneous model.

Practically all of the computational effort lies in calculating the aspherical potentials Φ¹_j over the surface projection of the scattering volume Γ; in this case, the area enclosed by points b, c, f, and e (Fig. 4). We may then readily calculate displacements and strains at arbitrary points in the Earth. In the following figures I present vector displacement and tensor strain components at 10 km depth over a 150 × 150 km² area centered on the 3-D heterogeneity.

Figs 7–12 present results calculated on a 3-D model with a weak heterogeneous region of viscosity η_h= 10¹⁸ Pa s. Figs 7 and 8(b) show right-lateral shear strain on vertical East-West trending planes and vertical displacement, respectively, evaluated several times after a synthetic strike-slip event. Fig. 8(a) shows the corresponding horizontal vector displacement at time t= 20 yr after the event. The presence of the low viscosity zone amplifies both horizontal displacements and strains by up to 33 per cent compared with the 1-D model. Given the much lower viscosity of the heterogeneous region, larger amplifications might be expected. The width of the weak region, however, is only 50 km or about 3 times the elastic plate thickness. The edges of this fairly narrow weak zone thus retard potentially strong horizontal amplification effects. The vertical motions predicted by the 3-D model, however, are about 100 per cent larger than those on the 1-D model (Fig. 8b). This arises from the great sensitivity of post-seismic vertical motions from a strike-slip faulting event to the depth-dependent viscosity structure. Pollitz et al. (2000) found that post-seismic vertical motions following strike-slip faulting on a vertical fault will generally be in the opposite sense to the coseismic motions if the ratio of lower crust to mantle Maxwell relaxation times is greater than about 4.5. In that situation one would predict post-seismic subsidence in those quadrants where coseismic uplift occurred. In the present case, on the 1-D model this ratio is 3.3, and a small amount of post-seismic uplift occurs in those quadrants (Fig. 8). On the 3-D model this ratio is 1.0, and the amount of post-seismic uplift during the early post-seismic epoch is greatly amplified as a result.

7

Strain component e_xy (x, y denotes distance measured towards due East and North, respectively) evaluated at depth 10 km and at several times following an impulsive seismic event on a synthetic strike-slip fault which extends from depth 5–15 km and has right-lateral slip U= 1 metre. The surface projection of the fault plane is indicated. The points b, c, f, and e identify the projection to the 10 km, depth isosurface of the 3-D heterogeneous volume of viscosity η_h= 10¹⁸ Pa s (see Fig. 4).

Open in new tabDownload slide

8

(a) Horizontal vector displacement evaluated at depth 10 km and time t= 20 yr following a strike-slip event. Viscosity in the 3-D heterogeneous volume is η_h= 10¹⁸ Pa s (see Fig. 4). The fault plane is marked as a gray rectangle. Although 1024 vectors are available on the 10 km depth isosurface, only every 25th point is plotted. The contour used in employing the representation theorem (eq. 33) is shown as a dashed circle. (b) Vertical displacement evaluated at depth 10 km and at several times following a strike-slip event. The vertical projection of the fault is indicated with a black line segment.

Open in new tabDownload slide

9

Strain component e_yy (y denotes distance measured towards North) evaluated at depth 10 km and at several times following an impulsive seismic event on a synthetic 45°—dipping thrust fault which extends from depth 5–15 km and has right-lateral slip U= 1 meter. The surface projection of the fault plane is indicated. The points b, c, f, and e identify the projection to the 10 km depth isosurface of the 3-D heterogeneous volume of viscosity η_h= 10¹⁸ Pa s (see Fig. 4).

Open in new tabDownload slide

10

(a) Horizontal vector displacement evaluated at depth 10 km and time t= 20 yr following a thrust event. Viscosity in the 3-D heterogeneous volume is η_h= 10¹⁸ Pa s (see Fig. 4). The fault plane is marked as a gray rectangle. (b) Vertical displacement evaluated at depth 10 km and at several times following a thrust event. The vertical projection of the fault is indicated with a black rectangular plane.

Open in new tabDownload slide

11

Same as Fig. 9, but thrust fault is located 15 km further to the north.

Open in new tabDownload slide

12

Same as Fig. 10, but thrust fault is located 15 km further to the north.

Open in new tabDownload slide

Figs 9 and 10(b) show cumulative North-South uniaxial strain and vertical displacement, respectively, at different times following a thrust-faulting event on a 45° south-dipping, East-West trending plane. Fig. 10(a) shows the corresponding horizontal vector displacement at time t= 20 yr after the event. In this case the effect of the 3-D aspherical region on the strain pattern is very pronounced, affecting not only the amplitude but also the spatial pattern. While in the 1-D case the central post-seismic region of compression is surrounded by regions of moderate post-seismic tension, these tension lobes are much stronger on the 3-D model. The differences between the 1-D and 3-D cases extend slightly beyond the heterogeneous region to just south of, where the 3-D model exhibits larger tensional strains than the 1-D model. (This is also seen from the horizontal displacement gradients just outside the aspherical region in Fig. 10.) The effect of the aspherical region on the vertical displacement pattern (Fig. 10) is more complicated. Within the aspherical region, the 1-D model produces post-seismic subsidence in the near-fault region and gentle post-seismic uplift elsewhere, but the 3-D model produces post-seismic uplift almost everywhere. Outside the aspherical region, both models predict subsidence, the greater amplitude being produced by the 3-D model. These patterns reflect the competing tendencies of accelerated relaxation on the 3-D model versus the ratio of crust-to-mantle Maxwell relaxation time, which equals 3.3 on the 1-D model and 1.0 on the 3-D model.

A final example is shown in Figs 11 and 12, which is identical to the previous thrust fault example except that the fault is located 15 km further north, just on the edge of the aspherical region. In the 3-D case there are now pronounced asymmetries in the overall deformation pattern. This is demonstrated for the horizontal displacement (Fig. 12a), but it is true also for the vertical displacement and horizontal strain patterns. There are no greatly amplified southward horizontal displacements in the region just north of the fault in Fig. 12(a), compared with the large southward displacements seen in Fig. 10(a); this reflects the different position of the thrust fault relative to the edge of the heterogeneous region in the two cases. In the 3-D case, the amplitude of the uniaxial strain (Fig. 11) and horizontal displacement patterns are greater in this example than in the example of Figs 9 and 10 in which the thrust fault was placed in the middle of the aspherical region. This is seen by the larger tensional strain lobe just north of the boundary as well as larger horizontal displacements south of the fault. This is caused by the diminished boundary effect of the background medium in this example: is 40 km south of the thrust fault in this example compared with 25 km south in the previous example, leading to more vigorous relaxation of the low-viscosity region south of the fault.

4.3 Convergence Of Born Series

The convergence of the Born series in eq. (23), defined by the recursion formula (24), depends on several factors including the size of the heterogeneous region, the contrast in viscoelastic structure with respect to the background medium, and the s-value (Laplace transform parameter) at which the heterogeneous potentials are to be evaluated. In general, taking a larger heterogeneous region, a larger contrast, or an s-value closer to the negative real s-axis will produce a Born series sum that tends towards divergence (unstable iteration of eq. 24). This tendency is clear by analogy with well-known examples in seismic wave propagation. A test of the convergence of the Born series in the first example (post-seismic response to strike-slip faulting) is presented in Fig. 13. It demonstrates that truncating the Born series at 3rd order is sufficiently accurate for the given aspherical viscoelastic structure at all times considered.

13

(a) Profiles of cumulative East-displacement following an impulsive strike-slip event (as in Figs 7 and 8) as a function of elapsed time at 10 km depth along a North–South profile which bisects the fault. Different curves on each snapshot represent the result of truncating the Born series (eq. 23) with different number of terms, i.e. from first to fourth order scattering. The solid curve represents the laterally homogeneous model. Numerals correspond to the order of scattering of each respective curve. (b) Cumulative East-displacement relative to that of the first-order scattering displacement.

Open in new tabDownload slide

In general, one needs to bear in mind the convergence of the Born series when selecting the window of valid time domain results. From this corresponds to times 0 < t < ∼150/s_max. As remarked in, the value of s_max to be used is arbitrary as long as it is larger than the largest inverse material relaxation time of the viscoelastic earth model. A larger value will result in a smaller time window of valid results; on the other hand, a larger value will require samples of the potential functions at s-values further removed from the negative s-axis (i.e. Fig. 3a). I find in numerical experiments that this always results in a more stable Born series summation, and in many situations it may be better to obtain a stable iteration at the expense of a narrower valid time window. (Other factors being equal, one can always increase the temporal range of valid results by taking more samples of the deformation field in the Laplace transform domain.)

5 Conclusions

The problem of post-seismic relaxation has been approached following an impulsive seismic event on an aspherical viscoelastic model using the theory of coupled viscoelastic normal modes. It is a perturbation treatment with an iterative solution. This approach provides a powerful computational framework for post-seismic relaxation problems. The numerical effort involved with accounting for laterally heterogeneous viscoelastic structure is proportional to the volume of the heterogeneous region. The problem is cast as a set of coupled integral equations in aspherical spheroidal and toroidal potential functions defined on the surface projection of the laterally heterogeneous volume onto the unit sphere. The approximation involved with the solution are the needs to (1) discretize the problem, (2) truncate the modal basis set, (3) iterate the solution of the discretized integral equations using, as in the examples presented here, an ascending series solution beginning with the solution on the 1-D viscoelastic model, and (4) perform a numerical inverse Laplace transform. The method of obtaining time domain results with the inverse Laplace transform, valid to a certain maximum time, requires a minimal amount of evaluations in the Laplace transform domain. The examples given here involve laterally heterogeneous domains of dimension 100 × 50 × 65 km³ and evaluation times up to 150 Maxwell relaxation times (of the weakest component of the earth model). The domain size used here is modest, with 32 768 cells to represent the laterally heterogeneous volume. The calculations described here required about 18 h on an Ultra 10 single processor.

It should be noted that the employed method is semi-analytic in both space and time. This carries several advantages over purely numerical treatments. A single computation yields time-dependent deformation at all times up to a certain maximum time (which depends on the amount of sampling taken in the complex Laplace transform domain) over the continuum of time after the forcing event, not at discrete samples of time. Similarly, it yields deformation on the entire spherical Earth, not only within the laterally heterogeneous volume. This advantage would carry greatest weight if deformation were required over a large region which could be characterized predominantly by a 1-D background model and only in part by a 3-D model.

Acknowledgments

I thank Marleen Nyst and Jim Savage for their comments on a preliminary draft. This paper benefitted from the constructive criticisms of two anonymous reviewers and the Associate Editor.

Appendix

Appendix A:

Given here are the explicit expressions for the excitation tensor and in terms of the source-observation point geometry of Fig. 1. These expressions describe the normal mode response to seismic source excitation derived by Gilbert & Dziewonski (1975). Anticipating evaluation at scattering points within Γ, the argument of the observation point is taken at. For a spheroidal mode,

(A1)

(A2)

For a toroidal mode,

(A3)

(A4)

In these expressions the X^m_l are Legendre functions, Δ is the angular distance between and, and φ is the azimuth of the minor arc from r₀ to measured positive counter-clockwise from due South. The constants a (l), b (l), and c (l) are

(A5)

(A6)

(A7)

Appendix

Appendix B:

The following derivation closely parallels the derivation given by Hudson & Heritage (1981) for scattering in elastodynamic problems. Let u(r; s) be a quasi-static displacement field and T(r; s) the corresponding stress tensor in the Laplace transform domain with independent variable s. We shall assume for simplicity that the earth model consists of welded elastic and viscoelastic layers between some lower boundary and Earth's surface. Following Pollitz (1997), we consider deformation prescribed by viscoelastic-gravitational coupling taking account of coupling of viscoelastic deformation with Earth's gravitational acceleration, but neglecting coupling with perturbations in gravitational potential. This approximation is sometimes referred to in the seismological literature as Cowling's approximation. Dropping the notation s, the equations of static equilibrium for a distributed body force f(r) are:

(B1)

(B2)

(B3)

(B4)

(B5)

where G is the gravitational constant, ρ₀ and g₀ are the density and gravitational acceleration on a reference spherically symmetric earth model, and c is the elastic tensor. Eqs (B1)–(B5) are to be solved subject to continuity of displacement and traction and the free-surface boundary condition

(B6)

and an appropriate lower boundary condition in Earth's interior. The following development is valid for the case of this boundary condition being either vanishing traction or vanishing displacement on a fixed interior surface. In general, one may define Ω_u and Ω_T (with local unit normal vector ) to be the surfaces on which boundary conditions u (r) = 0 and are applied, respectively.

We may decompose the following quantities into spherically symmetric and asymmetric components:

(B7)

(B8)

(B9)

(B10)

where u⁰ and T⁰ are defined by (B1)–(B5) and associated boundary conditions with (B1)–(B2) replaced by

(B11)

(B12)

with u⁰(r) = 0 for r∈Ω_u and for r∈Ω_T. Eqs (B1), (B5), (B7) and (B9)–(B10) give

(B13)

The aspherical part of T is provided by eqs (B2), (B7)–(B9) and (B12):

(B14)

Subtracting (B11) from (B13) and making use of (B2), (B8)–(B10) and (B14) yields

(B15)

(B16)

where the forcing term is given by

(B17)

with u¹(r) = 0 for r∈Ω_u and for r∈Ω_T. It follows immediately from (B7) that the boundary conditions on Ω_u for u⁰ and u¹ are compatible with the corresponding boundary condition for the total displacement field u. It can be shown that the boundary condition ) implies that on this surface. From (B14) this leads to, and together with and (B9), this shows that the boundary condition is satisfied by the solutions for u⁰, T⁰ and u¹, T¹.

Consider the contribution over a volume δV centered on a point r′. The force f¹ can be described as the sum of an equivalent moment tensor and single force:

(B18)

where

(B19)

(B20)

If the solution to (B11)–(B12) has the form of eq. (6) (), then (B15)–(B16) show that the solution for u¹, the aspherical component of u, also has the form of eq. (4) with local equivalent moment tensor and single force given by (B19)–(B20). Convolving over all equivalent forces distributed over the volume of the Earth and defining δM=mδV, δF=fδV, yields eq. (14).

Appendix

Appendix C:

Interaction kernels δω_j′j which appear in the scattering integrals (17) and (18) may be given in terms of the spheroidal mode displacement eigenfunctions U(r; s) and V(r; s) and toroidal mode displacement eigenfunction W(r; s), associated with spherical harmonic degrees l_S and l_T, respectively. Dropping the argument (r ; s) we define

(C1)

where the overdot denotes differentiation with respect to r. Let κ(r), μ(r), ρ(r), and η(r) be the elastic parameters, density, and viscosity on the laterally heterogeneous model. Then for a Maxwell viscoelastic fluid, relaxation in shear is represented by the transformed shear modulus

(C2)

Let δκ(r ; s) =κ(r) −κ₀(r) and δμ(r ; s) =μ(r; s) −μ₀(r; s) represent the lateral variations in the viscoelastic model, δρ(r) =ρ(r) −ρ₀(r) be the lateral variation in density, and g₀(r) be the gravitational acceleration. The interaction kernels take the form

(C3)

(C4)

(C5)

(C6)

(C7)

(C8)

(C9)

(C10)

(C11)

The integration is carried out from radius 0 to Earth's radius R.

References