Forward and inverse modelling of post-seismic deformation

We consider a new approach to both the forward and inverse problems in post-seismic deformation. We present a method for forward modelling post-seismic deformation in a self-gravitating, heterogeneous and compressible earth with a variety of linear and nonlinear rheologies. We further demonstrate how the adjoint method can be applied to the inverse problem both to invert for rheological structure and to calculate the sensitivity of a given surface measurement to changes in rheology or time-dependence of the source. Both the forward and inverse aspects are illustrated with several numerical examples implemented in a spherically symmetric earth model.


I N T RO D U C T I O N
It is well known that after an earthquake the Earth continues to deform over a period of months to years.There are three main processes that contribute to this time dependent deformation-afterslip, poroelastic rebound and viscoelastic rebound.Afterslip is continued aseismic movement on the original fault plane or neighbouring fault planes.Poroelastic rebound arises due to flow of fluid through rocks, which is driven by pressure changes associated with the coseismic motion.Lastly, there is further flow at depth due to viscoelastic relaxation in the Earth's crust and mantle.If we can isolate the deformation due to this final aspect, we can use measurements of the deformation to place constraints on the viscoelastic structure of the Earth.
However, the form of the Earth's rheology is not well known.The simplest viscoelastic rheology is that of a Maxwell solid, which has historically been considered sufficient to explain the form of post-seismic deformation and other viscoelastic processes such as post-glacial deformation.However, more recent studies have argued that some measurements of post-seismic deformation show evidence of a more complex rheology, such as a transient rheology with multiple decay times (e.g.Pollitz 2003) or a nonlinear rheology (e.g.Pollitz et al. 2001;Freed & Bürgmann 2004).The different orders of magnitude of the viscosities obtained from studies of post-seismic and post-glacial displacements gives further suggestion of this.
The study of post-seismic deformation involves two aspects.The first of these is forward modelling, where, given an earth model and earthquake source, the surface deformations can be calculated.The second, as mentioned above, uses measurements of this time-dependent deformation to constrain the structure of the Earth-the inverse problem.There are a number of challenges associated with both of these aspects.
First, it is common when forward modelling to solve the equations of motion in the Laplace transform domain.However, this method is not well suited to lateral variations in structure nor to radial structures that are not piecewise continuous (e.g.Fang & Hager 1995;Han & Wahr 1995;Boschi et al. 1999).As the Earth's rheological structure will be pressure and temperature dependent, quantities that vary continuously, it would be preferable to be able to permit continuous structure variation in the model.There is also recent evidence for the effect of lateral variations in structure on post-seismic deformation (e.g.Pollitz 2015).A further difficulty relates to the inclusion of self-gravitation.Some studies ignore its effect completely or include only an approximation to it (e.g.Pollitz et al. 2006;Wang et al. 2006).Another approximation that is sometimes made is that the Earth is incompressible (e.g.Boschi et al. 2000).
One of the main difficulties associated with the inverse problem is that of non-uniqueness.It is clear that due to the multiple processes contributing to the deformation and the wide variety of rheological structures possible, there will be many scenarios that yield essentially the same surface deformation.As an example of this, Freed et al. (2006) list the variety of mechanisms that have been suggested by different 846 O. Crawford et al. studies to have caused the time-dependent surface deformation after the 1992 Landers earthquake.This non-uniqueness emphasizes the need for evaluating which features of any model obtained are robust.
When performing an inversion from post-seismic data, most studies invert for a very small number of parameters and perform a gridsearch to find the best-fitting values of these parameters (e.g.Diao et al. 2014;Huang et al. 2016;Han et al. 2016).This limits the range of models that can be explored, and so raises questions of the physical significance of the features, such as discontinuities, of these models.A wider variety of models could be included by increasing the number of parameters, but this would also increase computation time.
Another approach is to use gradient-based optimization methods where the gradient of the misfit with respect to model parameters is calculated and used to lower the misfit iteratively (e.g.Nocedal & Wright 1999).A common way to calculate such gradients is to use finite difference methods, although this can also become computationally unfeasible for large numbers of model parameters.The adjoint method (e.g.Lions 1970) provides a more efficient way of calculating the gradient, and has been used for a number of geophysical applications (e.g.Bunge et al. 2003;Tromp et al. 2005;Li et al. 2011;Al-Attar & Tromp 2014).The exact gradient can be calculated with just one solution of the forward problem and one of the corresponding adjoint problem.
In this work, we present a method for modelling post-seismic deformation in a fully heterogeneous, compressible, self-gravitating Earth, which allows for both linear rheologies with multiple decay times and nonlinear rheologies.We also show how the adjoint method can be used to perform an inversion for earth models with continuous variation in structure and how the resolving power of the data can be quantified.Many aspects of this forward and inverse modelling are also applicable to the post-glacial problem, although this would further require inclusion of rotation and coupling of the deformation to sea level (e.g.Farrell & Clark 1976;Mitrovica & Milne 2003).
In Section 2, we present the theory relating to the forward viscoelastic problem, including both linear and nonlinear rheologies.In Section 3, we turn to the inverse problem and discuss the adjoint method and its application to calculating sensitivity kernels and performing viscosity inversions.We illustrate both the forward and inverse problems with several numerical examples using spherically symmetric earth models, but emphasize again that the theory allows for lateral variations in model structure.

T H E F O RWA R D V I S C O E L A S T I C P RO B L E M
We wish to formulate the equations of motion that describe post-seismic deformation in a self-gravitating, viscoelastic earth model.Our approach will be similar to that used in Al-Attar & Tromp (2014) to describe post-glacial deformation, but with the addition of a body force through the stress glut.Furthermore, while it is often considered sufficient to model the Earth's response as a Maxwell solid in the post-glacial case, we extend the equations to include other linear rheologies that are of interest in the post-seismic case.In Appendix A, we further describe how nonlinear rheologies could be included in a similar manner.

Equations of motion
Following Al-Attar & Tromp (2014) and making use of Dahlen & Tromp (1998), Dahlen (1974) and Tromp & Mitrovica (1999), we can formulate the equations of motion.We have an earth model that occupies a volume M ⊆ R 3 with external boundary ∂M.The model has a number of fluid and solid regions, the former denoted by M F and the latter by M S , that are separated by smooth, non-intersecting, closed surfaces called internal boundaries.The union of all boundaries, both internal and external, is denoted by , with the internal boundaries consisting of four subsets-SS , SF , FS and FF .Here, the first subscript denotes whether the region on the inside of the boundary is solid (S) or fluid (F), and the second subscript refers to the region on the outside of the surface.Neglecting inertial terms, the quasi-static momentum equation for a non-rotating, self-gravitating, hydrostatically pre-stressed, compressible and laterally heterogeneous earth model is where u is the displacement, ρ is the density, T is the Lagrangian-Cauchy stress tensor, is the reference gravitational potential, φ is the gravitational potential perturbation and is the stress glut (Backus & Mulcahy 1976).satisfies Poisson's equation, and, following Dahlen (1974), φ satisfies the modified Poisson's equation, where g is the gravitational acceleration and ∂ n is the directional derivative along the outward normal to the level surfaces of density.The stress glut is a body force that arises due to the localized failure of the assumed constitutive relation due to a seismic source (e.g.Dahlen & Tromp 1998).It is defined as Post-seismic deformation

847
where T model is the stress defined using the assumed constitutive relation, and T true is the physical or true stress.For many applications, it is sufficient to consider the earthquake to be a point source, for which the stress glut can be written as where M is the moment tensor, x 0 the hypocentral location and H(t) is the Heaviside step function.Finite sources can be built up by the superposition of multiple point sources.
We further need the boundary conditions for the system, which, again following Al-Attar & Tromp (2014) but generalizing to include the stress glut, are where + and − indicate whether a term is evaluated on the upper or lower side of a discontinuity respectively.We further require the gravity perturbation to tend to 0 as x tends to infinity, which we can write as lim x →∞ φ = 0. (2.16)

Linear viscoelasticity
In order to formulate the equations for the problem, we require a constitutive equation to specify the relationship between stress (T) and strain (e).In this section, we will consider only linear viscoelasticity, but demonstrate how the theory can be extended to include nonlinear viscoelasticity in Appendix A. For a classical, linear elastic solid, the constitutive equation is given by Hooke's law, where C is the fourth-order elastic tensor and e is given by e = 1 2 ∇u + (∇u) T . (2.18) In the case of an isotropic solid, eq.(2.17) simplifies to where κ and μ are the bulk and shear moduli respectively and d is the deviatoric strain, (2.20) Extending eq.(2.19) to an isotropic linear viscoelastic material using the Boltzmann superposition principle, we have where μ(t) is now the time-dependent relaxation function and we have neglected bulk viscoelasticity, which is thought to be a valid simplification (e.g.Wu & Peltier 1982).We could, however, include bulk viscoelasticity and, more generally, anisotropy through a simple extension of the internal variables method described below (Simo & Hughes 1998).Integrating the second term in eq.(2.21) by parts, we find (2.22)where μ 0 = μ(0) is, as we will see below, a measure of the elastic response of the material.For a wide class of linear viscoelastic materials, we can write μ(t) in the general Maxwell form (2.23) where each term in the sum represents a Maxwell element (Simo & Hughes 1998).We note that a relaxation function of this form allows only for linear rheologies with discrete decay times, and so we could not, for example, model an absorption band solid.However, the theory does include most linear rheologies of interest, and we will next demonstrate its equivalence to the formulation in the Laplace domain more commonly used in post-seismic literature (e.g.Piersanti et al. 1995;Pollitz 1997).
We recall that the Laplace transform, f (s), of a function f(t) is defined as where s is the transform parameter.It can be seen that the Laplace transform of eq.(2.21) is and the Laplace transform of eq.(2.23) is Using eq.(2.26), eq.(2.25) becomes (2.27) We can illustrate the equivalence with several simple examples of linear rheologies commonly used in post-seismic deformation.We shall first consider the case of a Maxwell solid.The relaxation function of such a material is given by (2.28) and so its Laplace transform is (2.29) Substituting eq.(2.29) into eq.(2.25) gives (2.30) Comparing eq.(2.30) to, for example, eqs (6) and (7) in Peltier (1974) shows that eq.(2.28) is indeed the relaxation function of a Maxwell solid.
We can further consider a Burgers body rheology, which is also commonly used in post-seismic literature (e.g.Pollitz 2003;Pollitz et al. 2006).This is the simplest transient rheology-it has two different time scales of decay.The relaxation function of a Burgers solid is with Laplace transform (2.38) We therefore see that a material with shear relaxation function given by eq.(2.23) behaves like a linear elastic solid for deformations that are fast relative to its decay times and like a Newtonian fluid for deformations that are slow relative to its decay times.Its exact behaviour on time scales between these two limits depends on the number of terms in eq.(2.23) and the relative magnitudes of the μ i s and τ i s.
We will now discuss how we can use the form of the relaxation function given in eq. ( 2.23) to model linear viscoelastic materials in the time domain.This is useful for numerical calculations and circumvents the issues with the Laplace transform method of modelling continuous variation in model parameters (e.g.Fang & Hager 1995;Han & Wahr 1995;Boschi et al. 1999).We define internal variables, m i , (e.g.Simo & Hughes 1998) such that (2.39) These are seen to satisfy and using this result, we can write eq.(2.22) as (2.41) As d is a trace-free tensor, the same is true of m i , and the deviatoric stress, T, is given by (2.42)

Rate formulation
In order to incorporate time dependence and a viscoelastic rheology into eq.(2.1), the momentum equation, we follow the method of Al-Attar & Tromp (2014).Differentiating eq.(2.1) with respect to time and using eqs (2.40) and (2.41), we obtain a rate formulation of the problem, We can similarly differentiate eqs (2.6)-(2.15) to obtain the rate formulated boundary conditions.These equations, along with eq.(2.40) and the initial conditions constitute the strong form of the problem.We now have a system of equations that implicitly give the time derivatives u, φ and ṁi at time t as a function of the current state u, φ, m i and .We can therefore calculate the time derivatives at a given time, and use them along with some time-stepping scheme to calculate the state of the system at a future time.This approach is equivalent to other time-domain methods in the geophysical literature (e.g.Hanyk et al. 1995;Zhong et al. 2003;Latychev et al. 2005), although it differs due to the introduction of internal variables.Furthermore, we have avoided explicitly introducing a numerical scheme for time-stepping, which allows for greater flexibility.

Weak form
In order to implement eq.(2.43) and the associated boundary conditions numerically using a finite-element-type method, we derive the so-called weak form of the problem.We can introduce test functions u , φ and m i , where we require that (2.46) and φ → 0 as x → ∞.Multiplying eqs (2.43) and (2.40) by these test functions and integrating by parts with respect to V, we obtain the weak form, (2.47) which for a symmetric stress glut can be written as ˙ : e dV = 0, (2.48) where e is the strain tensor given by e = 1 2 [∇u + (∇u ) T ]. (2.49) Henceforth, we will assume that the stress glut is symmetric and so we can use the weak form in eq.(2.48).As derived in appendix B of Al-Attar & Tromp (2014), A is given by We note that A is symmetric under interchange of (u, φ) and (u , φ ).It may be shown that all u, φ and m i that are solutions of this weak form are also solutions of the strong form, and vice versa.

Forward modelling in spherically symmetric earth models
To illustrate the forward method, we will consider our equations in spherically symmetric earth models.This allows the equations to be reduced to a system of 1-D equations in the radial coordinate.In order to do this, we expand the fields in generalized spherical harmonics (Phinney & Burridge 1973) which results in total decoupling of the radial expansion functions for each spherical harmonic of degree l and order m.We therefore need only construct a radial mesh, and then solve the radial equations using the Spectral Element Method (Patera 1984).
Further details of the generalized spherical harmonic expansions, particularly of the form of the source term, are given in Appendix B, but we state the basic equations here.In order to express eq.(2.48) using generalized spherical harmonics, we first recall the canonical basis vectors (Phinney & Burridge 1973;Dahlen & Tromp 1998), defined here relative to the basis vectors in spherical polar coordinates.The so-called contravariant components of the displacement, u, in this basis are (2.56) In order to expand these components using generalized spherical harmonics, we write where Y N lm are the generalized spherical harmonics defined in appendix C of Dahlen & Tromp (1998) and the summation is over integer values for 0 ≤ l ≤ ∞ and −l ≤ m ≤ l.We will find it useful to define the coefficients U lm , V lm and W lm such that where k = l(l + 1). (2.60) We can also expand the spatial derivatives of u and all other fields in terms of generalized spherical harmonics, and these expansions are given in Appendix B. In practice, we cannot sum to infinite l, and so truncate the series at l = l max , which is chosen to achieve convergence to a desired level of accuracy.
Our equations are implemented in an earth model with a solid elastic inner core, a compressible inviscid outer core, a viscoelastic mantle with a linear rheology of the form given in eq. ( 2.23) and an elastic lid.Our implementation allows for discontinuities between model layers, but also for continuous variation in model parameters within these layers.This is in contrast to the more common implementation in the Laplace domain, for which problems arise unless the model is piecewise continuous (e.g.Fang & Hager 1995;Han & Wahr 1995;Boschi et al. 1999).
At each time point, we use the current fields to calculate their time derivatives by solving eq. ( 2.48) and time step using an explicit second order Runge-Kutta integration scheme (e.g.Press et al. 1986).The time step is taken to be half the shortest decay time in the model, which is sufficient for numerical stability and to obtain the desired level of accuracy.
Unless otherwise stated, in the examples below our source is a thrust fault earthquake at 15 km depth, striking North and with a dip of 30 • .Our source is an approximation to a point source-its spatial distribution is a delta function in θ and φ, but a Gaussian with a standard deviation of 2.0 km in the radial direction.This is chosen for convenience, and we note that it would be possible to choose any other radial function.We use the elastic structure of the Preliminary Reference Earth Model (PREM; Dziewonski & Anderson 1981), with a 19 km thick elastic lid and a viscosity of 6 × 10 18 Pa s from below this layer to a depth of 400 km, 5 × 10 20 Pa s from a depth of 400 km to 670 km and 5 × 10 21 Pa s from 670 km to the core-mantle boundary.The shallowest viscosity is taken from Huang et al. (2016), and the deeper values are common in post-glacial simulations; however, the post-seismic displacements are unlikely to be sensitive to these deeper viscosities.In general, we use a Maxwell solid rheology, though we include an example using a Burgers body rheology.We note again that the theory described allows for other linear and nonlinear rheologies, although the latter necessarily introduces lateral heterogeneity into the problem due to the dependence of the effective viscosity on stress, a laterally varying field.Therefore, nonlinear rheologies are not currently implemented.More details on the theory and method in the case of a nonlinear rheology are given in Appendix A.

Comparison with analytic solutions
In order to benchmark the code, we compared our numerical solutions to those obtained analytically in a few specific cases.For a non-selfgravitating, spherically symmetric earth model, it is possible to derive an analytic solution for some simple models.This can be done using propagator matrix methods (Gilbert & Backus 1966;Woodhouse & Deuss 2007;Al-Attar & Woodhouse 2008).The source has the same form as described above, except that is a delta function in the radial direction.We obtained analytic solutions for the initial values of the U lm and V lm coefficients and time-dependent values of the W lm coefficients for a homogeneous elastic structure-the parameters from the base of the mantle to the surface are equal to the values at a depth of 19 km in PREM.The viscosity is 6 × 10 18 Pa s throughout the entirety of the mantle.These solutions were compared to the numerical solutions for the earthquake described above for many values of l and m, and were found to agree to within numerical precision, except within the vicinity of the source where the fact the numerical source only approximates a delta function increases the error slightly.Examples are shown in Figs 1 and 2.

Comparison with Laplace domain code
Our code was further compared to an existing Laplace domain code, based on Wu & Peltier (1982), which models post-glacial deformation.For this case, we must replace the stress glut term in eq.(2.48) with a term that represents a surface load; from Al-Attar & Tromp (2014), this term is for some time-dependent surface load σ (x, t).The surface load used was The analytic solutions are derived using propagator matrix methods.The source is as described in Section 2.5 and the earth structure as in Section 2.5.1.The solutions obtained by the different methods agree very well.for t ≥ 0, where l max = 256 and L = 100.Its value as a function of latitude is shown in Fig. 3.For this test, we used an earth structure with a 120 km thick elastic lid, a viscosity of 5 × 10 20 Pa s in the upper mantle and a viscosity of 5 × 10 21 Pa s in the lower mantle.Fig. 4 shows the radial displacement and the gravity perturbation calculated by the two codes initially and 10 000 yr after the load was added.We can see that the two codes agree very well, with any differences of order 1 per cent.

Forward modelling examples
In Fig. 5, we show the surface displacement at two different times for the earthquake and earth structure described above calculated using our model.Fig. 6 shows the depth dependent displacements for a slice through the Earth for the same earthquake.Fig. 7 shows the horizontal displacement at the epicentre for both a Maxwell solid and Burgers solid mantle rheology.The parameters of the Burgers body relaxation function (eq.2.31) are chosen such that the two rheologies have the same shear modulus (given by eq.2.35) and long-term viscosity (eq.2.38) and so, as the figure shows, the displacements will have the same initial value and long-term gradient.Using parameters similar to those found by Pollitz (2003), we have chosen μ 1 = μ 2 and τ 1 = 30τ 2 .

T H E I N V E R S E P RO B L E M
In the previous section, we discussed how we could model post-seismic deformation in an earth model with a given structure.We will now turn our attention to the inverse problem-using measurements of post-seismic deformation to make inferences about the structure of the Earth.Given a set of surface displacement data, we seek the earth structure that minimizes the misfit between this data and the calculated displacements.In practice, this can be done using iterative methods for which we require the gradient of the misfit with respect to a model parameter, and we will calculate these gradients using the adjoint method.We will also explore the sensitivity of particular measurements to changes in model parameters.This provides useful insight into the physics of the problem and into the resolving power of the data.

The adjoint method
The adjoint method enables derivatives of measurements with respect to model parameters to be calculated quickly, with only one solution of the forward problem and one solution of the corresponding adjoint problem required (e.g.Lions 1970;Tromp et al. 2005).In this section, we will briefly discuss the theory behind the adjoint method and illustrate it schematically.
We define an observable, J, of which we wish to calculate the derivative with respect to a particular model parameter.We could, for example, choose J to be a misfit function between observed and synthetic data, and could use the derivative as part of an iterative scheme to improve our model (e.g.Tape et al. 2007).We could, alternatively, choose J to be a particular surface measurement-in this case, the derivative would represent the sensitivity of the measurement to the chosen model parameter.J will typically be a function of the forward variables, but may also depend explicitly on the model parameters.We therefore introduce an objective functional, J(u, p), where u represents the variables of the forward equations (in our case, the displacements, u, the gravity perturbation, φ, and the internal variables, m i ) and p represents the model parameters.As the forward variables will be functions of the model parameters, we can define the reduced objective functional Ĵ ( p) = J (u( p), p). (3.1) We wish to calculate the derivative of this quantity with respect to p.We will denote this derivative by D Ĵ ( p), which is defined such that where •, • is an appropriate inner product.The linearized change in Ĵ ( p) due to a change δp in p is It is not efficient to calculate the derivative of Ĵ ( p) directly.We instead calculate the derivative of the objective functional, J(u, p), with respect to p, subject to the constraint that the forward variables, u are the solutions of the forward problem described by p.For the class of The radial displacements and geoid anomalies as a function of latitude computed using the approach of this study and a code based on Laplace domain methods (see Section 2.5.2).The green and red points are associated with the left axes, whereas the blue and pink are associated with the right.The force is an axially symmetric surface load shown in Fig. 3 and given by eq.(2.62).The comparison is shown at (a) 0 yr and (b) 10 000 yr after the load is first applied.
problems we wish to consider, the constraints can be written as a(u, p) = 0. Therefore, in order to calculate the adjoint equations, we introduce a Lagrangian functional defined as  provided where D p L is the partial derivative of L with respect to p, and similarly for u and u † .We can therefore say that the variation of Ĵ with respect to p is equal to the variation of the Lagrangian with respect to p, provided the Lagrangian is stationary with respect to perturbations of both the forward and adjoint variables.Evaluation of eq.(3.7) simply returns the forward equations for our problem through the weak form given in eqs (2.48), and evaluation of eq.(3.6) gives the corresponding adjoint equations.Therefore, in order to calculate the derivative of Ĵ with respect to a given model parameter, we must find the solutions of eqs (3.6) and (3.7), and substitute these into the derivative of the Lagrangian with respect to that parameter.

The adjoint equations for the post-seismic problem
We will now derive the adjoint equations for our post-seismic problem.The solutions to these and the forward equations can then be used to calculate the derivatives needed for the inverse problem.The Lagrangian for the post-seismic problem is where J is the chosen observable and the remaining terms are the time-integrated weak form of the post-seismic problem, as given in eq.(2.48) with A defined as in eq.(2.50).We note that we are ignoring any possible explicit dependence of J on a model parameter, which is sufficient for our purposes.We find that eq.(3.7) simply returns eq.(2.48) (with our test functions u , φ and m i replaced by δu , δφ and δm i ).Evaluation of eq.(3.6) gives where we have written the first-order perturbation of J with respect to u and φ as Here, ḣ is the Fréchet kernel of J with respect to u, and ḣ is the Fréchet kernel of J with respect to φ.It will now be useful to define the adjoint variables and adjoint tractions Integrating eq.(3.9) by parts with respect to time, and using the fact that it must hold for all times, we obtain where the forward variables are evaluated at time t whereas the adjoint variables are evaluated at time T − t.The adjoint variables satisfy the initial conditions which arise from the requirement that the boundary terms created by integrating by parts to obtain eq.(3.16) vanish.This constitutes the weak form of the adjoint problem, where δu, δφ and δm i act as the time-independent test functions.We note that this equation is identical in form to eq. (2.48) with the exception that the force term is now dependent on the form of J through the adjoint tractions h † and h † .This means we can use the same numerical code to solve both the forward and adjoint equations.

Sensitivity kernels
We now wish to use the solutions of the forward and adjoint equations to calculate how Ĵ changes in response to a change in a given model parameter.We will first consider the sensitivity of a particular surface measurement to changes in the model parameter and in this case, J will be the chosen measurement.We could, for example, choose J to be the displacement in a given direction at a time t s and at a point (θ s , φ s ) on the surface and so it will be given by where d is the unit vector in the direction of the displacement measurement, R is the radius of the Earth and T ≥ t s .The adjoint tractions will therefore be Calculation of sensitivity kernels gives insight into what can be learnt about the Earth from individual surface measurements.We will consider such applications further in Section 3.4.3.In the following subsections, we will derive the form of the sensitivity kernels for viscosity and the time-dependence of the source, but note that it would be simple to calculate the kernels for other parameters, such as the elastic structure.

Application of the adjoint method to viscosity sensitivity kernels
When calculating the sensitivity of surface measurements to viscosity, it will be useful to consider relative changes in viscosity.We therefore write the linearized change in Ĵ as (3.20)where η i = μ i τ i is the ith component of the viscosity, δη i is the change in this component and K η i is the derivative of Ĵ with respect to η i , which is commonly known as the sensitivity kernel.We will consider changes in η i due to changes in τ i only, and so assume the elastic structure does not change.Perturbing eq.(3.8) with respect to η i and using eq.(3.5), we find and so We note that eq.(3.22) gives the sensitivity to the ith component of the viscosity, η i = μ i τ i ; however, we can calculate the total viscosity sensitivity kernel by summing the individual contributions.We can therefore use the solutions of the forward and adjoint equations and eq.(3.22) to calculate the sensitivity kernel at all spatial positions in the Earth.However, for some applications, it will also be useful to calculate the sensitivity to viscosity as a function of depth only.We shall denote this radial sensitivity kernel by K η i and it is given by where S 2 is the unit two-sphere.

Application of the adjoint method to stress glut sensitivity kernels
As the stress glut is a time-dependent quantity, we will define the stress glut sensitivity kernel such that We note that, given the form of eq.(3.8), it is preferable to calculate the kernel with respect to ˙ rather than .Perturbing eq.(3.8) with respect to ˙ , we find and so where e † is the strain of the adjoint field (cf.Tromp et al. 2005).

Laterally varying viscosity kernels in spherically symmetric earth models
To illustrate the form of the kernels, we calculate the laterally varying kernel in a spherically symmetric earth model.As described in Section 2.5, we expand the fields in generalized spherical harmonics.Expansions of the stress fields and internal variables are given in Appendix B. In order to implement the adjoint tractions in eq. ( 3.19), we also require an expansion of the delta function in generalized spherical harmonics; discussion of this expansion for delta functions in the radial direction is given in appendix E of Al-Attar & Tromp (2014) and in the tangential direction in Appendix C of this work.
To calculate the kernel, we must (i) solve the forward equations given by (2.48) to calculate the spherical harmonic coefficients of the forward displacement field, U lm , V lm and W lm , and of the internal variables, m i ; (ii) solve the adjoint equations given by (3.16), with adjoint tractions from eq. (3.19), to calculate the spherical harmonic coefficients of the adjoint displacement field, U † lm , V † lm and W † lm , and of the internal variables, m † i ; (iii) calculate the spherical harmonic coefficients of the strain fields, d and d † ; (iv) calculate the spatial dependence of d, d † , m i and m † i by expanding the spherical harmonic representation; (v) calculate the sensitivity kernel using eq.(3.22).
We calculate both the forward and adjoint fields in coordinate systems defined such that their respective sources are at the North Pole.This greatly increases computational efficiency as the spherical harmonic coefficients are non-zero only for particular values of m. (For the forward fields, the coefficients are non-zero for −2 ≤ m ≤ 2 whereas for the adjoint fields, they are non-zero for −1 ≤ m ≤ 1 as the adjoint source is a vector quantity for point measurements.)The spatial dependence of the two fields is calculated using the angles between the actual locations of the sources and the North Pole.
For the source and Earth structure described in Section 2.5, we plot the kernel for two different displacements 20 yr after the earthquakethe vertical displacement at the epicentre (Fig. 8) and the displacement to the east at a point 0.1 • west of the epicentre (Fig. 9).

Radial viscosity kernels in spherically symmetric earth models
In order to calculate the radial kernel in a spherically symmetric earth model, we can write eq.(3.23) using generalized spherical harmonics.This expansion is given in Appendix B2.We again require the spherical harmonic expansion of the delta function, as described in Section 3.3.3.We see that in order to calculate the radial kernel, we must (i) solve the forward equations given by (2.48) to calculate the spherical harmonic coefficients of the forward displacement field, U lm , V lm and W lm , and of the internal variables, m i ; (ii) solve the adjoint equations given by (3.16), with adjoint tractions from eq. (3.19), to calculate the spherical harmonic coefficients of the adjoint displacement field, U † lm , V † lm and W † lm , and of the internal variables, m † i ; (iii) calculate the spherical harmonic coefficients of the strain fields, d and d † ; (iv) calculate the sensitivity kernel using eq.(B35).
In contrast to the calculation of laterally varying kernels, we require the spherical harmonic coefficients to be defined relative to the same coordinate system.We choose the coordinate system to be that which places the adjoint source at the North Pole, as this minimizes the number of m values for which the spherical harmonic coefficients are non-zero.This is similar to the method used by Nissen-Meyer et al. (2007) for calculating seismic kernels in an axisymmetric setting.
In Figs 10 and 11, we show several examples of radial sensitivity kernels for the earthquake source and earth structure described in Section 2.5.Fig. 10 shows the kernels for vertical and horizontal displacements at multiple times at the epicentre.Fig. 11 shows the kernels for the vertical and horizontal displacements at varying locations 20 yr after the earthquake.
In order to emphasize that the kernels give the linearized sensitivity of a particular measurement to change in viscosity, we compare the changes in displacement predicted by the kernel to the actual changes in displacement due to viscosity perturbations of different magnitudes.This comparison is shown in Fig. 12.We can see that the actual change in displacement is close to the linear kernel prediction for perturbations of up to approximately an order of magnitude, and that above this, the kernel is still of the correct sign, which is important for carrying out a gradient-based inversion.The linearized nature of the kernel is further shown by the dependence of the kernel on the viscosity model used to calculate it, as illustrated in Fig. 13.

Application of the adjoint method
We will now define J to be the misfit between given displacement data and displacements calculated from our forward model.We wish to minimize this misfit in order to find the earth model that provides the best fit to our displacement data.Suppose we have displacement data at N locations on the surface, with these locations given by x i for 1 ≤ i ≤ N. At each location, there are N i measurements that each give the displacement at a particular time in a particular direction.Thus u obs i j is the datum for the displacement at location x i and time t ij and in direction dij .We can then define the misfit between the displacement data and calculated displacements to be (3.28)where for simplicity, we ignore potential data errors.We note that we could include regularization in the inversion by adding a term to J with an explicit dependence on model parameter, but its contribution to the gradient can be calculated directly.Perturbing eq.(3.28) with respect to u, we find that the adjoint tractions are We see that the adjoint traction, ḣ † (t), in this case is simply the sum of the viscosity kernel tractions for the individual measurements weighted by the difference between the observed and calculated displacements.Therefore, the linearized change in the misfit, J, with viscosity is (3.30)where K i j η k is the kernel associated with the kth viscosity component for a measurement made at location x i and at time t ij of the displacement in the direction dij .The gradient of the misfit, g k , with respect to changes in the kth viscosity component, η k , is defined such that (3.31) Figure 11.The radial viscosity kernels, given by eq.(3.23), for the displacement at varying locations, 20 yr after the earthquake.The earthquake source and earth structure are as described in Section 2.5.The labels give the locations of the displacements in degrees relative to the epicentre of the earthquake.The kernels shown are for (a) the vertical displacement and (b) the displacement to the east.The spatial dependence of the surface displacement at this time for this earthquake is as shown in Fig. 5. and so We can also write the gradient of the misfit with respect to changes in viscosity as a function of depth only.This will be a function of the radial kernels, K i j η k , and is given by (3.33)

Viscosity inversion in spherically symmetric earth models
Given surface displacement data, we wish to find a radial viscosity profile that fits the data to within error.We can again use generalized spherical harmonics to expand the fields as described in Section 2.5.In order to perform an inversion, we (i) calculate the surface displacements for some initial model at the locations where there are data measurements, and calculate the misfit using eq.(3.28); (ii) calculate the radial sensitivity kernels for each measurement, as described in Section 3.3.4;(iii) calculate the gradient of the misfit function by evaluating the weighted sum of radial kernels given in eq.(3.33); (iv) calculate a model update using, for example, the conjugate gradient method (e.g.Nocedal & Wright 1999;Tape et al. 2007); (v) repeat steps (i) to (iv) until the value of the misfit fails to decrease by some suitable amount.In order to minimize the calculations required, and hence improve the computational efficiency, the adjoint equations are solved only twice at each iteration-once with a vertical traction and once with a horizontal.The location and orientation of the earthquake are then defined relative to the adjoint traction for each kernel calculation.
For several models, we have calculated synthetic data for the earthquake source described in Section 2.5, and then inverted for viscosity structure using the method detailed above.The locations, times and directions of the displacement data used are shown in Table 1.Unless otherwise stated, the model chosen to begin the inversion was that described in Section 2.5.Several examples of the actual viscosity structures and results of the inversions are shown in Figs 14-17, along with the corresponding displacement misfits, given by eq.(3.28), and viscosity Table 1.The locations (relative to the epicentre), times and directions of displacements used in carrying out the viscosity inversions.The displacements to the north are not used for the locations at 0 • latitude as these displacements are zero due to the geometry of the source.The actual viscosities used to generate the synthetic data and final inverted viscosities are shown in (a), the former using dotted lines and the latter using solid lines of the same colour.The corresponding displacement misfits, given by eq.(3.28), at each iteration of the inversion are shown in (b).The earthquake source and starting structure for the inversion are as described in Section 2.5.misfits where appropriate.The viscosity misfit after the ith iteration of the inversion is defined to be 1 2 where r 1 is the radius of the base of the mantle, r 2 is the radius of the top of the mantle, η act is the actual viscosity and η i is the viscosity obtained from the inversion after the ith iteration.First, we observe that the final models obtained in the cases shown in Fig. 14 are similar, and that both have a small final displacement misfit, despite one model having a discontinuous change in viscosity and one having a continuous change which occurs over approximately 10 km.We therefore conclude that, in this case, our data are unable to resolve sharp changes in viscosity, which is perhaps unsurprising given the physics of viscoelastic relaxation.
We have also investigated the ability of the data to resolve a discontinuous change in viscosity at different depths.In Fig. 15, we show the results of an inversion using synthetic data generated from a viscosity structure which jumps from 1.2 × 10 19 Pa s to 6 × 10 18 Pa s at three different depths-30, 40 and 50 km.We see that there is a substantial difference in the inverted viscosities in the case of the 30 and 40 km discontinuities, and both have a low viscosity misfit.However, the inverted viscosity in the case of the 50 km discontinuity is very similar to Figure 15.Viscosity inversions with a viscosity discontinuity at different depths.The actual viscosities used to generate the synthetic data and the final inverted viscosities are shown in (a), the former using dotted lines and the latter using solid lines of the same colour.The corresponding viscosity misfits, given by eq.(3.34), at each iteration of the inversion are shown in (b).The earthquake source and starting structure for the inversion are as described in Section 2.5.that of the 40 km discontinuity, and has a significantly higher viscosity misfit, despite all three models having a very low displacement misfit (4.1 × 10 −6 , 1.2 × 10 −5 and 7.0 × 10 −5 m 2 for the 30, 40 and 50 km discontinuities respectively, the largest of which corresponds to an average error per measurement of about 2 mm).We therefore conclude that in this case, the data struggles to resolve a change in viscosity at a depth of 50 km.This makes sense given the form of the sensitivity kernels in Figs 10 and 11, as the kernels are very small in magnitude at depths of 50 km and below.
In order to illustrate the dependence of the inversion results on the initial viscosity structure chosen, we have performed an inversion for the same viscous structure using two different starting models, as shown in Fig. 16.It can be seen that the final viscosity structures in the two cases are very different, despite convergence having occurred as the misfit is roughly constant for the last few iterations.The misfit in one case is notably higher, which suggests the inversion has found a local misfit minimum.Despite this, the misfits in both cases are still low enough to be virtually indistinguishable, particularly given a reasonable level of uncertainty in any real surface displacements.However, we also note that in one case, the viscosity change is very steep just below the elastic layer.Such a structure may be considered to be unphysical, and so the inversion could be improved by adding a regularization term.
All the inversions discussed above have used synthetic data with no errors.In practice, the data will not be error free and so we have investigated the effect of errors on a particular inversion.Having generated synthetic data for the model shown by the green dotted line in Fig. 14(a), we added random errors from a normal distribution with a standard deviation of 1 cm.In Fig. 17, the results for three inversions with different random errors selected from this distribution are shown.We see that the inclusion of errors results in model structures with higher viscosity misfits, although the difference is not obviously significant.Furthermore, we see that the final viscosity with the lowest displacement misfit of the three with data errors does not have the lowest viscosity misfit.As mentioned previously, adding a regularization term could also improve these inversions.
In this section, we have investigated a variety of viscosity inversions using post-seismic data.However, these examples are intended to illustrate the potential of the adjoint method and not necessarily to demonstrate exactly how an inversion could be carried out.We note in particular the dependence of the inversion results on the initial structure chosen as shown in Fig. 16, and the possibility for adding regularization to the inversion.We further emphasize that care must be taken in the interpretation of the inverted viscosity profile, particularly in the case of data errors, and we will consider a method for quantifying the resolution of the data in the following section.
Figure 16.Viscosity inversions using the same synthetic data but two different starting models.The actual viscosity used to generate the synthetic data, along with the starting viscosities and final viscosities in the two cases are in (a), the former using dotted lines and the latter using solid lines of the same colour.The corresponding displacement misfits, given by eq.(3.28), in the two cases are shown in (b).The earthquake source is as described in Section 2.5.

Resolution
We wish to quantify which parts of the model a given measurement, or set of measurements, can constrain, and how accurately it can do so given a level of uncertainty in our surface measurements.As we have seen, the linearized change in a surface measurement can be related to the change in viscosity, η, through (3.35)where u i is the ith surface measurement with sensitivity kernel K i , and I S is the set of radii within solid regions.Here we are assuming a Maxwell solid rheology, but the method could be extended by a generalization of what follows.We can write this relationship for all measurements as where δu is a vector of all the surface measurements and A is a linear operator from the model to data space constructed from the sensitivity kernels.Given N measurements, we can calculate a maximum of N model parameters.However, the kernels are not necessarily linearly independent (e.g.Gilbert 1971).A way to assess what can be resolved is to look at the singular value decomposition of A (e.g.Parker 1994).
A similar approach has been considered by Pollitz & Thatcher (2010) and Hines & Hetland (2016), and has also been used for post-glacial applications (e.g.Mitrovica & Peltier 1991).
It can be shown (e.g.Golub & van Loan 1983) that for any real matrix, there exist orthogonal matrices, L and R, such that (3.37) where is a matrix with diagonal components only.If we let l i denote the ith column of L and r i denote the ith column of R, it is easy to see that (3.38) (3.39) Figure 17.The effect on viscosity inversions of adding errors to the synthetic data.The actual viscosity and result of the inversion using the synthetic data is shown in (a) along with three inversions where the synthetic data have had random errors added from a normal distribution with a standard deviation of 1 cm.The corresponding (b) displacement misfits, given by eq.(3.28), and (c) viscosity misfits, given by eq.(3.34), are also shown.The earthquake source and starting structure for the inversion are as described in Section 2.5.

867
where λ i is the ith diagonal component of the matrix .We say that λ i , l i and r i are the ith singular value, left singular vector and right singular vector of A respectively.As L and R are orthogonal matrices, it is clear that l i , l j = δ i j , (3.40) (3.41)where •, • are the appropriate inner products.From eqs (3.38) and (3.39), we can further see that (3.43) In the case described by eq.(3.36), each left singular vector will be a vector of length N, the number of surface measurements, whereas each right singular vector will be a continuous function of radius.We will therefore write r i for the right singular vectors to indicate that they are continuous scalar functions.As AA † is simply a square matrix of dimension N, we can compute its eigenvalues, λ 2 i , and eigenvectors, l i , and then calculate the right singular vectors, r i , using eq.(3.38).The l i will provide an orthonormal basis for the data space, and the r i provide the corresponding orthonormal basis for a finite-dimensional subspace of the model space.We see that the measurements δu are only affected by relative changes in model that can be spanned by the r i basis.We can now quantify how well we can constrain each orthogonal part of the model space given by r i using the combination of data points given in the corresponding l i .
Suppose we have some surface displacement data and used the method described in Section 3.4.2 to find a solution to the inverse problem described by this data.Given that there will be errors associated with our data, we can ask how well constrained our model is, that is, how much we can perturb the model while keeping the surface displacements within the error bounds of our data.Let us consider a perturbation of the form δη η = br i (3.44) for a particular i, where b is a coefficient determining the size of the perturbation and η is the viscosity structure obtained in the inversion.Using eq.(3.39), the corresponding change in surface displacements will be where A is a linear operator constructed from the sensitivity kernels calculated with respect to our model, η.If this change in displacements leaves them within the error bounds of our data, then our data is unable to distinguish between the original and perturbed models.Therefore, in order for a given perturbation to be resolvable, the change in displacements must be greater than the error in the displacements.We can denote the errors in the displacement data by the vector , which we assume has zero mean and known covariance matrix, C, which we recall is defined such that We can further use this approach to say how well our model is constrained subject to perturbations of the form br i .Say, for example, that we wish to know which perturbations we can constrain sufficiently well such that the viscosity at any point in the model cannot change by more than a fraction f.In this case, we require where R i is the maximum value of the ith right singular value, r i .Combining this with eq.(3.49), we see that and so for perturbations of the form of eq.(3.44),only those with singular values for which eq. (3.51) holds are sufficiently well constrained.
As an illustration, we will consider the resolution of the example inversion in dark blue in Fig. 17.We can calculate the sensitivity kernels for each of the surface measurements in Table 1 for the final model obtained in this inversion and so construct the matrix A and calculate 868 O. Crawford et al.
Figure 18.The singular values in order of decreasing magnitude of the matrix A whose rows are given by the sensitivity kernels calculated for the measurements in Table 1 with respect to the final viscosity structure shown in dark blue in Fig. 17.
the matrix AA † .From eq. (3.42), we can see that the eigenvalues of this matrix will be the singular values squared and the eigenvectors will be the left singular values.Using eq.(3.38), we can then calculate the right singular vectors.The λ i are shown in Fig. 18 in order of decreasing magnitude and the right singular vectors with the four largest singular values are shown in Fig. 19.We see that the best constrained right singular vector (i.e. the one with the largest singular value) has its maximum amplitude at the top of the viscoelastic region and the amplitude decays as depth increases.Furthermore, as the singular value decreases, the associated right singular vector has a greater number of oscillations.
Let us assume, for example, that the standard deviation of the error in each surface measurement is 1 cm, and that the errors are uncorrelated.The trace of the covariance matrix will then simply be 0.01 √ 30 ≈ 0.05 as there are 30 displacement measurements.From Fig. 19, we can see that R i ≈ 0.03m −1/2 .If we wish to find which perturbations we can constrain such that the viscosity at any point varies by up to 10 per cent, we find from eq. (3.51) that λ i > 0.015, (3.52) and so, by examining Fig. 18, we can say our model is resolved to within 10 per cent subject to perturbations of the form of one of the right singular vectors with the seven largest singular values.

C O N C L U S I O N S
We have presented a new method for modelling post-seismic deformation in a compressible, self-gravitating and laterally heterogeneous earth model that can have a variety of linear and nonlinear rheologies as well as continuous variation in structure.This method permits the inverse problem to be solved using gradient-based optimization, and application of the adjoint method allows gradients to be calculated from a single solution of the forward problem and one of the corresponding adjoint problem.We have illustrated both the forward and inverse problems with several examples calculated in a spherically symmetric earth model.There are several possible avenues of future work.In this study, we have considered the deformation of the whole Earth.However, post-seismic deformation is largely localized to the vicinity of the seismic source and so we could increase computational efficiency by modelling the deformation in a plane layer geometry.One possible method would be to use a cylindrical coordinate equivalent of the generalized spherical harmonics introduced by Burridge (1984).It is not, however, straightforward to include self-gravitation in a plane layer (e.g.Rundle 1980), but its effect is arguably not vital in some situations.
Furthermore, the theory presented in this study is applicable to fully heterogeneous models.Implementation in a 3-D finite-element code, for example, would allow for both forward and inverse calculations to be carried out in laterally varying models.However, it is not clear if post-seismic data could usefully constrain any 3-D viscoelastic structure in the Earth.
We derived the form of the viscosity kernels and stress glut kernels, but only gave examples of the former.Using kernels for the time-dependence of the source, we could perform a joint inversion for afterslip and viscosity structure, and investigate the resolution in the case that the particular cause of the time-dependent displacement is unknown.It would also be possible to calculate the kernels for other model components, such as the elastic parameters or thickness of the elastic lid.
Figure 19.The right singular vectors, r i , with the four highest singular values, of the matrix A whose rows are given by the sensitivity kernels calculated for the measurements in Table 1 with respect to the final viscosity structure shown in dark blue in Fig. 17.
Many of the methods presented in this paper are also applicable to post-glacial deformation, although full application would require the addition of rotation and sea level into the problem.Investigation of the differing depth dependence of the viscosity kernels in the post-seismic and post-glacial cases, and inclusion of transient and nonlinear rheologies could lead to a possible explanation for the differing orders of magnitude of viscosity required to explain the surface displacements in the two cases.

B2 Radial Sensitivity Kernel
We wish to expand eq.(3.23), which is in generalized spherical harmonics.The adjoint coefficients are defined in terms of the complex conjugates of the generalized spherical harmonics.Using the spherical harmonic expansions of d and m i given by eqs (B11)-( B14) and ( B15)-(B18) respectively, we find that In practice, we cannot sum to infinite l and so would have to truncate the sum in eq.(C8).However, this will lead to ringing.Therefore, we instead expand the smoothed delta function as where L determines its width.

Figure 1 .
Figure 1.Comparison of the initial values of the displacement coefficients U lm , V lm and W lm for l = 5 and m = 2 for both analytic and numerical solutions.The analytic solutions are derived using propagator matrix methods.The source is as described in Section 2.5 and the earth structure as in Section 2.5.1.The solutions obtained by the different methods agree very well.

Figure 2 .
Figure2.Comparison of the analytic and numerical solutions for the time dependence of the value at the surface of the displacement coefficients W lm for l = 5 and m = 2, and l = 20 and m = −2.These coefficients have zero real part due to the geometry of the source.The analytic solutions are derived using propagator matrix methods.The source is as described in Section 2.5 and the earth structure as in Section 2.5.1.The solutions obtained by the different methods agree very well.

Figure 3 .
Figure 3. Axially symmetric surface load given by eq.(2.62) used to compute the displacements in Fig. 4.

Figure 4 .
Figure 4.The radial displacements and geoid anomalies as a function of latitude computed using the approach of this study and a code based on Laplace domain methods (see Section 2.5.2).The green and red points are associated with the left axes, whereas the blue and pink are associated with the right.The force is an axially symmetric surface load shown in Fig.3and given by eq.(2.62).The comparison is shown at (a) 0 yr and (b) 10 000 yr after the load is first applied.

Figure 5 .
Figure 5.The surface displacements due to the dip slip earthquake described in Section 2.5 at 15 km depth.The initial displacement is shown in (a) and the change from this displacement 20 yr after the earthquake in (b).The colour gives the vertical displacement and the arrows show the direction and size of the horizontal displacement.

Figure 6 .
Figure6.The displacements at depth due to the dip slip earthquake described in Section 2.5 at 15 km depth.The initial displacement is shown in (a) and the change from this displacement 20 yr after the earthquake in (b).The colour gives the magnitude of the displacement and the arrows show the direction and relative magnitude of the displacement.There is no displacement perpendicular to the plane of the figures due to the geometry of the earthquake source.

Figure 7 .
Figure7.The horizontal displacement at the epicentre due to the earthquake described in Section 2.5, for a Maxwell solid rheology and a Burgers solid rheology.The parameters for the latter are described in Section 2.5.3.With the Burgers body, there is a rapid change in displacement initially, due to the short decay time.At long times, as expected, the Burgers body and Maxwell solid exhibit similar behaviour.

Figure 8 .
Figure 8.The laterally varying viscosity kernel, given by eq.(3.22), for the vertical displacement at the epicentre, 20 yr after the earthquake.The earthquake source and earth structure are as described in Section 2.5.The value of the kernel is shown for (a) varying latitude and (b) varying longitude.

Figure 9 .
Figure 9.The laterally varying viscosity kernel, given by eq.(3.22), for the displacement to the east, 0.1 • (≈11 km) west of the epicentre, 20 yr after the earthquake.The earthquake source and earth structure are as described in Section 2.5.The value of the kernel is shown for (a) varying latitude, (b) varying longitude and (c) variation in both latitude and longitude at a depth of 21 km.

Figure 10 .
Figure10.The radial viscosity kernels, given by eq.(3.23), for the displacement at different times after the earthquake at the epicentre of the earthquake.The earthquake source and earth structure are as described in Section 2.5.The kernels shown are for (a) the vertical displacement and (b) the displacement to the east.

Figure 12 .
Figure12.Comparison of the change in displacements predicted by the sensitivity kernels to the actual changes due to viscosity perturbations of different magnitudes between 19 and 30 km depth.The displacement in the r direction is at the epicentre, whereas the displacement in the φ direction is 0.1 • west of the epicentre.Both are the displacements 20 yr after the earthquake.The lines show the decrease in displacements predicted by the kernels whereas the crosses shown the actual decrease in displacements.

Figure 13 .
Figure 13.An example of the dependence of the sensitivity kernel on the background model used.The radial viscosity kernels, shown in (b), for the vertical displacement at the epicentre 20 yr after the earthquake were calculated with the two different background models in (a).

LatitudeFigure 14 .
Figure14.Viscosity inversions using synthetic data generated from a structure with a discontinuous change in viscosity and a continuous change occurring over approximately 10 km.The actual viscosities used to generate the synthetic data and final inverted viscosities are shown in (a), the former using dotted lines and the latter using solid lines of the same colour.The corresponding displacement misfits, given by eq.(3.28), at each iteration of the inversion are shown in (b).The earthquake source and starting structure for the inversion are as described in Section 2.5.