Post-seismic rebound of a spherical Earth: new insights from the application of the Post–Widder inversion formula

The post-seismic response of a viscoelastic Earth to a seismic dislocation can be computed analytically within the framework of normal-modes, based on the application of propagator methods. This technique, widely documented in the literature, suffers from several shortcomings; the main drawback is related to the numerical solution of the secular equation, whose degree increases linearly with the number of viscoelastic layers so that only coarse-layered models are practically solvable. Recently, a viable alternative to the standard normal-mode approach, based on the Post–Widder Laplace inversion formula, has been proposed in the realm of postglacial rebound models. The main advantage of this method is to bypass the explicit solution of the secular equation, while retaining the analytical structure of the propagator formalism. At the same time, the numerical computation is much simplified so that additional features such as linear non-Maxwell rheologies can be simply implemented. In this work, for the first time, we apply the Post–Widder Laplace inversion formula to a post-seismic rebound model. We test the method against the standard normal-mode solution and we perform various benchmarks aimed to tune the algorithm and to optimize computation performance while ensuring the stability of the solution. As an application, we address the issue of finding the minimum number of layers with distinct elastic properties needed to accurately describe the post-seismic relaxation of a realistic Earth model. Finally, we demonstrate the potentialities of our code by modelling the post-seismic relaxation after the 2004 Sumatra–Andaman earthquake comparing results based upon Maxwell and Burgers rheologies.     

Analytical approach for the toroidal relaxation of viscoelastic earth

This paper is concerned with post-seismic toroidal deformation in a spherically symmetric, non-rotating, linear-viscoelastic, isotropic Maxwell earth model. Analytical expressions for characteristic relaxation times and relaxation strengths are found for viscoelastic toroidal deformation, associated with surface tangential stress, when there are two to five layers between the core–mantle boundary and Earth's surface. The multilayered models can include lithosphere, asthenosphere, upper and lower mantles and even low-viscosity ductile layer in the lithosphere. The analytical approach is self-consistent in that the Heaviside isostatic solution agrees with fluid limit. The analytical solution can be used for high-precision simulation of the toroidal relaxation in five-layer earths and the results can also be considered as a benchmark for numerical methods. Analytical solution gives only stable decaying modes—unstable mode, conjugate complex mode and modes of relevant poles with orders larger than 1, are all excluded, and the total number of modes is found to be just the number of viscoelastic layers between the core–mantle boundary and Earth's surface—however, any elastic layer between two viscoelastic layers is also counted. This confirms previous finding where numerical method (i.e. propagator matrix method) is used. We have studied the relaxation times of a lot of models and found the propagator matrix method to agree very well with those from analytical results. In addition, the asthenosphere and lithospheric ductile layer are found to have large effects on the amplitude of post-seismic deformation. This also confirms the findings of previous works.     

Coseismic and post-seismic signatures of the Sumatra 2004 December and 2005 March earthquakes in GRACE satellite gravity

The GRACE satellite mission has been measuring the Earth's gravity field and its temporal variations since 2002 April. Although these variations are mainly due to mass transfer within the geofluid envelops, they also result from mass displacements associated with phenomena including glacial isostatic adjustment and earthquakes. However, these last contributions are difficult to isolate because of the presence of noise and of geofluid signals, and because of GRACE's coarse spatial resolution (>400 km half-wavelength). In this paper, we show that a wavelet analysis on the sphere helps to retrieve earthquake signatures from GRACE geoid products. Using a wavelet analysis of GRACE geoids products, we show that the geoid variations caused by the 2004 December ( M _w= 9.2) and 2005 March ( M _w= 8.7) Sumatra earthquakes can be detected. At GRACE resolution, the 2004 December earthquake produced a strong coseismic decrease of the gravity field in the Andaman Sea, followed by relaxation in the area affected by both the Andaman 2004 and the Nias 2005 earthquakes. We find two characteristic timescales for the relaxation, with a fast variation occurring in the vicinity of the Central Andaman ridge. We discuss our coseismic observations in terms of density changes of crustal and upper-mantle rocks, and of the vertical displacements in the Andaman Sea. We interpret the post-seismic signal in terms of the viscoelastic response of the Earth's mantle. The transient component of the relaxation may indicate the presence of hot, viscous material beneath the active Central Andaman Basin.     

Post-seismic relaxation following the great 2004 Sumatra-Andaman earthquake on a compressible self-gravitating Earth

The   M _w γ 9.0  2004 December 26 Sumatra-Andaman and   M _w = 8.7  2005 March 28 Nias earthquakes, which collectively ruptured approximately 1800 km of the Andaman and Sunda subduction zones, are expected to be followed by vigorous viscoelastic relaxation involving both the upper and lower mantle. Because of these large spatial dimensions it is desirable to fully account for gravitational coupling effects in the relaxation process. We present a stable method of computing relaxation of a spherically-stratified, compressible and self-gravitating viscoelastic Earth following an impulsive moment release event. The solution is cast in terms of a spherical harmonic expansion of viscoelastic normal modes. For simple layered viscoelastic models, which include a low-viscosity oceanic asthenosphere, we predict substantial post-seismic effects over a region several 100s of km wide surrounding the eastern Indian Ocean. We compare observed GPS time-series from ten regional sites (mostly in Thailand and Indonesia), beginning in 2004 December, with synthetic time-series that include the coseismic and post-seismic effects of the 2004 December 26 and 2005 March 28 earthquakes. A viscosity structure involving a biviscous (Burgers body) rheology in the asthenosphere explains the pattern and amplitude of post-seismic offsets remarkably well.     

Post-seismic reloading and temporal clustering on a single fault

Geological studies show evidence for temporal clustering of large earthquakes on individual fault systems. Since post-seismic deformation due to the inelastic rheology of the lithosphere may result in a variable loading rate on a fault throughout the interseismic period, it is reasonable to expect that the rheology of the non-seismogenic lower crust and mantle lithosphere may play a role in controlling earthquake recurrence times. We study this phenomenon using a 2-D, finite element method continuum model of the lithosphere containing a single strike-slip fault. This model builds on a previous study using a 1-D spring-dashpot-slider analogue of a single fault system to study the role of Maxwell viscoelastic relaxation in producing non-periodic earthquakes. In our 2-D model, the seismogenic portion of the fault slips when a predetermined yield stress is exceeded; stress accumulated on the seismogenic fault is shed to the viscoelastic layers below and recycled back to the seismogenic fault through viscoelastic relaxation. We find that random variation of the fault yield stress from one earthquake to the next can cause the earthquake sequence to be clustered; the amount of clustering depends on a non-dimensional number, W , called the Wallace number defined as the standard deviation of the randomly varied fault yield stress divided by the effective viscosity of the system times the tectonic loading rate. A new clustering metric based on the bimodal distribution of interseismic intervals allows us to investigate clustering behaviour of systems over a wide range of model parameters and those with multiple viscoelastic layers. For models with   W ≥ 1  clustering increases with increasing W , while those with   W ≤ 1  are unclustered, or quasi-periodic.     

Viscoelastic relaxation and long-lasting after-slip following the 1997 Umbria-Marche (Central Italy) earthquakes

We combine Global Positioning System (GPS) measurements with forward modelling of viscoelastic relaxation and after-slip to study the post-seismic deformation of the 1997 Umbria-Marche (Central Apennines) moderate shallow earthquake sequence. Campaign GPS measurements spanning the time period 1999–2003 are depicting a clear post-seismic deformation signal. Our results favour a normal faulting rupture model where most of the slip is located in the lower part of the seismogenic upper crust, consistent with the rupture models obtained from the inversion of strong motion data. The preferred rheological model, obtained from viscoelastic relaxation modelling, consists of an elastic upper crust, underlain by a transition zone with a viscosity of 10¹⁸ Pa s, while the rheology of deeper layers is not relevant for the observed time-span. Shallow fault creep and after-slip at the base of the seismogenic upper crust are the first order processes behind the observed post-seismic deformation. The deep after-slip, below the fault zone at about 8 km depth, acting as a basal shear through localized time-dependent deformation, identifies a rheological discontinuity decoupling the seismogenic upper crust from the low-viscosity transition zone.     

Material versus isobaric internal boundaries in the Earth and their influence on postglacial rebound

Most previous earth models used to calculate viscoelastic relaxation after the removal of the Late Pleistocene ice loads implicitly assume that there is no exchange of mass across the mantle density discontinuities on periods of tens of thousands of years (the material boundary formulation). In the present study, simple incompressible models are used to determine the Earth's behaviour in the case where the density discontinuity remains at a constant pressure rather than deforming with the material (the isobaric boundary formulation). The calculation of the movement of the boundary is more rigorous than in earlier studies and uses the local incremental pressure calculated at the depth of the boundary and allows for the vertical deformation caused by the change in volume as material changes phase. It is shown that the buoyancy modes associated with the density discontinuities decrease in strength and increase in relaxation time analogous to what results when the density contrast is reduced. Also, two viscoelastic modes arise from an isobaric boundary, which is also predicted when there is a contrast in rigidity or viscosity across a material boundary. The difference in predicted radial deformation between the isobaric boundary model and the material boundary model is largest for long-wavelength loads for which the material incremental pressure at depth is largest. If the isobaric boundary model is appropriate for the treatment of the mineral phase changes in the mantle on glacial rebound timescales, then previous inferences of the deep-mantle to shallow-mantle viscosity ratio based on large-scale deformation (spherical harmonic degree < 10) of the Earth and including data from the early part of the glacio-isostatic uplift are too small.     

Stress transfer relations among the earthquakes that occurred in Kerman province, southern Iran since 1981

We explore the possible stress triggering relationship of the   M ≥ 6.4  earthquakes that occurred in Kerman Province, southern Iran since 1981. We calculated stress changes due to both coseismic sudden movement in the upper crust and the time-dependent viscous relaxation of the lower crust and/or upper mantle following the event. Four events of   M ≥ 6.4  occurred between 1981 and 2005, on and close to the Gowk fault, show a clear Coulomb stress load to failure relationship. The  2003 M = 6.5  Bam earthquake, however, which occurred approximately 95 km SW of the closest Gowk event, shows a very weak stress relation to preceding earthquakes. The coseismic static stress change at the hypocentre of the Bam earthquake is quite small (∼0.006 bars). The time-dependent post-seismic stress change could be 26 times larger or 7 times lower than that of coseismic static stress alone depending on the choice of viscoelastic crustal model and the effective coefficient of friction. Given the uncertainties in the viscoelastic earth models and the effective coefficient of friction, we cannot confidently conclude that the 2003 Bam event was brought closer to failure through coseismic or post-seismic stress loading. Interestingly, the southern Gowk segment with a similar strike to that of the Bam fault, experienced a stress load of up to 8.3 bars between 1981 and 2003, and is yet to have a damaging earthquake.     

Effect of the viscoelastic lithosphere on polar wander speed caused by the Late Pleistocene glacial cycles

Previous studies of the wander of the rotation pole associated with the Late Pleistocene glacial cycles indicate that the predicted polar wander speed is sensitive to the density jump at the 670 km discontinuity, the thickness of the elastic lithosphere, and the lower mantle viscosity. In particular, the M1 mode related to the density jump at 670 km depth has been shown to contribute a dominant portion of predicted polar wander speed for sufficiently small lower mantle viscosities. In this study, we examine the sensitivity of polar wander to variations in the viscosity of the viscoelastic lithosphere using simplified compressible Maxwell viscoelastic earth models. Model calculations for earth models with a viscoelastic lithosphere of finite viscosity indicate that the contribution of the M1 mode is similar to those associated with the density discontinuity at the core–mantle boundary (C0 mode) and the lithosphere (L0 mode). We speculate that this is due to the interaction between the M1 mode and the transient mode associated with the viscoelastic lithosphere, which reduces the magnitude of polar wander rates. Therefore, the M1 mode does not contribute a dominant portion of the predicted polar wander speed for earth models with a viscoelastic lithosphere of finite viscosity. In this case, predictions of polar wander speed as a function of lower mantle viscosity exhibit the qualitative form of an 'inverted parabola', as predicted for the J ˙₂ curve. We caution, however, that these results are obtained for simplified earth models, and the results for seismological earth models such as PREM may be complicated by the interaction between the M1 mode and the large set of transient modes.     

Stress diffusion following large earthquakes: a comparison between spherical and flat-earth models

A comparison is performed between predictions of coseismic and post-seismic deformations, based on two different earth models. The main purpose of this work is to clarify in a rigorous manner the role of sphericity in the assessment of both residual deformations and large-scale post-seismic motions. The two models considered in this study are commonly employed for predictions of post-seismic deformations. They exhibit the common feature of a layered viscoelastic structure but differ in their global geometrical features, the first being characterized by a flat-earth approximation and the second by a spherical geometry. The purpose of this study is twofold. First, we wish to assess the influence of sphericity in predictions of post-seismic deformations, especially in the far field of major earthquakes. Second, we wish to elucidate the potential role of stress migration in earthquake triggering over great distances.     

Guided waves in three-dimensional structures

A new formulation for the propagation of surface waves in three-dimensionally varying media is developed in terms of modal interactions. A variety of assumptions can be made about the nature of the modal field: a single set of reference modes, a set of local modes for the structure beneath a point, or a set of local modes for a laterally varying reference structure. Each modal contribution is represented locally as a spectrum of plane waves propagating in different directions in the horizontal plane. The influence of 3-D structure is included by allowing coupling between different modal branches and propagation directions. For anisotropic models, with allowance for attenuation, the treatment leads to a set of coupled 2-D partial differential equations for the weight functions for different modal orders.
The representation of the guided wavefield requires the inclusion of a full set of modes, so that, even for isotropic models, both Love and Rayleigh modes appear as different polarization states of the modal spectrum. The coupling equations describe the interaction between the different polarizations induced by the presence of the 3-D structure.
The level of lateral variation within the 3-D model is not required to be small. Horizontal refraction or reflection of the surface wavefield can be included by allowing for transfer between modes travelling in different directions. Approximate forms of the coupled equation system can be employed when the level of heterogeneity is small, for example the coupling between the fundamental mode and higher modes can often be neglected, or forward propagation can be emphasized by restricting the interaction to a limited band of plane waves covering the expected direction of propagation.     

Material versus local incompressibility and its influence onglacial‐isostatic adjustment

We present an analytical form of the layer propagator matrix for the response of a locally incompressible, layered, linear‐viscoelastic sphere to an external load assuming that the initial density stratification ϱ ₀( r ) within each layer is parametrized by Darwin's law. From this, we show that the relaxation of a sphere consisting of locally incompressible layers is governed by a discrete set of viscous modes. The explicit dependence of the layer propagator matrix on the Laplace transform variable allows us to determine the amplitudes of the viscous modes analytically. Employing Darwin's parametrization, we construct three simplified earth models with different initial density gradients that are used to compare the effects of the local incompressibility constraint, div ( ϱ ₀ u )=0, and the material incompressibility constraint, div  u =0, on viscoelastic relaxation. We show that a locally incompressible earth model relaxes faster than a materially incompressible model. This is a consequence of the fact that the perturbations of the initial density are zero during viscoelastic relaxation of a locally incompressible medium, so that there are no internal buoyancy forces associated with the continuous radial density gradients, only the buoyancy forces generated by internal density discontinuities. On the other hand, slowly decaying internal buoyancy forces in a materially incompressible earth model cause it to reach the hydrostatic equilibrium after a considerably longer time than a locally incompressible model. It is important to note that the approximation of local incompressibility provides a solution for a compressible earth model that is superior to the conventional solutions for a compressible earth with homogeneous layers because it is based on an initial state that is consistent with the assumption of compressibility.     

Some remarks about the degree-one deformation of the Earth

The degree-one deformation of the Earth (and the induced discrepancy between the figure centre and the mass centre of the Earth) is computed using a theoretical approach (Love numbers formalism) at short timescales (where the Earth has an elastic behaviour) as well as at long timescales (where the Earth has a viscoelastic or quasi-fluid behaviour). For a Maxwell model of rheology, the degree-one relaxation modes associated with the viscoelastic Love numbers have been investigated: the M_o mode does not exist and there is only one transition mode (instead of two) generated by a viscosity discontinuity.
The translations at each interface of the incompressible layers of the earth model [surface, 670 km depth discontinuity, core-mantle boundary (CMB) and inner-core boundary (ICB)] are computed. They are elastic with an order of magnitude of about 1 mm when the excitation source is the atmospheric continental loading or a magnetic pressure acting at the CMB. They are viscoelastic when the earth is submitted to Pleistocene deglaciation, with an order of magnitude of about 1 m. In a quasi-fluid approximation (Newtonian fluid) because of the mantle density heterogeneity their order of magnitude is about 100 m (except for the ICB, which is in quasi-hydrostatic equilibrium at this timescale).     

Coseismic rotation changes from the 2004 Sumatra earthquake: the effects of Earth's compressibility versus earthquake induced topography

Polar motion is modelled for the large 2004 Sumatra earthquake via dislocation theory for an incompressible elastic earth model, where inertia perturbations are due to earthquake-triggered topography of density–contrast interfaces, and for a compressible model, where inertia perturbation due to compression-dilatation of Earth's material is included; density and elastic parameters are based on a multilayered reference Earth. Both models are based on analytical Green's functions, propagated from the centre to the Earth's surface. Preliminary and updated seismological solutions are considered in elucidating the effects of improving earthquake parameters on polar motion. The large Sumatra thrust earthquake was particularly efficient in driving polar motion since it was responsible for large material displacements occurring orthogonally to the strike of the earthquake and to the Earth's surface, as imaged by GRACE gravity anomalies over the earthquake area. The effects of earthquake-induced topography are four times larger than the effects of Earth's compressibility, for l = 2 geopotential components. For varying compressional Earth properties and seismic solution, modelled polar motion ranges from 8.6 to 9.4 cm in amplitude and between 117° and 130° east longitude in direction. The close relationship between polar motion direction, earthquake longitude and thrust nature of the event, are established in terms of basic physical concepts.     

Temporal and spatial variations of post-seismic deformation following the 1999 Chi-Chi, Taiwan earthquake

We use GPS displacements collected in the 15 months after the 1999 Chi-Chi, Taiwan earthquake  ( M _w 7.6)  to evaluate whether post-seismic deformation is better explained by afterslip or viscoelastic relaxation of the lower crust and upper mantle. We find that all viscoelastic models tested fail to fit the general features in the post-seismic GPS displacements, in contrast to the satisfactory fit obtained with afterslip models. We conclude that afterslip is the dominant mechanism in the 15-month period, and invert for the space–time distribution of afterslip, using the Extended Network Inversion Filter. Our results show high slip rates surrounding the region of greatest coseismic slip. The slip-rate distribution remains roughly stationary over the 15-month period. In contrast to the limited coseismic slip on the décollement, afterslip is prominent there. Maximum afterslip of 0.57 m occurs downdip and to the east of the hypocentral region. Afterslip at hypocentral depths is limited to the southern part of the main shock rupture, with little or no slip on the northern section where coseismic slip was greatest. Whether this results from along strike variations in frictional properties or dynamic conditions that locally favour stable sliding is not clear. In general, afterslip surrounds the area of greatest coseismic slip, consistent with post-seismic slip driven by the main shock stress change. The total accumulated geodetic afterslip moment is  3.8 × 10¹⁹ N m  , significantly more than the seismic moment released by aftershocks,  6.6 × 10¹⁸ N m  . Afterslip and aftershocks appear to have different temporal evolutions and some spatial correlations, suggesting that aftershock rates may not be completely controlled by the rate of afterslip.     

Gravitational viscoelastic relaxation of eccentrically nested spheres

We present a semi-analytical solution to the 2-D forward modelling of viscoelastic relaxation in a heterogeneous model consisting of eccentrically nested spheres. Several numerical methods for 2-D and 3-D viscoelastic relaxation modelling have been applied recently, including finite-element and spectral-finite-difference schemes. The present semi-analytical approach provides a model response against which more general numerical algorithms can be validated. The eccentrically nested sphere solution has been tested by comparing it with the analytical solutions for viscoelastic relaxation in a homogeneous sphere and in two concentrically nested spheres, and good agreement was obtained.     

Numerical models of mantle convection with secular cooling

A 2-D time-dependent finite-difference numerical model is used to investigate the thermal character and evolution of a convecting layer which is cooling as it convects. Two basic cooling modes are considered: in the first, both upper and lower boundaries are cooled at the same rate, while maintaining the same temperature difference across the layer; in the second, the lower boundary temperature decreases with time while the upper boundary temperature is fixed at 0°C. The first cooling mode simulates the effects of internal heating while the second simulates planetary cooling as mantle convection extracts heat from, and thereby cools, the Earth's core. The mathematical analogue between the effects of cooling and internal heating is verified for finite-amplitude convection. It is found that after an initial transient period the central core of a steady but vigorous convection cell cools at a constant rate which is governed by the rate of cooling of the boundaries and the viscosity structure of the layer. For upper-mantle models the transient stage lasts for about 30 per cent of the age of the Earth, while for the whole mantle it lasts for longer than the age of the Earth. Consequently, in our models the bulk cooling of the mantle lags behind the cooling of the core-mantle boundary. Models with temperature-dependent viscosity are found to cool in the same manner as models with depth-dependent viscosity; the rate of cooling is controlled primarily by the horizontally averaged variation of viscosity with depth. If the Earth's mantle cools in a similar fashion, secular cooling of the planet may be insensitive to lateral variations of viscosity.     

Spectral-finite-element approach to two-dimensional electromagnetic induction in a spherical earth

We present a spectral-finite-element approach to the 2-D forward problem for electromagnetic induction in a spherical earth. It represents an alternative to a variety of numerical methods for 2-D global electromagnetic modelling introduced recently (e.g. the perturbation expansion approach, the finite difference scheme). It may be used to estimate the effect of a possible axisymmetric structure of electrical conductivity of the mantle on surface observations, or it may serve as a tool for testing methods and codes for 3-D global electromagnetic modelling. The ultimate goal of these electromagnetic studies is to learn about the Earth's 3-D electrical structure.
Since the spectral-finite-element approach comes from the variational formulation, we formulate the 2-D electromagnetic induction problem in a variational sense. The boundary data used in this formulation consist of the horizontal components of the total magnetic intensity measured on the Earth's surface. In this the variational approach differs from other methods, which usually use spherical harmonic coefficients of external magnetic sources as input data. We verify the assumptions of the Lax-Milgram theorem and show that the variational solution exists and is unique. The spectral-finite-element approach then means that the problem is parametrized by spherical harmonics in the angular direction, whereas finite elements span the radial direction. The solution is searched for by the Galerkin method, which leads to the solving of a system of linear algebraic equations. The method and code have been tested for Everett & Schultz's (1995) model of two eccentrically nested spheres, and good agreement has been obtained.     

Inversion for laterally heterogeneous upper mantle S-wave velocity structure using iterative waveform inversion

In contrast to previous work, which treats the Earth's lateral heterogeneity as an infinitesimal perturbation to a spherically symmetrical starting model, we conduct iterative linearized waveform inversion for the Earth's laterally heterogeneous structure. We use the Direct Solution Method (DSM) (Geller et al. 1990a) to calculate synthetic seismograms and their partial derivatives for a laterally heterogeneous earth model. We invert surface-wave data from the IDA and GEOSCOPE networks. We expand the lateral heterogeneity of rigidity in spherical harmonics up to angular order number 8 and use three parameters to specify the depth dependence of each harmonic, giving us a total of 240 unknowns. The short-wavelength lateral heterogeneity (s = 4, 6 and 8) in the deeper part of the upper mantle obtained by this study differs significantly from M84A. The relative improvement in the variance reduction as compared with model M84A is about 20 per cent for the IDA data and more than 100 per cent for the GEOSCOPE data.     

Surface vertical displacements, potential perturbations andgravity changes of a viscoelastic earth model induced by internal point dislocations

Using the viscoelastic correspondence principle, we utilize the surface coseismic spheroidal deformation fields (i.e. vertical displacements, potential perturbations and gravity changes) of SNREI earth models caused by four typical types of point dislocation, derived by Sun & Okubo (1993 ), to deduce the fundamental formulas for spheroidal fields relevant to viscoelastic earth models. In computations, we employ a strike-slip dislocation on a vertical plane buried at the bottom of the lithosphere to estimate the maximal viscous relaxation responses to this kind of source that possibly exist on the surface of the earth. We take the seismic moment as 10²²  N  m, which is characteristic of an average large earthquake. The numerical results demonstrate that, if we take the viscosity as 10¹⁹  Pa  s in the asthenosphere, and 10²¹  Pa  s in the other mantle layers, the rates of surface vertical displacements and gravity changes within about 2.5° for the 10 postseismic years are respectively 1.5–8.1  cm  yr⁻¹ and 4.0–14.9  μgal  yr⁻¹ : the viscous relaxation for this mantle viscosity profile proceeds much faster than for a constant mantle viscosity of 10²¹  Pa  s.