zotero/storage/TP3TLKER/.zotero-ft-cache

High Resolution Modelling of Post-Glacial Rebound and Sea-Level Change
Anthony Patrick Purcell
A thesis submitted for the degree of Doctor of Philosophy of the Australian National University

The work in this thesis was carried out while I was a full-time student at the Research School of Earth Sciences at the Australian National University, between March 1991 and April 1997. Except where mentioned in the text, the research here is my own and is believed to be original. No part of this thesis has been submitted to any other university or similar institution.
thony Purcell, September 1997.

Acknowledgements
I would like to take this opportunity to thank those who have made this thesis possible. The patience, advice, encouragement and support of my long-suffering supervisor, Prof. Kurt Lambeck, were invaluable, and I am deeply indebted for his insight and perseverance. During the course of my studies I was fortunate enough to share offices with Dr. Paul Johnston and Dr. Dan Zwartz, both of whom share an extraordinary enthusiasm for science and a flair for keeping life interesting, I am very grateful to both for their comments and support, which were freely given and always constructive. Special thanks also goes to Dr. Claudine Stirling for her encouragement and for being there to vent to, even though she well and truly beat me to the printers. I owe many thanks to Herb McQueen, Jean Braun, Cathy Smither, Georg Kaufmann, Malcolm Sambridge, Brian Kennett, and Patrick Wu, all of whom took the time to help me wrestle with problems theoretical and practical. Vivien Gleeson, Clementine Krayshek and Janine Dolton helped with administrative detail and encouragement, and saved me more angst than I care to contemplate. The companionship and support of my fellow students is also greatly appreciated, particularly Melita Keywood, Adam Kent, Leesa Carson, Julie Quinn, Yusuke Yokoyama, and Kevin Fleming.
Outside of the school, the understanding and tolerance of my family and friends has been nothing short of outstanding, for all they have given and have endured I extend my warmest thanks. In particular, I have been extremely fortunate to count Anthony Bushell, Andrew Hide, Lisa, Alyson, Gwen Morse, Neville Smythe, Brett Evill, Kenneth Beaton, and Takako Toda-Knight among my friends, and each is in no small part responsible for the completion of this work. Finally I thank God for all that he has granted, most especially the love and support of my parents, my grandmother, my sister Kerrie, and my brothers, Michael, Thomas and Paul.

ABSTRACT
The history of deformation of the earth's surface due to the removal of the Holocene and latest Pleistocene ice sheets is recorded in the global sea-level record. Recent advances in collection and dating of sea-level data, and computing technology have made possible high-resolution, high-precision numerical modelling of deglaciation events. Currently available mathematical formulations for the earth's response to surface loading are ill-suited to such high-resolution modelling however, suffering from numerical instability or excessive computational cost. In this thesis, the suitability of established techniques for modelling short wavelength surface load problems is closely examined. Achieving high resolution with the conventional spherical harmonic scheme has previously been shown to be prohibitively expensive computationally, the numerical instability or physical inappropriateness of associated harmonic analysis procedures is demonstrated here. The earth is then modelled as a flat, semi-infinite half-space and a new formalism developed that is exceptionally stable at depth, based on the wave propagation technique of seismology. This formalism is extended to include the effect of pre-stress and dilatation and their relative effects considered in detail. In the visco-elastic case, pre-stress advection is demonstrated to be critical in forcing deformation toward equilibrium with the load, and dilatation is shown to be largely negligible provided prestress is included. The new formulation is then tested numerically to establish its stability as part of a large super-position problem. The stability of numerical inversion schemes for the Fourier and Laplace transforms is also established, and their impact on the implementation of the wave propagation technique found to be small provided sufficient care is taken in choosing the inversion procedure to be used and their various input parameters.

Contents

Introduction

1

1 Deformation of a Spherical Body by Surface loading

1.0 Introduction

8

1.1 The Field Equations

10

1.1.1 The Initial State

11

1.1.2 The Incremental Field Equations

11

1.2 Boundary Conditions

14

1.3 Constitutive Equations

16

1.3.1 Constitutive Equation for an Elastic Body

16

1.3.2 Constitutive Equation for a Maxwell Visco-Elastic Body

16

1.3.3 Integral Transfarms

17

1.3.4 Applying the Correspondence Principle to a Maxwell Fluid

18

1.4 Properties of Spherical Harmonic Functions

20

1.4.1 Spherical Harmonic Functions

20

1.5 Calculating the Deformation of an Elastic Sphere

23

1.5.1 The Spherical Harmonic Formulation for an Elastic Solid

23

1.5.2 Spherical Harmonic Formulation for a Liquid Core

27

1.5.3 Boundary Conditions

28

1.5.4 Love Numbers

31

1.5.5 Loads ofDegree 1 or 0

32

1.6 Spherical Cap Harmonic Analysis

34

1.6.1 Solving Laplace's Equation

34

1.6.2 Theory ofSturm-Liouville Boundary Value Problems

36

1.6.3 Resolution ofSurface Features

38

1.6.4 Spherical Cap Harmonic Analysis

39

1.7 Discussion

46

2 Deformation of a Stratified Semi-Infinite Elastic Body

2.0 Introduction

47

2.1 Formulation of the Problem

49

2.2 Basic Propagator Matrix Analysis

50

2.2.1 Solving the Equilibrium Equations

51

2.2.2 Application to a Stratified Body

53

2.2.3 Boundary Conditions

54

2.2.4 Discussion

56

2.3 Conventional Propagator Matrix Analysis

57

2.3.1 Calculating the Exponential ofa Matrix

59

2.3.2 Application to the Problem ofan Elastic Medium ·

63

2.3.3 Application to a Stratified Body

66

2.3.4 Boundary Conditions

66

2.3.5 Discussion

68

2.4 Wave Propagation Analysis

69

2.4.1 Reflection and Transmission ofDeformation Fields

72

2. 4. 2 Composite Reflection and Transmission Matrices

7 6

2.4.3 Application to a Stratified Body

78

2. 4. 4 Introducing the Source ofDeformation

82

2.4.5 Boundary Conditions for a Buried Source ofDeformation

84

2.4. 6 Buried Sources and Receivers

88

2. 4. 7 Su,face Loads and Receivers

91

2. 4. 8 Discussion

92

2. 5 Numerical Stability of Propagator Matrix Techniques

93

2.5.1 Su,face Deformation due to a Square Load: Analytical Solutions 93

2.5.2 Numerical Stability ofPropagator Matrix Techniques

94

2. 6 Discussion

99

3 Numerical Implementation of Integral Transform Techniques

3. 0 Introduction

100

3 .1 The Fast Fourier Transform

102

3. 2 The Laplace Transform: Numerical Evaluation of the Mellin Integral

107

3.2.1 Comparison ofNumerical Techniques for Evaluating

the Mellin Integral

113

3. 3 Discussion

120

4 The Effect of Pre-Stress Advection and Dilatation

4.0 Introduction

122

4.1 Including Pre-Stress Advection and Dilatation in the Flat-Earth

Formulation

124

4.1.1 Incorporation into the Propagator Matrix Technique

124

4.1.2 The Importance ofPre-Stress Advection in the Elastic Case

129

4.1.3 Pre-stress Advection and Dilatation in the Visco-Elastic Case 132

4. 2 The Case of an Incompressible Half-Space

141

4. 3 Discussion

143

5 Conclusion

144

Bibliography

Appendix A: Propagator Matrix Entries

I Introduction

11111
Page 1

INTRODUCTION

This thesis is concerned with theoretical aspects of the earth's response to surface loading on time scales of 103 yr to 105 yr . The motivation for this investigation is to
provide a high accuracy formalism for modelling the Earth's response to cycles of glacial loading and unloading, from which the Earth's viscosity structure and response parameters may be inferred. A secondary motivation, not examined in this thesis, is the development of a high resolution inversion scheme to infer a detailed history of the deglaciation process from observations of the associated rebound.
Our understanding of the physical properties of the earth is severely limited by the fact that direct measurements can only be made at or near the earth's surface. Careful observation in a diverse range of fields has however made it possible to infer constraints on many of the processes at work in the earth's interior and their various interactions, though there is still much that we do not know with any certainty.
The earth's response to loading cycles of high frequency (less than a few years) is well established, the earth behaving essentially as an elastic body with parameters inferred from the analysis of seismic travel time data. At these frequencies the earth is shown to have a liquid core of radius 3 480 km with a semi-solid 1 250 km inner core, and to undergo a rapid change in elastic parameters from a depth of 410 km until a sharp discontinuity at 670 km marks the boundary between upper and lower mantle.
Apart from seismological phenomena, the most significant body of available observational data comes from the interaction of the continental plates. Plate tectonics is the largest scale geophysical phenomenon whose effects are directly observable at the surface and is therefore fundamental to our understanding of the earth as a whole. Any theory concerning the physical processes at work in the earth's interior must be consistent with the observed behaviour and properties of the continental plates.
On the time-scales of continental drift ( ~ 106 yr ) much of the earth may be treated as
a fluid under an elastic or visco-elastic lithosphere, the motion of the continental plates being attributed to convection of the mantle material beneath. The nature of this mantle convection is of great importance to the field of earth science and is very strongly dependent on the viscosity of mantle material. If it were possible to determine the rate of mantle convection and whether it takes place throughout the entire mantle or occurs separately in the upper and lower mantle we would be able to place important constraints on the properties and behaviour of mantle material, such as the degree of mixing that occurs throughout the mantle, temperature distribution, the origin and composition of mantle plumes, the nature of hotspots, and the genesis of material from mid-ocean ridges. These properties themselves have far reaching consequences in a wide range of

Introduction

Page 2

geophysical and geochemical researches from the nature and behaviour of the earth's magnetic field through to the origin of the moon.
Mantle viscosity is constrained by the response through time of the earth under a surface load. The deformation that results form this loading is a function of both the rheology of the earth and the nature of the load, so that if the loading history is well understood then the relaxation process provides important data on the properties of the earth over the time-scale of the load. The results of such studies may not be directly applicable to the problem of mantle convection which takes place over particularly long
time-scales ( ~ 106 - 107 ), but it may be possible to extrapolate such data to longer time-
scales. At intermediate frequencies the earth cannot be satisfactorily modelled as either an
elastic or fluid body and a theoretical treatment is required that spans these two limiting cases. This thesis attempts to develop an appropriate mathematical formalism, addressing in particular, methods that have the potential for very high spatial resolution. These techniques may then be applied to glacial rebound data to constrain both the history of the deglaciation process and the visosity profile of the earth.
A large body of high-quality rebound data is preserved in the sea-level record of the Holocene and Late Pleistocene. Geomorphological evidence indicates that over the last 2 million years the earth has been undergoing periodic glaciation and de-glaciation. The glaciation process seems to take place over a period of approximately 90 000 years and is followed by a rapid deglaciation which is complete after about 15 000 years. There then ensues an interglacial period of around 15 000 years before the cycle is repeated.
During the last glacial maximum, which occurred approximately 23 000 years ago, there were substantial glaciers over North America, Fennoscandia, Antarctica, Great Britain and the Barents Sea. These glaciers represented extremely large volumes of ice, the Laurentide ice sheet for example was about 3 000 metres thick and covered a large portion of North America. The melting of the glaciers and concomitant addition of meltwater to the oceans therefore resulted in substantial changes in the moment of inertia of the earth, and the shape of both the geoid and the earth's surface.
When a large load is placed on the surface of the earth it will displace material in the mantle and so deform the earth's surface and gravitational field . If the distribution of the load is subsequently altered the displaced material will move towards equilibrium with the new load distribution, resulting again in a change in the shape of the earth's surface and
the geoid. The relative position of the earth's surface and the shape of the geoid at times
during and since the last period of deglaciation are preserved in the global sea level record and geoid observations. Direct geoid measurements are only available for a small fraction

.,4

I

Introduction

111111
Page 3

of the period since the last glacial maximum so that the vast majority of detailed postglacial rebound data are obtained from observations of sea level.
The change in sea level around the globe resulting from the melting of latest Pleistocene ice sheets is a function of both the change in load and the rheological properties of the earth. Sea-level is measured relative to the position of the shoreline so that the effect of all up-lift and subsidence events will be preserved in sea-level signals. For many sites it is possible to remove the contribution to sea-level change due to tectonic processes and other loading events however, and observations from such sites may then be used in conjunction with the theory developed for the problem of glacial rebound and sea-level change to constrain both the rheological properties of the earth and the melting
history of the ice sheets. The theoretical development for the deformation of the earth under surface loading
is quite involved but has been extensively developed by a number of different investigators (Haskell 1935, McConnell 1965, Peltier 1974, Cathles 1975, Farrell & Clark 1976). In the simplest case sea level is altered simply by the volume of meltwater that is added to the oceans, this so-called eustatic sea level change is quite a large effect contributing 105-140m to global sea level depending on the ice model used. Farrell & Clark (1976) demonstrated that even if the earth is assumed to be completely rigid, postglacial sea level change will be affected by changes in the gravitational attraction of the ice sheets and gravitational self-attraction of the water in the ocean basins.
If we allow the earth to deform under surface loads then calculation of sea level change is further complicated by the earth's response to the change in distribution of the surface load as the ice sheets are removed and water is added to the ocean basins (Bloom 1967). The change in geometry of the ocean basins due to the volume of meltwater added and loading effects result in the integral equation for sea level change being implicit. That is, sea level change at a particular point is dependent on the amount of sea level change at that point so that an iterative procedure is required to find a solution (Peltier & Andrews 1976). Johnston (1993) showed that for realistic earth models the effect of the ice load on the geometry of the ocean basins is accurately determined after two iterations.
A change in surface load results in a redistribution of mass within the earth to achieve equilibrium with the new load. The resulting change in the earth's gravitational field must be considered since both the mean shape of the surface of the ocean (which will be an equipotential function) and the deformation caused by a given surface load will vary if there is a change in the earth's gravitational field (Farrell & Clark 1976).
Various models for the earth's rheology have been considered, their complexity increasing as data, computational resources, and theoretical development allowed. Haskell (1935) and others modelled the earth as being completely viscous while Love (1911, 1927) and Longman (1962) developed the theoretical techniques necessary to

Introduction

Page 4

model the deformation of an elastic sphere. McConnell (1965) considered the deformation of a stratified, semi-infinite half-space with Maxwell rheology and Peltier (1974) developed the formalism for a radially symmetric spherical Maxwell body.
Having developed the theory for the deformation of a self-gravitating visco-elastic sphere under a surface load past investigators have used the resulting forward modelling procedures to invert post-glacial rebound data and constrain both the rheology of the earth (e.g. Cathles 1975, Peltier 1976, Nakada & Lambeck 1987, Gasperini & Sabadini 1989, Mitrovica & Peltier 1991, Johnston 1993) and the melting history of the Holocene ice sheets (e.g. Lambeck 1993a,b, Quinlan & Beaumont 1982).
Determining the melting history of the ice sheets since the _last glacial maximum is a difficult problem in its own right. The extent of the glaciers as a function of time can often be determined from various geomorphological indicators but the thickness of the ice sheets cannot usually be determined directly except at a few isolated points where nanatuks were formed. Isotopic records broadly indicate the total amount of water that was removed from the oceans to form the ice sheets while palaeoclimatic data give some indication as to when and where the climate was most favourable for the formation and thickening of glaciers, there are of course also mechanical constraints on the shape of the glaciers but still many details of the melting history of the holocene ice sheets are relatively poorly known. For example, the vast majority of the land on which the Barents sea ice sheet rested is now submerged, obliterating much of the geomorphological evidence of its extent and progress, this is also a problem in many other areas such as the Irish sea, the North sea, the gulf of Bothnia, and the coastal fringes of the Laurentide and Antarctic ice sheets. Even where the geomorphological evidence is intact it is usually only possible to determine the direction and time of transit of the most recent ice sheet, all evidence of previous ice sheet transitions being overwritten.
The rapid improvement of the distribution and quality of sea-level observations throughout the Holocene and latest Pleistocene has made detailed modelling of ice sheets possible for the first time (Lambeck 1993b, Lambeck et. al. 1996). The vast majority of the formalism developed for the problem of glacial rebound has however been devoted to forward modelling, so that attempted inversions of the sea-level record to constrain the deglatiation process have been largely restricted to an iterative procedure of forward modelling and adjustment of the ice model. The detail that may be achieved with such a scheme is limited by the resolution of the modelling procedure being used which in some cases has proven inadequate.
Throughout the British Isles the sea-level record is generally of very high quality, particularly on the east coast of Scotland where past shorelines are preserved along the shores of firths and river valleys, in particular the Forth estuary, the valleys of the Tay and Earn Rivers, and the firths of Beauly, Cromarty, Inverness and Moray. In this

..4

Introduction

Page 5

region the indicators of past sea-level mostly take the form of preserved shorelines extending substantial distances along inlets and estuaries, allowing very accurate
determination of gradients. One feature of particular interest is the relationship between the Main Late-glacial,
High Buried, Main Buried and Low Buried shorelines. Age constraints from pollen and carbon-dating place these shorelines at approximately 10 300, 10 100, 9600, and 8700 years BP respectively (Haggart 1986, Firth & Haggart 1989) and their elevations indicate that sea-level in this region oscillated over that period. The nature of these oscillations is uncertain since available sea-level observations from other regions during this period are dominated by either addition of meltwater or post-glacial uplift, the sheer magnitude of these contributions is enough to completely mask any small amplitude features so that it is impossible to judge whether this is a global, regional, or local feature.
One possibility is that the oscillations are due to rebound associated with the Loch Lomond stadia!, an extended climatic deterioration from 12 000 years BP onward during which a small ice sheet extended back into the upper valleys of the Forth and Tay rivers and down into the Clyde estuary (Sissons 1974, 1981). The readvance seems to have reached its maximum extent between about 11 000 and 10 500 years BP before retreating once again to vanish completely by 9 000 years BP.
It is equally possible however that they are a feature of the global sea-level, and if this is the case, their nature and behaviour could be significant in understanding the climate of the late Pleistocene and the mechanism of sea-level change. To distinguish between these two possibilities requires a technique for accurately modelling the rebound due to the ice sheet over these time periods. The lateral extent of the ice sheet in this period is significantly smaller than for earlier phases of the deglaciation process and attempts to model its effect have so far been unsuccessful due to poor resolution of the global modelling procedure for surface features of this wavelength.
The purpose of the current work is to investigate the various techniques by which rebound data may be used to infer details of the surface load at very high resolutions, and particularly techniques appropriate to incorporation into a direct inversion scheme.
The power and flexibility of the spherical harmonic scheme of Cathles (1975), and Farrell & Clark (1976) is such that it should not be completely discarded without first examining the possibility of extending it to the high-resolution case. Chapter 1 reviews the formalism of loading problems and particularly the analytical techniques developed for the case of a spherical body. The performance of this global spherical harmonic scheme is generally poor when applied to short wavelength loads however, since the power spectra for earth models are heavily biased toward harmonics of low degree, and high resolution features such as the Loch Lomond readvance require extremely high degree approximations to achieve even moderate definition.

Introduction

Page 6

Spherical harmonic analysis is a common tool in modelling the potential of the earth's magnetic field and several techniques (Haines 1985, De Santis 1991, 1992) have been developed to extend the global formalism to regional scale problems. The formalism of these techniques will be developed and their suitability to high resolution rebound problems closely examined.
For short to medium wavelength features such as the Fennoscandian and British ice sheets the sphericity of the earth may be neglected and it may be modelled as a flat semiinfinite body (Wolf 1984, 1985a, Amelung & Wolf 1994). Chapter 2 details the analytical development of three distinct techniques for modelling the deformation of a flat, stratified, elastic half-space in three-dimensional Cartesian coordinates.
Theoretically, the results of these techniques are equally valid but in practice the formalism of propagator matrix procedures is very susceptible to numerical instability with depth. This is particularly undesirable if the technique is to be used to model a larger ice sheet through superposition since the larger ice sheet will stress the earth to far greater depths than the component loads. Chapter 2 closes with a discussion of the relative analytical and numerical properties of each of the flat earth procedures, and a comparison of each against the analytical results for a uniform half-space.
Propagator matrix procedures rely implicitly on the use of Fourier transforms, and the extension to the visco-elastic regime is achieved with the use of the Laplace transform. The final results of flat earth modelling procedures are therefore obtained only after numerical inversions from the Fourier and Laplace transform domains. Numerical inversion of integral transforms carries significant risk of computational error and inversion routines and associated parameters must be chosen carefully to guarantee that the accuracy and stability of the propagator matrix procedures are not undermined by either of the inversions required.
The Fourier transform is a powerful and widely used numerical tool that is generally well understood but requires particular attention when applied to propagator matrix procedures of this type. In contrast, inversion of the Laplace transform is a very complex issue, particularly for compressible earth models, and though the collocation technique of Schapery (1962) provides good results for a large variety of earth models, it is prone to introduce small amplitude oscillations over long time-scales and there has previously been no means of comparison for its results for rheologies where the normal analysis scheme fails (which includes a large class of plausible earth models).
Chapter 3 examines the implementation of both Fourier and Laplace transform inversion schemes and possible sources of error for each, as well as providing several means of comparison for the collocation and normal mode techniques for compressible and incompressible earth models.

Introduction

Page 7

The analytical complexity of the propagator matrix techniques developed in chapter 2 is initially reduced by neglecting pre-stress advection and dilatation. These terms are indispensable if the results of our modelling are to have any physical validity however, and are reintroduced into the forrnalism for the propagator matrix procedures in Chapter 4. The effect of each term singly and combined is considered in the visco-elastic regime, the analysis is also extended to the the incompressible case, and finally the procedure is tested for stability as part of a superposition scheme.
Chapter 5 summarises the results of the previous chapters and discusses their consequences and potential applications. The chapter and thesis close by considering possible refinements that may be made in future implementations to the techniques developed here.

,4

Deformation of a Spherical Body by Su,face Loading

Page 8

Chapter 1
DEFORMATION OF A SPHERICAL BODY BY SURFACE LOADING

1.0 Introduction Deformation of the earth's surface due to loading by the ice sheets and changes in
loading of the ocean basins is an important factor in calculating the sea-level change associated with the exchange of ice and water volumes. Accordingly, an extensive formalism has been developed to deal with this problem. For loads of small lateral extent, however, this global modelling procedure may not be particularly appropriate or efficient and it may be necessary to either adapt the procedure or develop new ones to accurately model the high-resolution case. In this chapter, I will briefly review the relevant formalism of continuum mechanics and the global modelling procedure itself with particular attention to its applicability to small-scale problems, and introduce the Spherical Cap Harmonic Analysis (SCHA) procedure developed by Haines (1985) and the Adjusted Spherical Harmonic Analysis (ASHA) scheme of De Santis (1992). Although these procedures were initially developed to model the regional behaviour of the earth's magnetic potential they may both be adapted and used with the global forrnalism to model local deformation.
The first step in developing a mathematical formalism for the deformation of the earth is to infer some of the a priori physical properties of the earth at the time scales we will be considering. The most obvious of these properties is of course the earth's shape. Neglecting ellipticity, topography, and diurnal rotation we may consider the earth as a static, spherical body and model its response to surface loading accordingly once we have determined the nature of its rheology.
Seismic travel time data indicate that the earth's elastic properties vary with depth and that there are layers within the earth throughout which these properties vary smoothly and some depths at which the physical properties change quite sharply (usually thought to correspond to chemical boundaries or mineralogical phase transitions). For mathematical convenience we will neglect lateral heterogeneities and assume that the various rheological parameters are a function only of depth. This allows us to model the earth as a stratified sphere whose rheological properties are prescribed throughout a series of layers. In some formulations the gradients of density and the elastic moduli are taken to be zero which is obviously not the case in the real earth although any function of depth may be closely approximated by a step-function if the layering is taken to be fine enough. Neglecting the local density gradient has a number of consequences for the validity of our modelling procedures and the interpretation of our results (see for example Kaufmann 1991 ,

,4

Deformation ofa Spherical Body by Surface Loading

Page 9

Johnston 1993) but a detailed discussion of these considerations is beyond the scope of the current work.
On seismic time scales the crust and mantle behave as elastic solids while the outer core acts as a liquid and the inner core as some sort of slurry. On longer time scales (of about 106 yrs) the mantle seems to convect as a liquid. Rebound data indicate that for deformations over glacial time scales (of 103 - 105 years) there is an initial 'elastic'
component of deformation fallowed by an extended period of viscous creep indicating some sort of viscoelastic rheology (Lambeck, Smither & Johnston 1997). In fact we will assume that the earth (with the exception of the fluid core) is a Maxwell viscoelastic body since for such a body the constitutive equations have the simplest possible form consistent with the properties discussed above. Petrological and laboratory evidence indicates that mantle material actually has a more complicated non-linear rheology (see for example Ranalli 1982) but model results obtained assuming the earth to be 'effectively' a Maxwell viscoelastic body have so far proven sufficiently accurate to obviate the need to consider more complicated models (see for example Lambeck 1993a,b, Peltier & Andrews 1976, Wu & Peltier 1983).
Mindful of these assumptions, in this chapter we will review the theory for the deformation of a radially stratified, spherical, Maxwell viscoelastic body under surface loading from its foundations in continuum mechanics, paying particular attention to the Love number derivation (bearing in mind that the SCHA technique uses associated Legendre functions of non-integral degree). We will then discuss the global technique's applicability to short wavelength problems and perform a similar derivation and analysis of the SCHA and ASHA procedures.

~

Deformation of a Spherical Body by Su,face Loading

Page 10

1.1 The Field Equations For the moment we consider a general body subject to some surface or body force.
We will assume that the resulting deformation is small and occurs slowly and steadily, ensuring the system is in equilibrium at all times. We will also assume that the system is closed to any other external thermodynamic contributions and that there are no chemical changes in the system itself.
Under these assumptions the continuum form of Newton's second law (otherwise known as the inertia equation) is:

V . t + p<O)f =p<O)ii

(1.1.1)

where the (0) superscript indicates the initial value of the relevant function (in this case density but this convention will be adopted throughout for other functions) prior to the application of the force under consideration. The double dot over the u term indicates second derivative with respect to time.
Equation (1.1.1) is the Lagrangian form of the inertia equation, where t is the
Piola-Kirchoff stress tensor, f is the body force per unit mass, p is the density, and u is the displacement. We will restrict our attention to the Lagrangian formulation rather than the Eulerian since it is most easily adapted to the modelling of glacial loading problems.
These equations are valid for a reference frame that rotates with the earth. The body forces in this instance include gravitation, centrifugal force, and the Coriolis force. We will neglect the latter two effects since diurnal and annual contributions are unlikely to have a significant impact on deformation occurring on a glacial time frame (c. 103 - 105 yrs), so that we are effectively considering a stationary sphere.
The second of the field equations is the continuity equation which enforces conservation of mass within the earth:

·v-(pu) = 0

( 1.1.2)

The last of the relevant field equations is that governing the earth's gravitational field. Gravitational acceleration is taken to be the gradient of the gravitational potential </J which satisfies the equation:

V 2 ¢=-4nGp

(1.1.3)

where G is Newton's gravitational constant.

~

Deformation of a Spherical Body by Su,face Loading

Page 11

§1.1.1 The Initial State We assume that the body we are considering is initially in hydrostatic equilibrium.
This is not entirely the case in the earth today as is demonstrated by the existence of gravity anomalies, but should be a reasonable assumption for the most part. If we also assume that initially the only force acting on our body is self-gravitation then the initial field satisfies the relations:

lO) = - I p(O)
V p (O) = p (O) V q/0)
V 2 </Jco) = - 4rcGpC0)

( 1.1.4a) (l.1.4b) (l.l.4c)

where pis pressure and I is the Kronecker delta tensor (see Wolf 1991).

§1.1.2 The Incremental Field Equations Rather than modelling the total response of the body to a force it is useful to
consider only perturbations or increments from the initial state. The various properties of a body may be expressed as functions of either current position (the Eulerian formulation) or initial state (the Lagrangian formulation). It is easier to express boundary conditions using the Lagrangian framework but both are capable of expressing perturbations in the quantities we are interested in.
We will adopt the incremental notation of Wolf (1991) given in table 1.1 below:

Property particle position displacement function initial value material increment local increment

Lagrangian Formulation
X
z(x,t)
u(x, t) =l(x, t)-x
f(x, t)
fc0\x) =f(x , 0) fUi>(x) =f(x, t)- f c0)(x) f CL1)(x) =f(x, t)- fc0\z (x, t))

Eulerian Formulation x(l, t)
l
u(x, t) u(x, t) =l -x(l, t)
F(l, t)
F c0)(z) =F(z, o) p (o>(z) =F(l, t)- p(O)(x(l, t))
p(L1)(z) =F(z, t)- pc0)(z)

Table 1.1: Notation for local and material increments in Eulerian and Lagrangian forms.

Deformation ofa Spherical Body by Su,face Loading

Page 12

Given a function f the material increment f c8>(x) is the change in the function for the particle X while the local increment f c.1)(x) is the change in the function at the
position currently occupied by X . Assuming infinitesimal perturbations the material and local increments may be related using a Taylor series expansion:

f (O) = f (L1) + u . V f (O)

( 1.1.5)

which is accurate to first order. Substituting this expression into the inertial and gravitational equations ( 1.1 .1) and
(1.1.3), and noting that for the problem of glacial rebound the acceleration term on the right hand side of the inertia equation may be neglected, yields the incremental forms of
the equations:

u) V. t(O) + V ( p (O) g (O) . + p (L1) g CO) + p <O)g<t1) = 0

(1 . 1.6)

V 2 </Jct.) = - 4nGp<t1)

( 1.1.7)

where g is the acceleration due to gravity and is given by:
g=V<p

( 1.1.8)

and is taken to point upwards. The local increment in density is given by the incremental form of the continuity
equation (1.1.2):

u) p <L1) = - V .(p<0)

(1 . 1.9)

This term assumes a particular convenient form inside a region of uniform density where it becomes the product of the original density field and the divergence of the deformation (also called the dilatation). I will hereafter refer to the component of (1.1.6) containing
p(t1) as the internal buoyancy term.
The second term in equation (1.1.6) is the pre-stress advection term which represents the fact that before the current stress field was applied the body under consideration was already under stress due to its own weight as a result of its hydrostatic initial condition. The pre-existing stress field (the pre-stress) advects with the deformed material within the body and is super-imposed by the new stress field.
The pre-stress and internal buoyancy terms are both crucially important in assuring the physical meaningfulness of our modelling, particularly for deformations of the scale that typically occur in analyses of glacial rebound. The nature and magnitude of their

Deformation ofa Spherical Body by Su,face Loading

Page 13

effect will be discussed extensively when we come to model the deformation of a flat semi-infinite half-space in chapter 4.
The incremental field equations take a simpler form inside a liquid since an inviscid material cannot support any shear stresses and there is no dilatation since a fluid can instantaneously move to dissipate any increase in pressure (except in the case where the increase is hydrostatic). Since there is no dilatation there will be a local density gradient (radial in the case of a sphere, vertical in the case of a flat semi-infinite half-space). From (1.1.9) we see that for a spherical body the local change in density is given by:

=- - - = 0p(O))

P(L1)

( dr U -p (O)' U

r

r

(1.1.10)

where we have used the prime(') to denote differentiation with respect tor. Using (1.1.6), (1.1.7) and (1.1.10), noting that there is no dilatation, and setting
t <O'J = 0 yields the incremental field equations for an inviscid body such as the core:

V-u =0

_v( U) + p (O) g <O) .

p (O)' Ur g <O) _ p (O) g CL1) = 0

V 2 q/..1) = 4,rGp(O)' Ur

( l.1.11) (l.l.12) ( l.1.13 )

Deformation ofa Spherical Body by Su,face Loading

Page 14

1.2 Boundary Conditions Throughout an elastic body the absence of cavitation, cracking and slip requires that
the deformation, gravitational potential, and normal stress be everywhere continuous. At a layer boundary, the partition between two different rheological regimes, there may be a discontinuity in the density or elastic moduli of the material. Regardless of any such
inhomogeneity, given the vector normal to the boundary, n, the quantities ui, q/'1) , and
L,i n;°l tij"l , and n<0J - ( V ¢/M - 4;rGp<0lu) must be continuous across the boundary. The
third of these boundary conditions is obtained by integrating the inertia equation (equation (1.1.6)) and letting the thickness of the pillbox approach zero while the last is given by repeating this integration, applying the convergence theorem, using the form for the internal buoyancy term given in (1.1.9), and letting the thickness of the pillbox again approach zero (see Cathles 1975, p 17). Doing so yields:

i l ( 0 = V · (V¢ (6.) - 4n-Gp<0lu)dv = V¢ (6.) - 4n-Gp<0lu) · da
J: =dAn ·[V9<"l -4n-Gp<0lu

(1.2.1)

Letting the height of the pillbox approach zero causes the area of the side of the cylinder to vanish and we get the boundary condition we require. Note that for a phase boundary the boundary conditions are different since (1.1.9) no longer applies the effect of this contribution has been considered for the spherical case by Johnston, Lambeck & Wolf (1996) but is beyond the scope of the present work.
We will assume that these boundary conditions hold at all internal boundaries except the core-mantle and free surface boundaries which we will discuss shortly.
It should be noted that in the case of a stratified sphere, the normals to the layer boundaries are simply radial vectors so that the third of the boundary conditions given above corresponds simply to continuity of the radial components of shear stress, tre and trrp, and the radial stress itself, trr. Analogously, for a flat, stratified half-space the condition guarantees continuity of the vertical components of shear stress, txz and tyz , and the vertical component of stress tzz .
At the core-mantle boundary the situation is complicated (Dahlen 1974, Dahlen and Fels 1978) by the presence of a thin boundary layer inside which the dilatation is non-zero and the radial and tangential components of deformation change rapidly. If this layer is neglected then the deformation appears to be discontinuous across the core-mantle boundary. The perturbation to the gravitational potential and its gradient are continuous and smooth through the boundary layer so that the potential is continuous across the boundary. Since the core-mantle boundary layer represents a chemical transition and is

~

Deformation ofa Spherical Body by Surface Loading

Page 15

therefore governed by (1.1.9), the effective discontinuity in the gravity perturbation is a function of the effective discontinuities in density and radial deformation. The mathematical implementation of these boundary conditions will be discussed in greater detail during the spherical harmonic formulation in section 1.5.
For a flat semi-infinite half-space additional internal boundary conditions are given by a simple physical principle: all perturbations must tend to zero as distance from the source of deformation tends toward infinity.
At the surface we will restrict our attention to the case of an external force applied normal to the surface of the body for which there is no traction, in particular a surface load whose initial force per unit mass would be gf), the gravitational acceleration of the body at the surface (in both the spherical and the flat earth cases this value will be constant over the free surface). Outside the body we also have that p<!J.) = 0 and the local increment in gravitational potential satisfies Laplace's equation (V2 q}i1) = O). Our other surface boundary conditions are given by observing that perturbations to gravity and potential must tend towards zero with increasing distance from the body. The mathematical formulation of the surface boundary conditions for a spherical body will be discussed in more detail in section 1.5.3.

,..

Deformation of a Spherical Body by Surface Loading

Page 16

1.3 Constitutive Equations §1.3.1 Constitutive Equation for an Elastic Body
For deformations with sufficiently small gradients of displacement in a Cartesian coordinate system the various components of strain, Eu , are defined:

Eu =

1
2

(dduxi

+<dIux1; )

1

(1.3.1 )

ui where the represent the corresponding components of displacement (see for example
Fung 1965). The material increment in stress , t (o) , is related to the strain by the constitutive equation. For an elastic solid the constitutive equation has the form:

i(

0 )

=

.,1,(V

· u

)I+

2µE

(1.3 .2)

II

where .,l andµ are the Lame parameters of the body. Noting from ( 1.3.1 ) that

3
V · U = L £kk we may rewrite (l.3.2) in component form: k =l

A( t~o) = ~ Js kk) + 2µ £ii

( 1.3 .3)

3
L where we have adopted the Einstein summation convention £kk = £kk . k= l Setting i = j in (1.3.3) and summing, we see that:
tli)= (3A + 2µ) £kk

(1.3.4)

which shows the effect of the bulk modulus of the body K = A+ 2µ/3

§1.3.2 Constitutive Equation for a Maxwell Visco-elastic Body We wish to model the deformation of the earth which initially has shear strength like
an elastic body but under continuous stress loses its rigidity and behaves something like a Newtonian fluid. For such a material the constitutive equation has the form:

tli)) ( · ) i ~o) + µ7J (t~o) - 8-11-3 = A 8iJ £kk + 2µs· u

( 1.3.5 )

where T] is the viscosity of the fluid. Setting i = j in (1.3.6) and summing yields a system directly equivalent to (1.3.4):

i~) = (31 + 2µ ) Ekk

(1.3.6)

~

Deformation ofa Spherical Body by Su,face Loading

Page 17

§1.3.3 Integral transforms
A useful tool in the solution of differential equations of the form of (1.3.5) is the integral transform. Using integral transforms translates a differential equation into a corresponding algebraic equation which may be solved using standard techniques and the resulting expression inverted to yield the desired function. In the time domain for example, it is often more convenient to mathematically model the behaviour of the Laplace transform of a function in the frequency domain rather than the function itself and then invert to yield the desired result in the time domain, particularly to model the behaviour of a Maxwell fluid (see for example Biot 1965, Peltier 1974). Similarly, in the spatial domain Fourier and Hankel transforms are used to simplify the problem being considered, particularly when modelling the behaviour of an elastic solid (see for
example, Sneddon 1951, McConnell 1968).
The Fourier and Laplace transforms of a function f and their corresponding
inverses are defined (see for example Spiegel 1968):
r r :F(J(x)) =J(v) = f(x)e- ;vxc1x .r-l (J(v)) = 2~ J(vkvxdv = f(x)

f L(i(t)) = ](s) = f(t)e- "dt

(1.3.7)

r .l-

1(J( s ))

=

-1.
21tl

liil!,
T•

Jc-i+Tff

1( s)

est

ds

=

f(

t)

where the constant c used in the bounds of the last integral is chosen to be larger than
Re (sa for all the singular points, Si' of J(s) .
These transforms have the very desirable property of transforming differentiation into an algebraic operation (again see for example Spiegel 1968):

.r(t)=ivJ(v)

L(dt) =s](s)- J(o)

(1.3.8)

where it should be noted that for any incremental quantity f , we havef(0) = 0 .
J , The singular points or poles of the transform function , are the points in the
complex plane at which the function becomes undefined (that is, its magnitude becomes
arbitrarily large as it approaches the pole). A pole si of J is said to be of order mi if mi
is the smallest integer such that }LIT};J(s)(s - si)m; is finite ; it is said to be an essential
singularity of J if there is no finite value of mi for which this condition holds. The corresponding residue, ri , of J is then defined:

~

Deformation ofa Spherical Body by Suiface Loading

Page 18

~ f( )( r1-= I

1

\

• sh••ms;

d mI·-1 (
dsm;- I

~ s

s - s1.) mi)

(1.3.9)

{,

Residues are particularly important in the calculation of the inverse Laplace

transform since the most convenient technique for evaluating the inverse integral given in

(1.3.7) is by closing a path to include all of the poles of the function (assuming that they

are distributed through a bounded region), and then invoking the residue theorem (see for

example Kreyszig 1983). The residue theorem states that the integral of a function around

a closed curve is equal to the sums of the residues of the poles of the function that lie

inside it (this assumes that they are distinct and that none of the singularities are essential).

In the case of the inverse Laplace transform the function est has no poles so all of

the poles of the integrand must be due to J(s). Using a Bromwich path (see for example

Krylov and Skoblya 1977) to close the curve of integration means that the inverse integral

is simply the sum of the residues of the integrand, which are themselves related to the

residues of J(s) (the integral around the closure approaches zero for sufficiently large

values of T by the Jordan Lemma).

§1.3.4 Applying the Correspondence Principle to a Maxwell Fluid Laplace transforming (1.3 .5) and (1.3 .6) yields the alternative form for the
constitutive equations for a Maxwell visco-elastic body:

t~(LJ~)(S ) -

£\t~(s) i

= S2+µsµIn

(~-{ )£11 S

~jt3\is))=

2µ~*(S)(~£11-{S)-

~jt3\ls))

(1.3.10)

r1f)(s) = (3J + 2µ) t\ls)

(1.3.11)

wheres is the Laplace transform variable. Substituting (1.3.11) into (1.3.10) yields the transformed constitutive equation for a
Maxwell fluid:

t~8l(s) = l*(s)( ~jEkk(s)) + 2µ.(s)Eis)

(1.3.12)

which may be re-written:

f(8)(s) = l*(s)(V ·ii(s))I + 2µ*(s)E(s)

(1.3.13)

111
This is directly analogous to (1.3.2) so that the constitutive equation for a Maxwell fluid may be transformed to a form identical to that of an elastic solid with Lame parameters dependent on sand related to the original values by the equations:

...

Deformation ofa Spherical Body by Surface Loading
µ~*(s)= s µ+sµ/ry i'(s)= AS+%(A+2hµ)

Page 19 (1.3.14) (1.3.15)

Comparing (1.3.6) with (1.3.4) we see that the bulk modulus of the material is unchanged and that the constitutive equation for a Maxwell viscoelastic fluid is equivalent to that of an elastic solid with s-dependent Lame parameters.
Numerically inverting from the frequency domain back to the time domain is however not so straightforward in the case of the Laplace transfarm. From the form of the inverse transform given in (1.3.7) we see that a direct inversion is only possible if the transformed function is known to have a particular form for which an inverse has been calculated analytically or if all the singularities of the transformed function and its residues at the poles are known.
From the farm of the Laplace transform in (1.3.7) we see that if the function under
consideration, f(t), is changed dramatically over some appropriately small region then the
value of the transform function, ](s) , will not be significantly affected. The value of the
inverse transform at any point is therefore very strongly dependent on the exact nature of
the transform function, small changes in the behaviour of J(s) may result in significant
changes in f(t) (see for example Krylov & Skoblya 1977). This property of Laplace transforms makes the problem of numerically calculating the inverse transform particularly delicate and every effort must be made to ensure that any numerical error in the calculation of the transform function is as small as possible.
Since the exact form of a function and the nature and position of all its singularities are both rather difficult to determine numerically the problem is often simplified by making a series of assumptions about the form of the transform function (see for example Mitrovica & Peltier 1992). By calculating the value of the Laplace transformed
deformation, ii(s) , of the body for a number of different values of s we can approximate the s-dependence of u(s) and invert to get the original time dependent function u(t). The
problem of determining the deformation of a Maxwell visco-elastic body can therefore be simplified to that of determining the deformation of an elastic solid under a surface load for a range of different Lame parameters and inverting to the time domain. The numerical inversion of the Laplace transform will be discussed in greater detail in section 3.2.

...

Deformation ofa Spherical Body by Su,face Loading

Page 20

1.4 Properties of Spherical Harmonic Functions Mathematically modelling the behaviour of a spherical body is most conveniently
done in spherical polar coordinates (see for example Farrell & Clark 1976). One of the most natural tools for such an analysis are the spherical harmonic functions. We shall
follow the development of Johnston (1993)
I
m

§1.4.1 Spherical Harmonic Functions

The normalised spherical surface harmonics are a set of orthogonal harmonic

I
I
I

functions on the surface of sphere. They arise naturally in the mathematical analysis of

gravitational and magnetic potentials or any conservative field.

I

I
r' j

By separation of variables we see that Laplace's equation in spherical polar

coordinates is

j .

f) f t.f(r, 0, IP)= V2f

1 = ,2

da r

(r2ad!r)

+

1 sin

0

da0

( .
sm

8

ad0

a2

+

1 sin2

0

d<p2

(1.4.1)

where r is radial distance from the origin, 0 is co-latitude, and cp is longitude. This
equation has linearly independent solutions of the form:
fnm(r, 0, <p) = (anm r" + bnm r-(n + 1l)[Anmcos (m<p) + Bnm sin (m<p)] Pnm(cos 0) (1.4.2)

where anm, bnm, Anm, and Bnm are arbitrary constants, and Pnm is the associated Legendre function of degree n and order m (see for example Abramowitz & Stegun 1972) which satisfies Legendre's associated differential equation:

:i) Jx((1

-x2 )

+ (n(n + 1)- 1:~2)y =0

(1.4.3)

The fnm given in (1.4.2) are finite, continuous and single-valued at all points inside a sphere only if all bnm are zero and both m and n are positive integers. In this instance the Pnm are in fact polynomials (see for example Spiegel 1968) and have the property that
for \m \> n , Pnm(x) =0 . Restricting our attention to the case bnm = 0 , m and n integers
Im I such that < n we may define fully normalised spherical surface harmonics:

••1

ynm( e,

'P)

=

Knmpnm( cos

0)

cos sin

lmmclcpp

n>m >O O>m>-n

( 1.4.4)

where the constant term Knm is defined

....

r

Defonnation ofa Spherical Body by Surface Loading

111111
Page 21

I
r( 2
(2n + 1)(2- 8mo) n - Im I+ 1) '

Knm=

r(n+lml+ 1)

( 1.4.5)

By the theory of Sturm-Liouville boundary value problems (see for example Boyce

& DiPrima 1977) these functions satisfy the relation:
r f' d0 sin 0 d({J Yn,mi( 0, ({)) Yn,mi( 0, ({J) = 4n- 8n, n, 8m, m,

( 1.4.6)

which makes them particularly convenient when trying to approximate a function over the surface of a sphere using least squares (see for example Kreyszig 1983).
Another important property of Legendre polynomials (which correspond to the case
m = 0 ) is their generating function (see for example Spiegel 1968). Given two scalar
[-1, quantities, t and x, with x lying in the region 1] , the following equivalence holds

✓

1

1 -2tx

+

t2

00
= n~O Pn(x)tn

(1.4.7)

From the theory of harmonic functions we have that any suitably well-behaved

function (i.e. continuous, differentiable and non-singular almost everywhere) over a

closed, regular region may be approximated to arbitrary accuracy in that region with an

harmonic function (see, for example, Gilbarg & Trudinger 1983). So that over the

surface of a sphere any such function may be written as an infinite sum of surface

harmonics, i.e. there exist coefficients Fnm such that:

r f' q,i)r )irp_ d0 sin 0 d(f) (F~0, q,i)-F(0, = 0

( 1.4.8)

where the functions FN have the form:

N n
Ftv( e, cp) = n~O m~n Fnm ynm( e, cp)

(1.4.9)

In fact using the method of least squares to minimise (1.4.8) for each value of N and
'

1«:

applying (1.4.6) we see that the Fnm are given by the relation:

Fnm=

1
4n

JfoJr

•
d0s1n

0

2
Jfo1t

dcp

Ynm(e,

<p)F(e,

cp)

(1.4.10)

and are notable for being independent of N.

Deformation ofa Spherical Body by Surface Loading

Page 22

A spherical surface harmonic of degree n is defined as any function of the form:

Yn( 0, cp) = mtn cnm Ynm( 0, cp)

(1.4.11)

where the coefficients cnm satisfy the constraint:
n
m£"=.".J-'n cn2m = 1

( 1.4.12)

The gradient of a spherical harmonic function may be written:

ae VYn( 0,

cp) = r0 ay

1p + rsin 0

ayn
cfqi

(1.4.13)

(see for example Spiegel 1968), where 0 and {p are the unit latitudinal and longitudinal
vectors respectively and we have used the fact that surface spherical harmonic functions do not vary radially and so have no radial component to their gradient.
From (l.4.2), (1.4.3) and (1.4.4) we have the following properties of surface spherical harmonic functions:

(r' v2 Yn(0, <p)) = 0

(1.4.14)

si~ l0(sin t~) v2 Yn(0, <p) = ~ 0

0 + 1_ ~2y~

- - n (n + 1) Yn(0, <p)

a (a V2(0oY0n)

= -

n(

n + r2

l)

a Y0 n

+

l
sin 2

e

oY0n + 2cot e0o2qYin)

(1.4.15) (l.4.16)

Deformation ofa Spherical Body by Surface Loading

Page 23

1.5 Calculating the Deformation of an Elastic Sphere In spherical polar coordinates the most general vector field (and in particular the
displacement field for a given surface loading over a spherical body) may be written:

u(x) = f U(x) + VV(x) + f X W(x)

(1.5.1)

where f is the unit vector in the radial direction. The three components of the field are known respectively as the spheroidal, poloidal and toroidal modes. In the case of a radial load with no shear stress the toroidal component of the displacement field must be zero.

§1.5.1 The Spherical Harmonic Formulation for an Elastic Solid Since any finite, differentiable scalar function may be written as a sum of
components of the form of (1.4.2), using the notation of (1.4.9) we may substitute this result into (1.5.1) to show that the displacement field of a spherically symmetric body may be written as an infinite sum of the form:

.'tJ~r u(r, 0, <p) =

[r un(r) Yn( 0, <p) + rvn(r)VYn( 0, <p)]

(1.5.2)

where a is the radius of the sphere being considered, in our case the earth. The functions un(r) and vn(r) are, respectively, the radial and tangential spherical harmonic coefficients of deformation.
Being a suitably well behaved scalar function, the local increment in potential may be written as an infinite sum of the form:

r qi~)(r, 0, <p) =nto (~ </>n(r) Yn( 0, <p)

(1.5.3)

From (1.5.2) and the form for the divergence of a vector field in spherical polar

coordinates (see for example Spiegel 1968) we see that the dilatation, L1 = V· u, may be

written as the infinite sum:

f (L)n l n n(n = L1 = V•u

n -- o a

y n[Un' +

+r 2 Un -

+ 1)
----r'--~ V n

(1.5.4)

where, as in equation (1.1.10) we have used the prime (') to represent differentiation with respect to r .
Similarly, we may substitute (1.5.2) into the spherical polar representations of the various components of strain (see for example Goodbody 1982) to show that they may also be expressed as infinite sums of the form:

Deformation ofa Spherical Body by Surface Loading

to (~r Err=

Yn(u~ + ~ un)

~ =
Ere

l
2

n~O

(L)n
a

ddYen

(urn

+ vn,

+

n -
r

I Vn)

00

L (a) ar ( r ) r Er<p =

1 2 Sin 0 n= O

I..

nci7jni

Un

+

, Vn

n +

I Vn

2
00
E00 =n~ ( ~ )n rI ( Un y n+ Vn 0O(YJ2n)

"111111
Page 24 (l.5.5a) ( l.5.5b) ( 1.5 .Sc) (l.5.5d)

a2r 00

a 1 )

(L y e n=O a r u cC<'p<p -- ."/"-'"J"

n -

n

n + Vn cot

0

o-:\yn

+ S.ilvln2

0

o-:\

m
't'

n
-

(l.5.5e)

~ a2

E0
<p

-_

n.=/-JO

(-ar)n

- rsivnn-0 (- O(f)y- 0n0

-

cot

0

Oyn)
ci7ji

( 1.5 .Sf)

The radial component of the strain tensor, f · £, contains no toroidal component and since the material comprising the body has a linear rheology, the radial component of the stress tensor must also have no toroidal component. The radial component of the material increment in stress may therefore also be written as an infinite sum of spherical harmonic .
components:

JJ~r( f-t(SJ =

Tm(r) Yn(0, <fJ)f + rT11n(r)V Yn(0, <fJ))

(1.5.6)

In practice we are interested only in the problem of calculating the deformation of a Maxwell visco-elastic body which, as we saw in the previous section, may be reduced to that of calculating the deformation of an elastic body for a range of Lame parameter values in the frequency domain and then inverting back to the time domain. From (1.5.5), (1.5.6), and the constitutive equation for an elastic body (1.3.2), we see that the spherical harmonic coefficients of stress must satisfy the relations:

A(u' u- v) T = m

n + n+r 2

n

n(n

r +

1 )

n

+ 2µ(u'n +nrun)

(l.5.7)

T0n = µ (rUn + vn, + n -r I vn )

(1.5.8)

I

Deformation of a Spherical Body by Surface Loading

..... Page 25

Substituting ( 1.5 .3) into the incremental gravitational equation ( 1.1.7), using (1.1.9), and noting that in a stratified body the initial density field does not vary laterally
(so that in the spherical case p C11) = -pC0) V · u - p(0)' ur ), yields the following second order
differential equation governing the spherical harmonic coefficients of gravitational
potential:

1) / 2 1) + 2(n

r ,+('+_______;_-~,+(

'f' n

'f'n

=4nGp(0)

u' n

p(o)' +p-(O)un

+

n+r

n(n +

u n

-

-

-r -

v

n

(1.5.9)

This may be converted into two first order differential equations by introducing a new quantity, q (see for example Longman 1962, 1963, Peltier 1982 and Wu & Peltier 1982), associated with the perturbation in gravitational acceleration. The form for this perturbation is given by the incremental form for the equation for gravitational acceleration (1.1.8) and may be expressed in terms of spherical harmonic components:

g C~) =f · V q/~) =nio (~r gn(r) Yn(0, (p)

(1.5.10)

The form for the spherical harmonic components of the related quantity q are then:

= = qn

¢r n

-

g n

4

n

G

p(0)

u n

,+( 'f'n

+

n

+ r

l

'i+f'in -

4

n

Gp(0)

u n

(1.5.11)

which, when substituted back into (1.5.9) yields:

q'n

=

n(n r+

1) (¢rn

-4nGp(o)v) n

+

4nGp(o)

urn

-

n+r

1 q n

(1.5.12)

The form for this equation is different from that given by Longman (1962, 1963), and

(~r Wu & Peltier (1982) because a factor of

has been removed from all spherical

harmonic coefficients. Applying the general form for the divergence of a tensor in spherical polar

coordinates (see, for example, Goodbody 1982) to the material increment in stress given

in (1.5.6), using the identities (1.5.7) and (1.5.8), and applying the properties of

spherical harmonic functions given in (1.4.11-14), results in:

f.(v ./8))= t (L)nY n = o a

IT' +nT -
n m rm

n(n+l)T

r

0n

+

4µ(n-l)
..'-?

Un

+

-4µr

,
U n

+

2µn(n+l)
-----=-.,-.2-___:_

V n

(1.5.13)

Deformation ofa Spherical Body by Surface Loading

Page 26

t (L)ndYn ))= 2

{J.(v./ r 8

T' +n+3T Tm 2µ(1-n) 2µ, 2µ(n +n-I)

n =o a d 8 en

r

0n + r +

..,2.

Un -

Un -

.2

V n

(1.5.14)

These expressions may be substituted back into the incremental form of the inertial equation ( 1.1.6) to give:

1) 4 1) Tr'n

+

nr y m

-

~ n (n-r+-

T 0n

+µ-(~ n-,2-

-

u n

+4-rµu

'
n

+p(0) ( q n

+4nGp(0)u n

_

<Pi )
r__.!:

f I) 0 + p(o) un( g(o) - g(o)') + n(nr+ vn( 2f - g(o) p(o)) =
(1.5.15)

r r = , n + 3

Trn 2µ( l - n) 2µ , 2µ(n 2 + n - l) p(o) g(o) p(o)

T0n + r T 0n + +

r2

Un -

Un -

r2

V n - r Un + -;::- <pn 0

(1.5.16)

where we have used equation (1.5.11) to substitute for the perturbation in gravitational acceleration, and g(o) is the initial value for the magnitude of the acceleration due to gravity which we may see from (1.1.8) is given by g(0) = - <p(o)' for a radially symmetric earth. This may be substituted into the initial gravitational equation (1.1.4c) to show that:

= g(o)' 4nG p(o) _ 2 g(o)

(1.5.17)

This expression may in turn be substituted back into (1.5.15) and (1.5.16) to produce a

closed set of differential equations comprising (1.5.7), (1.5.8), (1.5.11), (1.5.12),

(1.5.15), and (1.5.16).

This form for the system of equations is almost identical to that given for the fully

non-adiabatic case by Johnston (1993), whose derivation we have largely followed, and

is analogous to that given by Longman (1962).

The system of equations may be re-written in Runge-Kutta form by defining the

vector:

r Yn = ( Un, vn, Tm' Ten, 1/>n, qn

(1.5.18)

which satisfies the differential equation:

Deformation ofa Spherical Body by Su,face Loading

Page 27

dyn = An(r)yn
dr

(1.5 .19)

The entries of An may be calculated from the formulae given above:

-

2A
/Jr

-

r n

1 r

An(n+ 1)
/Jr 1-n
r

1
/3

0

0

0

0

1
µ

0

0

r- 4 4p(D)g(o) n(n +I) (p(o)g(o) r-2r) - 4µ - n n(n +I)

A -j ,2

r

,2

/Jr r r

n -

(o) (o)
Pg;-

2 r

~

( n(n+l)(r+µ)-2µ

)

-fr

_n~3

p(o) r p(o)
- -r

-p(o)
0

4nGp(0)

0

0

0

-

n+l r

1

4nGp(0) r

n(n + 1)4nG p(0)

-

r

0

0

n(n +
r

1)

-n-+r l

J

(1.5 .20)

r= where f3 = K + 4 µ/3 and 3K µ/f3.

§1.5.2 Spherical Harmonic Formulation for a Liquid Core Applying the spherical harmonic coefficient notation to the incremental core
equations (1.1.11), (1.1.12) and (1.1.13) yields the following system of differential equations for the spherical harmonic coefficients in the core:

<tfn=qn-(nt I+ 4,rgG(op) (o)) <Pn

l o)) 2
, = I n (n + 1) _ 4nG p(0)] )

l) -( ( n +

4nG p(

qn

,2

g(o)

<Pn

r2 + g(O) qn

(1.5.21) (1.5.22)

Tm =0
Ten= 0

(1.5.23) (1.5.24)

Ir•,

</Jn
Un= g(o)

(1.5 .25)

v = n

r n(n+

1)

u' +
n

n(nn++21)

u n

(1.5.26)

r

Deformation ofa Spherical Body by Surface Loading

"11111111
Page 28

These equations do not apply if the increase in pressure is hydrostatic (i.e. for the case where n = 0 ) where the dilatation and material increment in the stress tensor are both

non-zero. Equations (1.5.21) and (1.5.22) may be solved as a Runge-Kutta system with two
independent solutions. One of these solutions produces a singularity at the origin when substituted into (1.5.25) and so may be neglected. Assuming that for small values of r the solutions take the form of a power series gives the following form for the bounded
solution:

</Jn = CI

qn =Cr(~ -4:irGp(o)) g(o)

(1.5.27)

where C1 is a constant to be determined from the boundary conditions at the surface via the core-mantle boundary. Although the form for gn given in (1.5.27) is not bounded at the origin, the local increment of gravitational acceleration, which from (1.5.10) has the form gn r'1, is well-defined for all values of r for n > I (again the case n = 0 must be
considered separately). The behaviour of the solutions in the mantle and at the surface is insensitive to the
properties of the inner core so that the simplifying assumption that the inner core is comprised of an inviscid fluid will have a negligible effect on the predicted deformation of the mantle and lithosphere.

§1.5.3 Boundary Conditions Our assumption of radial symmetry means that the boundaries between layers with
different rheological properties will be a series of concentric spheres centred for
convenience on the origin. As was discussed in section 1.2 the boundary conditions within an elastic body
consist of the requirement of continuity in the radial stress, gravitational potential, and displacement throughout the body, and the gravitational boundary condition given in equation (1.2.1). The first three conditions translate to continuity of the corresponding spherical harmonic coefficients defined in equations (1.5.2), (1.5.3) and (1.5.6). If we then consider the continuity of the spherical harmonic components of the quantity q
defined in equation (1.5.11) we see that:

[qn(r)]: =[¢n(r)r- 1]:-[gn(r)J:-4nG[p(?)(r)un(r)]:

(1.5.28)

which is equal to zero from equation (1.2.1)

Deformation ofa Spherical Body by Surface Loading

Page 29

At the core-mantle boundary however the situation is slightly more complex. The fact that the core is incapable of supporting shear-stress leads to an apparent discontinuity in displacement and radial stress across the boundary. This apparent inconsistency between the governing equations in the core and mantle can be reconciled by inserting a thin boundary layer just below the core-mantle boundary in which the dilatation is nonzero and both the radial and tangential components of deformation change rapidly. We still have smoothness and continuity of the gravitational potential, however, and the perturbation to gravitational acceleration is given by the discontinuity in density and the radial deformation as in equation (1.5.28) (see for example Dahlen 1974, Dahlen & Fels
1982). Under the conditions given above for the boundary layer, the incremental form of
the field equation (1.1.6) within this region will be dominated by the divergence of the stress tensor and the pre-stress advection term. Neglecting the other terms in the equation, integrating and using the fact that the stress perturbation below the boundary layer is zero gives the boundary condition for radial stress, giving us a complete system
of core-mantle boundary conditions:

Un (b+)

=

<Pn (b-)
g(O)( b)

+

C3

(l.5.29a)

vn(b+) = C2

(l.5.29b)

Tm(b+) =Tm(b-) + p(o)(b-)g(o)(b)[un(b+)- un(b-)] =p(o)(b-)g(o)(b) CJ (l.5.29c)

T0n(b+)=O +) </Jn ( b = </Jn ( b-) qn(b+) = qn(b-)

( 1.5.29d) (l.5.29e) (l.5.29f)

The constants C2 and C3 represent the apparent discontinuity in the radial and tangential components of displacement across the boundary layer. Their values and the value of C1 will be determined by the boundary conditions applied at the surface of the
sphere. We will restrict our attention to the case of a unit load applied to the surface of the
body, larger loads may be modeled by scaling the results of our analysis accordingly. Such a load will exert no horizontal traction at the free surface and we have continuity of radial stress, yielding the surface boundary conditions:

Deformation ofa Spherical Body by Surface Loading

Page 30

T0n(a)=O
Tm (a) = - g(0) (a)

(1.5.30) (1.5.31)

Outside the body p(~) = 0 so that the local increment in potential satisfies Laplace's equation. Substituting into (1.5.9), integrating, applying the condition that the increment in potential must approach zero as radial distance tends to infinity, and using the definition of gn yields the following equations outside the body (i.e. for r > a):

',r1n /=_2nr+l"'"rn

(1.5.32)

gn=n+rl,1'r\n

(1.5.33)

To calculate the local increment in density at the Earth's surface under a load of

magnitude Yn (0, (fJ) , we assume that the added mass of the load lies inside a layer of

negligible thickness just above the earth's surface. We then calculate the change in

density for a small pillbox of material of height Afi and volume ~V = ~hM (where M

is the area of the ends of the pillbox) that includes the added mass but lies mostly below

the earth's surface. The mass of the load material is, to first order, Yn (0, (fJ) M . Letting

the volume of the pillbox become arbitrarily small but keeping the mass added constant

we may substitute back into (1.1.9) to show that the local increment in density is given:

t P~L!)=- y' · (p<O) Un) + Yn (0, tp)

(1.5.34)

Substituting this into the gravitational equation, (1.1.7), allowing the dimensions of the pillbox to approach zero, integrating, and applying the divergence theorem gives:

i ((v t )dv -4nGdA yn (0, q:,) = V · ¢ -4nGp (O)Un

[(v unI = dAn · ¢).-4nGp(O)

(1.5.35)

We can use equations (1.5.11) and (1.5.32) to rewrite (1.5.35) entirely in terms of spherical harmonic coefficients. Doing so yields:

gn(a-) + 4n-G p(o)(a-)un(a)- ; ¢n(a-) = qn(a-) =-4nG

(1.5 .36)

Equations (1.5.30), (1.5.31) and (1.5.36) comprise the surface boundary conditions for a load on the surface of a sphere and may be combined with the boundary

1

ij

Deformation ofa Spherical Body by Surface Loading

...... Page 31

conditions at the core mantle boundary given in (1.5.27) and (1.5.29). There will be three linearly independent vector solutions and the solution we want will be a linear combination of these three. For example let Yi be the solution to the system for the case
when = C1 ~J then our solution vector y may be written y =Ciyi .

§1.5.4 Love Numbers Love numbers are a useful tool in modelling the response of a spherical body to
loading by surface harmonic functions. They are a set of independent, dimensionless, depth-dependent quantities defined in terms of the gravitational potential of the surface
r load. If the gravitational potential of the load may be written as a sum of spherical
harmonic components lf/n (r)( ~ Yn (0, <p) then the Love numbers (see for example Love
1927) satisfy the relations:

t ;)~ u(r, e, ¢) =

W}n(r) Yn( e, cp)f + ln(r) V Yn( e, cp)]

(1.5 .37)

,i.

i

'~ 1

tl(r, 0, ¢) = "t)1 + kn(r))lfln(r)(n Y"(0, cp)

(1.5.38)

where the kn 1/fn Yn term represents the contribution of the perturbation in gravitational
potential due to the load.
(r, The gravitational potential at a point 0, <p) due to a distribution of mass with
spherical harmonic components Yn (0', <p') over the surface of a sphere of radius a is:

.to ~r " r Yn((f, cp')dcp' Jo lfln(r)( Yn( 0, cp) = Ga2 d(/ sin (f Ja2+ r2-2 rcos a 0

(1.5.39)

where a is the angle between the rays from the origin that pass through the points
(a, 0', <p') and (r, 0, <p) . The form of the cosine of a may be deduced from the law of
cosines for a spherical triangle (see for example Spiegel 1968) and is given in equation
(1.5.44) below. Substituting (1.4.6) for the denominator in the integrand in equation (1.5 .39) and
applying the orthogonality relation for spherical harmonic functions yields:

Deformation of a Spherical Body by Surface Loading
~r f i'' ~r 1//n(r)( Y" (0, ¢)=Ga d(f sin (f dip' Y" ((f, t/!').t ( Pn.(cos a)
r i'' = Ga(~)" d(f sin (f dip' Yn((f, t/!')Pn(cos a)

Page 32

(L)n a m) = 24Jnr +Gl a

y (0 n ,"f-'

(1.5.40)

By direct comparison with (1.5.2) and (1.5.3) we see that the Love numbers may be related to the spherical harmonic coefficients of deformation and gravitation by the equation:

hn(r)
ln(r)
l +kn(r)

_ 2n + 1
- 4naG

g(0 )(r)un(r) g(0)(r)vn(r)
¢n(r)

(1.5.41)

§1.5.5 Loads of degree 1 or 0 Surface loads of degree 1 require separate treatment (Farrell 1972) since for such a
load the earth's centre of mass is translated. The combined centre of mass of the earth and load remains stationary however, and surface deformation can only be measured relative to the earth's centre of mass. This is not a problem for loads of other order since the centre of mass of such loads coincides with the centre of mass of the earth.
For the degree 1 problem the differential equations have a non-trivial solution to the homogeneous boundary value problem (the case where at the surface of the sphere
T m= 0 and. P(-6.) = -V P u · ( (0) ) , and at the core-mantle boundary the Ci are all taken to be
zero). This solution, y f, has the form:

T

YH= \ ~ l_ 0 0 g(o)( r) 4n G P(O)(r) 2g(o)(r)

I

7' T' > ,

+ r ' - - - ~ .,,.. ~

(1. 5.42)

Given a solution to a particular non-homogeneous boundary value problem, yf, we may construct a new solution to the boundary value problem, y f = yf + c y f , where c is
an arbitrary constant. We may use this extra degree of freedom to choose c such that the displacement of the centre of mass of the load-earth system is zero. For the degree 1 boundary value problem any solution to the system of differential equations (1.5.1 9) that satisfies any two of the boundary conditions (1. 5.30), (1. 5.31 ), and (1. 5.36') will
. automatically satisfy the third.

Deformation ofa Spherical Body by Surface Loading

Page 33

Given an arbitrary degree 1 surface harmonic defined:

1

Y1(0,

<p)=

~
n --I

C1nYin(e,

'P)

(1.5.43)

there is a point (a, fl, <p*) on the surface of the sphere such that Yi( 0, <p) = 13 cos a
where a here represents the same quantity as defined for equation (1.5.39). The form
for the cosine may be deduced from the law of cosines for spherical triangles:

ql) cos a= cos 0cos fl + sin 0sin fl cos (<p -

(1.5.44)

which in turn allows us to express the coefficients of Y1( 0, cp) in the following form:

c1 _1 = sin <p* sin 0*

= C1 o COS fl

c1 1 = cos cp* sin fl

(1.5.45)

So that the load and its associated deformation are symmetric about the radial position
vector r* of the point (a, fl, cp*) .
The shift in the centre of mass of the earth must balance the shift in the centre of mass of the load which is given by (Farrell 1972):

U com = ? Y ( 0, <P) 1/11 (a) - </Ji (a) 1

(1.5.46)

A
where r* is the unit vector in the direction of r *. To calculate the Love numbers relative to the centre of the earth rather than some
generalised origin we must subtract the shift in the earth's centre of mass from the displacement coefficients. Doing so, and retaining the notation we introduced earlier,
yields the expressions (again due to Farrell 1972):

h 1 = h~ + k~

l1 = Zf + k~

k1 =0

(1.5.47)

which define the displacement due to the degree 1 load. A degree O load corresponds to a force applied uniformly over the surface of the
body. The resulting surface deformation will correspondingly be uniform and entirely due to hydrostatic compression of the material within the body. We are interested primarily in modelling surface loading of the earth by ice sheets and their associated waterloads, which, because of conservation of mass of the ice and water within the system, will never have a component degree O so that a detailed analysis of this case is unnecessary in this ,particular instance.

Deformation ofa Spherical Body by Su,face Loading

Page 34

1.6 Spherical Cap Harmonic Analysis As we discussed in sections 1.4 and 1.5, in a conventional spherical harmonic
analysis we obtain solutions to Laplace's equation by translating into spherical polar coordinates, the solutions are sums of components of the form given in (1.4.2). These components are finite, continuous, and single-valued over the surface of a sphere centred on the origin only if the quantities n and m are both positive integers.
In this section we will examine the theoretical foundations of this spherical harmonic analysis technique with slightly more rigour with a view to adapting the procedure to the high resolution problems we wish to consider. In particular we will employ surface harmonics of non-integral degree to model appropriately well-behaved scalar functions over a spherical cap close to the pole. The price of this technique is loss of orthogonality and increased computational cost per spherical harmonic coefficient. The benefit is a dramatic reduction in the number of coefficients required to achieve a given resolution.
In our analysis of the deformation of a sphere under a load we represent the loading
function, and response as a sum of surface harmonics , Ynm . From the theory of harmonic functions (see for example Gilbarg & Trudinger 1983) we see that we can accurately model any suitably well behaved scalar function by summing such terms to sufficiently high degree. Changing the region of interest from the entire surface of the sphere to a small spherical cap centred at the pole does not affect the validity of our mathematical formulation, simply the details of its implementation.

§1. 6.1 Solving Laplace's Equation In solving Laplace's equation in spherical polar coordinates (equation 1.4.1) we
employ separation of variables, that is we assume that the function in question, f , has the
form f(r, 8, cp)=R(r)e(e)tt(cp) , where R is taken to be a function only of r ,

completely independent of 0 and cp , and (9 and 6 are taken to be dependent only on 0
and cp respectively. Substituting this form for f back into Laplace's equation, (1.4.1 ), multiplying by r2 and then dividing through by f yields:

R" R') G" 0 8) iitt") (r2

+ 2r
R

+

_ 1
sin

e (

sin e+(9' cos

+

_ \
sin

e (

= o

( 1.6.1)

where in this case the prime (') denotes differentiation by the corresponding variable. The first term in the equation is a function only of r but each of the other grouped
terms is independent of r, so that for them to sum to zero they must be constant. That is , there exists a constant Xi such that:

~ .

I11.·.

Deformation ofa Spherical Body by Surface Loading

'11111111
Page 35

t} _1
Slil 0

(B"sin0+e<9'cos0)+

1 (tJ")=
sin2 0

-(r-R" +2rR')=-x

R

1

(l 6 2) • .

Using an analogous line of reasoning we see that there exists a second constant

quantity,

x 2

,

such

that:

(1t}3") =x2

( 1.6.3)

which may be substituted back into (1.6.2) to yield:

~ e" sin +EJ' cos 0 + ( X2 + X ) e =o

sin 0

sin20 1

( 1.6.4)

By substituting 17 = cos 0 into this equation we obtain the governing differential equation
for 17:

)~~)+( + f11 ((1 -112

vr(vr 1)- 1~172 )e=O

( 1.6.5)

where Xi = V1( v1 + 1) and X2 =- di . This is simply Legendre's associated differential
equation which has as its solutions the associated Legendre functions of the first and

second kind of degree

V1

and order

V2 ,

Pv vi11) 1

and

Qv vJ11) 2

respectively (see for

example Abramowitz & Stegun 1972). Associated Legendre functions of the second

kind, Qv vJ 17) , are undefined at the poles and are therefore neglected when modelling 2
physical quantities.

If the

order

V2

is

imaginary

or non-integral

then the

solution,

tJv , 2

of equation

(1.6.3) is either aperiodic or has a period that is not an integral fraction of 2n , so that

tJ v2 ( <p) =t= tJv2 ( <p + 2n), which is clearly undesirable since it implies that tJv2 is multi-

valued on the surface of a sphere. We therefore require that v2 be an integer (due to the

symmetry of the trigonometric functions we can without loss of generality assume v2 1s

positive) so that equation (1.6.3) has solutions of the form:

tJv2 ( <p) = A v2 cos V2 <p + B v2 sin v2 <p

( 1.6.6)

where Av and Bv are arbitrary constants determined by boundary conditions. These

2

2

n] . functions are clearly mutually orthogonal over the interval [0, 2

When the form for Xi used in (1.6.5) is substituted back into (1.6.2) we see that

the equation has solutions:

Deformation ofa Spherical Body by Su,face Loading

Page 36

R (r) = C rv1+ D r-Cv1 + n

VJ

VJ

V1

(1.6.7)

The second term in this expression is undefined at the origin when Re( v1) >-1 while the

first term is undefined for Re (v1) < 0 so that for most of the quantities we model we

require Re( v1) > 0

and Dv = 0. 1

With appropriately defined constants our solutions to

equation (1.6.1) will therefore be sums of components of the form:

!vi u, (r, 0, cp) = Cul r"1( Av, cos V2C)J + Bv,sin V2 cp)Pvlv,(cos e)
+ = CVJ rv1(A V2 Y VJ V2 B V2 Y VJ ) -V2

(1.6.8)

The usefulness of this formulation is largely dependent on the orthogonality of
terms of this type over the region [0, .1r] x [0, 2n] . This orthogonality is a result of the
orthogonality of both the trigonometric functions on [0, 2n] (this establishes mutual
orthogonality for different values of v2 ) and the associated Legendre functions on [O, .1r]
(which establishes mutual orthogonality for different values of v1 ). The orthogonality of trigonometric functions whose periods are integral fractions of
2n is an elementary identity of calculus and forms the basis of Fourier analysis (see for example Kreyszig 1983). The orthogonality of the associated Legendre functions is a similar though rather more complicated issue we will examine in slightly greater detail to utilise some of its more general properties.

§1.6.2 Theory of Sturm-Liouville Boundary Value Problems A complete set of solutions to equation (1.6.5) is generated by solutions to the
associated boundary value problems. Requiring that our solution be continuous and finite
n] v throughout the region [0, constrains 1 to be a positive integer. We also apply two e boundary conditions, one at the pole, = o:

eV1 vi(0) = 0 for V2 =I= 0 aev10
ae =0 8= 0

(1.6 .9)

and at the equator, e= n12 :

Deformation ofa Spherical Body by Su,face Loading

Page 37

v2(~) evl

= 0 for V1 + V2 odd

aeV1 V2

=0

ae e= 1rh

for V 1 + v2 even

Equation (1.6.4) can be recast into the form:

e) - e e e ~ /

1)2

( EJ' sin

X sin 2
C, r, L)

=- 1

(1.6.10) (1.6.11)

restricting our attention to the interval [0, 1r12], this is the traditional form for a Sturm-
Liouville problem with eigenvalue X1 (see for example Boyce & DiPrima 1977) with corresponding separated boundary conditions given by equations (1.6.9) and (1.6.10).
For the purposes of this discussion we may take the order of the the Legendre functions as being fixed since orthogonality between surface harmonic functions of different order is inherited from the trigonometric terms. The order of the functions is
therefore taken here to have a fixed positive integer value, v2 •
We now define two sets of solutions to equation (1.6.11 ), for a given value of V2 , we will denote by 'Ev2 the set of functions for which (1.6.11) holds, with associated
integers v1 such that v1 + v2 is even, the corresponding set of functions for which v1 + v2 is odd we will denote (jv2 • It should be noted that each of these sets has its own
distinct set of boundary conditions that apply to their respective member functions. From the theory of Sturm-Liouville eigenvalue problems (see again Boyce &
* DiPrima 1977), given any two distinct eigenvalues of equation (1.6.11 ), Xi O Xi 1 , then
provided the corresponding eigenfunctions evl OV2 and eV1 1V2 both lie in either 'Ev2 or Yv (i.e. as long as both functions satisfy the same set of boundary conditions) they are
2
orthogonal on the interval [0, 1r12].
Orthogonality between the two sets of eigenfunctions, 'Ev2 and Yv2 , is not guaranteed by Sturm-Liouville since functions from one set do not have common boundary conditions with the members of the other. In this instance orthogonality results
from the symmetry properties of the two sets. The element functions of 'Ev2 are either symmetric or anti-symmetric about 0 = 1r12, depending on the parity of v2 (see for
example Abramowitz & Stegun 1972), if v2 is even then they will be symmetric,
otherwise anti-symmetric. The members of {jv2 have a converse relation, they are anti-
e symmetric about = 1r12 when v2 is even and symmetric otherwise. That is:

Deformation ofa Spherical Body by Surface Loading
ev1v,(;- e) = (-1)"1eVIv,(; + e)

Page 38 (1.6.12)

so that the two sets of functions will have opposite symmetry relations. This in tum
implies that the integral from Oto n of the product of two given eigenfunctions, one

from 'Ev 2

and the second from

Yv2 , will be identically zero.

This is very similar to the

relationship between the trigonometric functions over [0, 2n] though orthogonality in that

case is more easily deduced from direct analytical computation of the integrals.
Orthogonality of the associated Legendre polynomials over the region [o, n] for

fixed values of v2 is therefore a consequence of their symmetry properties, and in tum grants the surface harmonics complete mutual orthogonality over the surface of a sphere,

[0, Jr] X [0, 2n] .
The most significant consequence of mutual orthogonality is that when we attempt

to approximate a given function over the surface of a sphere with a sum of surface harmonics, the coefficient of each harmonic component is independent of the degree of

the approximation and takes a quite convenient form given in equation (1.4.9).

The formalism of Sturm-Liouville boundary problems is worth considering in some

detail in this case since it illustrates very clearly the effects of altering the boundary

conditions. In the case of spherical cap harmonic analysis we do this not by applying

different conditions but by changing the boundaries themselves as will be discussed in

section 1.6.4.

§1.6.3 Resolution ofSurface Features
The resolution of a surface harmonic of degree n, Yn , is given by the distance between consecutive zeroes, so that if we want to model a feature of diameter £ on the surface of the sphere whose radius is a then we will need to sum to degree N given by:

N= na £

(1.6.13)

For this reason 2£ is often called the minimum representable wavelength of a surface
harmonic of degree N . From the form of equation (1.4.7) we see that in order to calculate the coefficients
of a spherical harmonic approximation FN we must perform (N + l) (N + 2) integrals of
the form given in equation (1.4.8). If, for example, we wanted to model a feature of diameter 100km on the surface of the Earth (a = 6371km) we would have to sum to degree 200. If we use standard spherical harmonic procedures to achieve this accuracy we will need to calculate some 40601 coefficients (remembering that we need coefficients

Deformation ofa Spherical Body by Suiface Loading

Page 39

for both the sine and cosine terms), which makes high resolution modelling using this technique quite expensive computationally.

§1.6.4 Spherical Cap Harmonic Analysis The Spherical Cap Harmonic Analysis (SCHA) technique was first developed by
Haines (1985) as a means of modelling variations in a regional geomagnetic field. In this procedure we restrict our attention to some small area, A, on the surface of the sphere, in
particular a cap of angular radius ~ which we may without loss of generality assume to
be centred on the pole e= O .
We intend to use this technique to model deformation of the earth within a spherical cap due to an ice sheet centred around the pole. We will therefore (see for example Johnston 1993) require spherical harmonic expansions for the ice sheet, and the earth's rheology (in terms of Love numbers). These are all well behaved scalar functions and can be approximated to arbitrary accuracy using surface harmonics, but by changing the region of interest we also change the nature of the Sturm-Liouville problem for the associated Legendre functions.
In our cap analysis it is natural to use boundary conditions analogous to those applied in the spherical case:

evl v2(~) = 0 for V1 + V 2 odd

aeVl V2
ae

=0

0= ~

for v1 + V2 even

(1.6 . 14)

In general however, these boundary conditions will require V1 to be non-integral. Since we are modelling functions that are everywhere finite we would normally avoid using Legendre functions of non-integral degree since they become indeterminate for various values of 0. We are only interested in the area inside the cap however, and the quantities we are modelling may be taken to be zero outside A , so that we can use Legendre functions of non-integral order provided their singularities are far enough removed from the centre of the cap (the pole of the sphere).
If we let m be a fixed positive integer value for V2 , then from the theory of SturmLiouville problems we have that the eigenvalues of the corresponding boundary value problem are real and form a countably infinite sequence. If we denote the k -th eigenvalue for which there exists an eigenfunction that satisfies equations (1.6.11 ), (1.6.9) and (1.6.14) by nk then our boundary conditions become:

ir

~

Deformation of a Spherical Body by Su,face Loading

11111111
Page 40

lePn,m(cos e) I•·<= 0 fork+ m even
Pn~m(cos ~)= 0 fork+ m odd

(1.6.15)

where we should note that in general the exact value of nk is dependent on the order m ,
and that in the standard spherical harmonic analysis scheme we have ~ =~ and nk = k .
The advantage of this type of analysis is that for small values of ~ the values of nk are quite widely spread as we see from the approximate relations (De Santis 1985):

(k nk z 27r~ + 0.5)-0.5
2~Jr
nk+ I -nk::::;

(1.6.16)

So that the number of coefficients that need to be calculated in order to model a

surface feature to a given resolution is greatly reduced. For example, as we have already

seen, if we want to model a feature of radius 100 km we would have to expand to at least degree 200. If we were to use spherical cap analysis over a cap of angular radius ~ then

we

would need to

expand our functions

to

degree

nk max

where

kmax

is

the

smallest integer

such that:

kmax > 200.5 ,2r~ - 0.5 z 128 ~ -0.5

(1.6.17)

The latitudinal resolution of a given surface harmonic expansion is determined by the maximum degree of the expansion while the longitudinal resolution is a function of
the maximum order of the expansion. In a conventional spherical harmonic analysis these
two quantities are the same. But more generally, the order of expansion required to resolve a feature of diameter£ (as in equation (1.6.13)) at colatitude 0 is given by:

mmax(e) ~ Jra ~in 0

(1.6.18)

In general though the features being modelled will be closer to the centre of the cap than
the edge, otherwise the assumption that deformation and stress due to the load are zero at the edge of the cap may not apply. This being the case it is normally adequate to set mmax = kmax , the maximum value taken by the index k .
The number of coefficients that need to be calculated to achieve a given resolution is therefore substantially smaller than in the conventional spherical harmonic analysis for the same resolution. Modelling a 100 km diameter feature inside a cap of angular radius 10°

Deformation ofa Spherical Body by Su,face Loading

Page 41

gives kmax = 23 yielding a total of 601 coefficients, less than two percent of that required

to achieve comparable resolution in the conventional analysis. The price of this reduction in the number of coefficients is two-fold. Firstly, the

associated Legendre functions of non-integral degree are simply more expensive to calculate numerically than the associated Legendre polynomials of integral degree,

requiring significant numerical effort even to determine appropriate values of nk . More

importantly, moving to the Legendre functions of non-integral degree results in a loss of

2n]. [o, ;] [o, orthogonality between functions

from

f£v 2

and

Yv 2

over the region

x

We

still have orthogonality within each set, since this is still a Sturm-Liouville problem but

the symmetry property that gave us a completely orthogonal system is no longer

guaranteed, a result of the singularities outside the cap. As a result of this there is a

significant increase in the computational cost of calculating the coefficients of the surface

harmonics since many of the cross-multiplication terms in the least squares problem are

non-zero. In this instance, the difference between our spherical harmonic approximation and

the function being approximated is represented by the integral in equation (1.4.8):

r f' )r d0 sin 0 dq, (FN(0, q,)-F(0, q, = ']JN

(1.6.19)

where the approximating function, FN, is given:

N

n

FJ_ 0, <p) = n~ min F! Ynm( 0, <p)

( 1.6.20)

The coefficients of the approximating series are chosen to minimise the quantity 'DN 1n equation (1.6.19), so that for each value of n and m the coefficients of our spherical harmonic expansion are given by the relation:

l"L a'D: = 2 N d0 sin 0 127t d<p Ynm( 0, <p) Yim( 0, <p)

F nm

1 =0 o

O

r f' -2 d0 sin 0 dq, Y,m(0, q,)F(0, q,) =0

(1.6.21)

In the case considered in equation (1.4.9) only one of the integrals in the summation in this expression is non-zero. Removing orthogonality between surface harmonics introduces the linear system above, making the coefficients dependent on the degree, N, of the approximation, and prone to numerical error due to the instability of the system of linear equations.

Deformation of a Spherical Body by Suiface Loading

Page 42

3
2
1
0
-1
-2
-3
-3 -2 -1 0 1 2 3
Figure 1.1: Results for the SCHA approximation of a square ice sheet of side 400 km and 400 m thick (indicated by the dashed line). Units on axes are degrees , contour spacing is 100 m.
The formalism given earlier for calculating the Love numbers of a spherical elastic body does not depend on them being of integral degree so that it is valid to apply the formalism of the SCHA technique to the problem .of surface deformation. The increased analytical and computational complexity of the procedure raises the possibility of reduced numerical stability however and its suitability to very high resolution problems remains untested.
Figure 1.1 illustrates the results of SCHA analysis used to model a 400 m thick square ice sheet of side 400 km. The cap used had an angular radius of 10° and the harmonic approximation was continued to order 20, maximum degree 185. Attempts to go to a higher order approximation failed as the least squares problem became subject to catastrophic numerical instability.
The resulting computational saving with the change from the conventional to the SCHA technique was undeniably significant, but the fit of the resulting approximation was unsatisfactory, particularly at the centre of the load itself. The instability of the SCHA technique at high degrees severely reduces its usefulness for modelling very high resolution surface features, while the performance of the global spherical harmonic scheme when applied to regional scale rebound problems is more than adequate (Lambeck

Deformation of a Spherical Body by Surface Loading

Page 43

1993a,b). The SCHA procedure therefore offers no advantage over the global scheme for modelling rebound processes.
A second adaptation of the spherical harmonic procedure to regional problems is the Translated Origin Spherical Cap Analysis (or TOSCA) scheme of De Santis (1991 ). Calculation of the Love numbers for an elastic body is however very dependent on the radial symmetry of the earth model, which is destroyed by moving the origin closer to the pole of the sphere, so that this analysis scheme is not valid for rebound problems on a
stratified earth. A third alternative procedure that addresses many of the limitations of the SCHA
technique is the Adjusted Spherical Harmonics Analysis (or ASHA) scheme of De Santis (1992). In this approach, as with the SCHA procedure, we restrict our attention to a cap of angular radius ~ , centred on the pole of the sphere, latitude is then scaled so as to project the cap onto a sphere which may then be treated using orthodox spherical harmonic techniques to achieve an approximate solution.
In his original treatment De Santis maps the cap onto a hemisphere since the boundary conditions for magnetic potential cannot be reduced to a point. While this transformation from cap to hemisphere allows DeSantis to use Legendre functions of integral degree, they must still be treated as two distinct sets of orthogonal functions satisfying different boundary conditions with no orthogonality holding between the sets. Determining the coefficients of the spherical harmonic expansion in this instance therefore still requires that a large least squares problem of the form given in equation (1.6.21) be solved.
In the case of deformation due to a load of small lateral extent we may apply a simple zero boundary condition to the edge of the cap, so that in our application the boundary may be mapped onto the pole and the zero boundary condition applied there. This allows us to map the cap onto a sphere and apply the full mutual orthogonality of the Legendre polynomials of integral degree, making it ideal in principle for use in this class
of application. The ASHA technique relies on the approximation 0 = sin 0 holding throughout the
region of the cap, which is valid for caps of angular radius less than 14° degrees , and
reasonably accurate for radii of up to 22° . Substituting this approximate relation back
into equation (1.6.4) yields:

2
d (9 de2

+l0_

dd (09

+

(

Xe22 +X1 )

B=O

( 1.6 .22 )

f Introducing a new scaled variable for latitude, 0', defined such that, 0' = s0 = 0 , then
yields:

"

•

Deformation of a Spherical Body by Su,face Loading

'11111111
Page 44

d2e + J_ de + ( X2 + X1 ) e = 0
d ff2 0' d 0' 0'2 s2

( 1.6.23)

The solutions to this equation are finite and continuous over the surface of the sphere only
1) when X1 =s2 k(k + and X2 =-m2 for k, m positive integers. An approximate
solution to equation (1.6.23) is then given by the associated Legendre polynomial
Pt(cos ff). De Santis acknowledges that the new latitude 0' is no longer always small
enough for the approximation ff= sin ff to hold so that equation (1.6.23) is not strictly equivalent to equation (1.6.4), but uses the associated Legendre polynomials as approximate solutions to achieve reasonable agreement with the SCHA technique. Mapping from a cap onto a sphere is a powerful analytical tool that would allow us to very conveniently perform calculations over small regions with great accuracy and speed but its validity needs to be closely examined. In the case of modelling the earth's geomagnetic potential the problem is largely a mathematical one of approximating a given function by surface harmonics with only a few constraints placed by the physical processes involved. The problem of calculating the deformation of a spherical body under a surface load is however rather different, the physics of the processes involved must be faithfully reproduced in the mathematical model for the results obtained to have
any significance. Figure 1.2 illustrates the effect of magnifying the region of interest. The mapping
from cap to hemisphere has two immediate and completely artificial consequences, it greatly exaggerates the effect of sphericity and alters the lateral scale of the problem.

a)

b)

Figure 1.2: Illustration of the effect of the ASHA technique, the spherical cap (fig. a) is mapped onto the sphere (fig. b) with consequent lateral expansion of the load, exaggerated sphericity, and anisotropy.

Deformation ofa Spherical Body by Suiface Loading

Page 45

The stresses due to a surface load are significant to depths comparable to the lateral extent of the load, translation from a cap onto a sphere greatly increases the apparent width of the load and therefore the depth to which it will stress the earth, so that the ASHA technique immediately produces unphysical behaviour when applied to the problem of surface deformation. The most significant effect of the mapping however, is that it renders the body in question anisotropic, severely complicating the form of the constitutive equation, and invalidating much of the formalism developed for the Love
number calculations. As with the TOSCA scheme, the mathematical convenience of the ASHA procedure
is undermined by the physical inappropriateness of its formulation . The price of the promised numerical convenience is that the procedure requires fundamental modifications to the mathematical formalism of the problem, modifications that render the system analytically intractable.

Deformation ofa Spherical Body by Surface Loading

Page 46

1.7 Discussion The mathematical formalism for calculating the deformation of an elastic sphere
under a surface load has been the subject of extensive development and now provides an extremely powerful tool for the calculation of glacial rebound and associated sea-level change. The nature of this procedure is such however, that modelling the contribution of high resolution surface features is computationally very expensive, so that incorporation into a direct inversion procedure is not entirely feasible.
Surface harmonics are commonly applied whenever an approximation of a conservative field is required over the surface of a sphere, two particular examples being the modelling of the earth's geomagnetic and gravitational potentials. The failure of the conventional spherical harmonic technique at high resolution has prompted workers in these fields to develop several alternative variations for use in regional modelling such as the SCHA and ASHA techniques examined here.
The SCHA scheme of Haines (1985) does not offer significant improvement over the global spherical harmonic procedure for high resolution or regional scale modelling and suffers from significant numerical instability at very high degrees.
Fitting a series of surface harmonics to observational data is a predominantly mathematical problem, while numerical modelling of rebound processes relies heavily on the physical accuracy of the underlying theory. In the case of both the TOSCA and ASHA procedures of De Santis (1991, 1992), analytical tools that may validly be applied to the former class of problem undermine the theory developed for rebound problems on
a spherical body. High resolution modelling of surface loading therefore remains inconvenient and
expensive in the spherical case, despite the undoubted power and flexibility of the analytical techniques discussed in this chapter. One of the most natural adaptations of the techniques discussed here is to neglect the effect of sphericity entirely and model the earth as a flat, stratified, semi-infinite Maxwell body. For loads of small lateral extent the effect of sphericity will not be significant, this approach has been employed by many workers (see for example McConnell 1965, 1968, Jeffreys 1976, Wolf 1985a,b) and offers significant advantages over the spherical formulation for high resolution modelling.
In the next chapter we will adapt the basic theory of loading problems to the case of a flat earth, and introduce three distinct techniques for modelling the deformation due to a surface load. Using the correspondence principle discussed in section 1.3, a Maxwell body may be modelled mathematically as elastic with frequency-dependent Lame parameters so that we may restrict our analysis in the first instance to the deformation of a stratified, elastic half-space.

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 47

Chapter 2
DEFORMATION OF A STRATIFIED SEMI-INFINITE . ELASTIC BODY
2.0 Introduction If we restrict our attention to an ice sheet whose radius is very small compared with
that of the earth then we may take the earth to be a semi-infinite solid, the upper boundary of which is a flat plane. In such an analysis it is usually more convenient to approximate the complex geometry of an ice sheet by a set of superimposed loads of reasonably simple geometry (such as cylinders or rectangular prisms) , just as in our global analysis we approximate the load by superimposing surface harmonics. It is obviously more useful, though slightly more difficult, to consider rectangular loads since they may more conveniently be used to approximate complex load geometries. Previous workers have favoured the 2-dimensional case of cylindrical loads (Nakiboglu & Lambeck 1982, Wolf 1985a,b), though some attempts have been made to model ice sheets by superposition of prism loads over a square grid (Quinlan & Beaumont 1982).
Our aim ultimately is to model the deformation of a visco-elastic body under a surface load but, as discussed in section 1.3, by invoking the correspondence principle we may use the theory for an elastic body to derive the deformation of a corresponding Maxwell visco-elastic half-space. This chapter is devoted to the development of a number of procedures for modelling the deformation of a stratified, semi-infinite elastic half-
space. In the case of a semi-infinite elastic half-space, the majority of the deformation
caused by a given surface load occurs above a depth approximately equal to the lateral extent of the load. The smaller the load, the smaller the depth to which the load significantly deforms the elastic body. When we superimpose a set of small loads to model a much larger surface load we superimpose their deformation and stress contributions at depths much larger than the lateral extent of the individual loads. While individually these contributions may be negligible they will become significant upon summation and it is therefore important that the deformation and stress at depth due to these individual loads are calculated as accurately as possible.
In this chapter I will examine three approaches to the problem of calculating the deformation of an elastic half-space by a surface load. All three techniques use Fourier transforms to translate the problem into one of solving a set of linear equations over a region of transform space and then inverting to the spatial domain.

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 48

The first and oldest of these techniques is an early approach developed by McConnell (1965) that relies on explicitly solving the equilibrium equations throughout the body, while the second is a slightly more conventional propagator matrix approach discussed by Cathles (1975) that uses McConnell's direct solution technique to apply the boundary condition in the bottom layer of the stratification but propagates the solution to the surface using a slightly more elegant and theoretically convenient technique from the theory of differential equations. The third is a variation of this traditional propagator matrix technique developed by Kennett (1981) originally designed to model elastic wave propagation through a stratified medium.
I will adapt each of these techniques to three dimensional Cartesian coordinates (the most convenient regime for examining the deformation due to a rectangular load), and examine their numerical stability and suitability to the problem we wish to discuss, with particular attention to the accuracy of the calculations at depth. The final section of the chapter will be devoted to a detailed discussion and comparison of the numerical performance of each of the techniques introduced and an analysis of their numerical stability. Our ultimate goal of course is to use the most appropriate numerical technique as part of a formal inverse procedure.
We will start with a reformulation of the problem of the surface deformation of an
elastic body using Cartesian coordinates.

Deformation of a Stratified Semi-Infinite Elastic Body

Page 49

2.1 Formulation of the Problem

t

We will consider a homogeneous, semi-infinite, elastic body with boundary x3 = 0 ,

with x3 increasing in the downward direction (into the solid). We will denote the

component of displacement in the direction of the axis ei by ui and the various

components of material incremental stress by t&8'J . The form of the components of stress

are given by the constitutive equation for an elastic body:

auk auj) (dUi t1(1..8) =A811-~ox+k,, µ

~+~

ox; ox,

}

I

(2 .1.1 )

with repeated indices indicating summation as per the Einstein convention and ~J representing the Kronecker delta function.
As in chapter 1 (equation 1.5.4) we define the dilatation, L1 , to be the divergence of the deformation field:

A
Ll

--

auk
cfx;

--

daXuli

+

au2
Tx;

+

carux3;

(2.1.2)

If we once again assume that the acceleration terms are small enough to be neglected (as discussed in section 1. 1) and neglect the effect of pre-stress and internal buoyancy (to be reintroduced in chapter 4), then the equations of motion may be written (e.g. Love 1927):

(;., + µ) ~~- + µ V2u1 = - pJ;
J

(2.1.3 )

for all values of j, where V2 denotes the Laplacian operator, the Jj are the various components of the body force per unit mass, and /4 and µ are the Lame parameters of the
elastic body.

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 50

2.2 Basic Propagator Matrix Analysis

I

This analysis is based on obtaining direct analytical solutions to the governing system of differential equations throughout the body and applying the appropriate

boundary conditions as discussed by McConnell (1965). We now follow his analysis and

transform the problem into that of solving a set of 3 simultaneous first order differential equations using Fourier transforms . We Fourier transform all of the quantities in

equations (2.1.2) and (2.1.3) with respect to x1 and Xi using transform variables v1 and

v2 respectively.

Let Xi be the Fourier transform of ui, let Ll denote the Fourier transform of Ll ,

and let r ij be the Fourier transform of tijo> . Then make the following definitions:

Z 1 =-iX1 • Z2 =-iX2 • Z3 =-X3 • P=(A+µ)~

'

'

'

(2.2. 1)

T13 = - i r 13 ; T23 = - i r 23 , T33 = - r 33
11·

where i in this case is defined by the relation i2 = -1 . We assume in this case that
there is no body force acting inside the body we are considering (i.e. that the components
of the body force, iJ in (2.1.3), are everywhere zero throughout the body). Then, upon
transforming equations (2.1.2) and (2.1.3), substituting the quantities defined above into the resulting expressions and rearranging we see that within the body:

[a33 -(v? +vi)]z1 = - v}I'

(2.2 .2a)

[a33 -(v? +vi)]z2= - vf
l [a33 - (vf + vf) Z3= µ1 doiP;

(2.2.2b) (2 .2.2c)

z 2z2 a v1 1 + v

az3

P

+ X3 = - A + µ

a2

where

033

represents the second order linear differential operator

~ .
UX3

(2.2. 3)

_-4

.1 1

11

Defonnation of a Stratified Semi-Infinite Elastic Body

Page 51

§2.2.1 Solving the Equilibrium Equations

We want to evaluate the stress and deformation terms throughout the body by

solving this set of simultaneous differential equations and then inverting from the

r

transform domain. We begin by defining the quantity k0 = 1 + ~, rearranging (2.2.3)

and substituting into equations (2.2.2) to yield:

!!: [a33 -((ko + l)v~ +vi)lz,-kov,v2Z2 =kov1

(2.2.4a)

I
I
I

!!: -(vf }v~}l [c)33

+ (ko + 1

k v, ko = Z 2 - 0 V2 Z,

V2

(2.2.4b)

a r [(ko + 1)

33

-

2 ( V1

+ V22)] Z 3

-_

-

ko

Vi

dZ 1
dX3

+ V2 ddxZ32]

(2.2.4c)

1r

where . Rearranging (2.2.4b), substituting into (2.2.4c) and simplifying gives:

I z,] v1

a2 2 ; dx3

- (vf

+ vi)

= 1 {(vf + v}) [vf + (k0 + 1)v:] Z2 ko V2

(2.2.5)

2

4

ax/ ax/ - [(ko

+

2 )

2 VI

+ 2(ko

+

1 )

2] V 2

a

z2

+ ( ko

+ 1) 0 z2

Equations (2.2.4a) and (2.2.4b) may be combined to yield the following equivalence:

2

2

Iaa~ a~ V2

+ + z 21

-

(

2 V1

V22) Z 1]

= V1 [ a

z 22

-

(

2 V1

V22) Z2]

(2.2.6)

which may be incorporated back into (2.2.5). Collecting like terms, simplifying and removing unnecessary constant factors yields the biharmonic equation:

a4 a o + + + z 2
- -4

-

2( vi2

2) a2 z 2

v2

2

( vi2

2)2

-

v2 z2 -

X3

X3

(2.2.7)

We may immediately see that equation (2.2.7) has solutions of the form:

-A

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 52

Z 2

-
-

B eax3

+ Cv2 X3 eax3

+ £

e-ax3

+

Fv 2

X 3

e-ax3

(2.2.Sa)

where B , C , E , and F are arbitrary constants to be determined by the boundary

ljll
i

conditions, and a is defined such that: a2 = v? + v} .

Using the definition of P, substituting (2.2.Sa) back into (2.2.4a) and (2.2.4c),

and then solving the resulting differential equations, yields the additional results:

Z1 =Aeax3 +Cv1X3eax3 +De-ax3 +FV1X3e-ax3

+[(~:~)- z3 =-(v,A 1v2B)eax,

ax3]Ceax3

+ (v,D; V2E) e-ax3 + [(~:~) +ax3] Fe-ax'

(2.2.Sb) (2.2.Sc)

p = - 2ax3 [Ceax3 + Fe-ax3]

(2.2.Sd)

where A and D are also arbitrary constants. Substituting these expressions back into (2.1.1) and using (2.2.1) yields the forms for the components of stress:

(2 2 r )1 T13

µa\I

v2 1

+

v2 2

)

eax3

A

+ V1V2 eax3 B

+

av 1

ax3 -( A+µ µ

eax3 C

(2.2.9a)

\

- (2 2 v2 1

+ v2) e- ax3 D 2

-

V V e- ax3 E 1 2

-

av

I

[

ax 3

+

(

A+µ µ

)]

e-

ax3

F

T23

(v v µa

\ J

v v 1 2

eax3

A

+

2 1

+ 2

2 )

eax3 B

2

+

2

av 2

[

ax 3

-(

A+µ µ

)]

eax3

C

(2.2.9b)

(v v -

V V e- ax3 D 1 2

-

2 1

+ 2

2 )

e-

ax3 E

2

-

2 av 2

[

ax 3

+

(

A+µ

µ

)]

e-

ax3

F

~

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 53

a[(A µ) - ax] T33

=

-

2µ <V1 eax3 A

+ v. eax3 B 2

-

A++2µ

eax3 C
3

(2.2.9c)

lJ )l + V e- ax3 D + V e- ax3 £ +

1

2

ax
3

+

(

Aµ+µ

e- ax3 F

r If we let w be the deformation and stress vector ( Z,, Z2, Z3, T13, T23, T33 and let
y be the vector (A, B, C, D, E, Fr then we have My= w where the entries of M are
given by equations (2.2.8) and (2.2.9) above and more explicitly in section A.1. If we
set v2 = 0 in M the resulting matrix is equivalent to that given for the two dimensional
case by McConnell (1965), as would be expected.

§2.2.2 Application to a Stratified Body We will now assume that our elastic halfspace is stratified into N uniform layers as
illustrated in Figure 2.1, where the Lame parameters for the n- th layer are An and µn ,
and the lower boundary is the plane x3 = hn . We will take the upper boundary of the first layer to be the surface of the body x3 = h0 = 0, while the bottom layer is taken to be semiinfinite so that hN = 00 •

= x3 = h0 0
X3 = h 1 X3 = h2

/4i , µI
~'~

X3 = hN-1

AN, µN

Figure 2.1: Side on view of a layered elastic half-space. The bottom of the n-th layer has
depth x3 = hn, the top of the first layer is the upper boundary of the half-space x3 = 0.

_-A

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 54

l

r Using the analysis outlined in section 2.2.1 we see that at any point inside the n-th
layer (i.e. if hn- I < X3 < hn) w(x3) = Mn(x3)Yn , where Yn = ( An, Bn, c"' D"' En, Fn and

the entries of Mn(x3) are obtained by substituting x3, An, and µn into A.1. It should be
remembered that the entries of the matrices Mn are in fact functions of v1 and v2 and that we have temporarily suppressed this dependence for ease of notation.
As discussed in section 1.3 the boundary conditions between layers are simply continuity of the various components of stress and displacement, this continuity is inherited by the Fourier transforms of these functions as a matter of course.

Mathematically, this may be expressed

Mn(hn)Yn =Mn+ 1(hn)Yn + I

(2.2. l0a)

and may be rewritten:
Ii;
Yn = [(Mn(hn)f Mn+1(hn)lYn+I

(2.2. l0b)

Applying (2.2.l0b) recursively we see that at the surface:

w(o) = Mi(o)[D (Mn(hn}t Mn+1(hn)l YN
= LyN

(2.2.11)

Once we have determined YN we can use (2.2.10) to evaluate the components of displacement and stress at any depth within the solid. The lateral variation of these quantities will be given when we invert from the transform domain.

§2.2.3 Boundary Conditions Physically, we require that the various components of stress and deformation due to
the surface load obey the regularity condition. That is, that they be finite and continuous inside the body and further, that they approach zero as distance from the load becomes arbitrarily large. In the bottom-most layer of the stratification, x3 > hN- i , depth may become arbitrarily large, we therefore require that the magnitude of the various components of stress and displacement decrease with depth within this layer (this need not necessarily be the case near the surface of the body but it is certainly true inside the

~

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 55

bottom-most layer of the stratification). Applied to our formulation this requirement

means that AN, BN, and CN must all be equal to zero.

At the surface the vertical stress, t~f , will be specified by the normal loading

r i

function, and the shearing stresses, t~f and t1f will both be zero. Hence equation
(2.2.11) becomes:

'/

11!

Z1

0

Z2

0

=L Z3

0

0

DN I

0

EN

(2.2.12)

I
i

\ <I>(V1, V2)

FN

II:

where </J is the Fourier transform of the loading function. We can rearrange (2.2.12) by
collecting all of the unknowns into the one vector. This yields:
II:

0 0

-1 0 0 l14 115 116 Z1 0 - 1 0 l24 125 126 Z2

0

- 0 0 - 1 /34 /35 /36 Z3

0

0 0 0 /44 /45 146 DN

0

0 0 0 /54 155 156 EN

</>(vi, v2)

0 0 0 164 165 166 FN

=KU

(2.2.13)

The problem has now been reduced to solving this set of simultaneous linear equations. In order to evaluate the entries of the vector u for particular values of v1 and v2 we need to first evaluate <p and the entries of K at this point, substitute these values into (2.2.13) and then solve. The resulting values of Dn, En, and Fn can then be used along with equations (2.2.10) and (2.2.11) to determine the components of stress and displacement at any depth within the body.
For a uniform load of magnitude p applied over the rectangular area S defined to be
the set of points (x1, X2, x3) such that - a< x1< a ; - b <Xi< b ; x3 = 0 , the Fourier transform, <p , of the loading function may be found in any table of Fourier transforms (see for example Spiegel 1968) and is given by:

_ 4 p sin (av1) sin (bv2)

</>(v1, v2)

V1 V2

(2.2.14)

~

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 56

In practice it will be necessary to solve (2.2.13) over a grid of values of v1 and v2 and then invert the results using a numerical technique (see for example Press et. al. 1986) to give the desired function values over a corresponding grid in the space domain. This being the case it is more appropriate to use a numerically derived value for the Fourier transform of the loading function than the theoretical value given in (2.2.14).
§2.2.4 Discussion The procedure outlined above has several distinct advantages. One of its most
attractive features is that it is relatively uncomplicated theoretically, which makes it convenient to implement numerically. Also, the formulation of the problem is one that makes it relatively easy to calculate the stress and displacement at depth, a desirable property when interpreting the rheological implications of a given loading history.
It also however has several shortcomings. One of the most significant of which is that the entries of the matrices Mn in equation (2.2.1 0a) include both exponentially growing and decaying terms making matrix manipulation numerically unstable at large depths. Also of some concern is the fact that the propagation procedure involves numerically calculating the inverses of a series of matrices (or at least solving the associated linear systems at each layer boundary) which makes the technique computationally expensive and particularly prone to numerical instability given the large variation in the magnitude of the entries at depth. A detailed review of the numerical performance of this procedure and discussion of its limitations will be given in section 2.5.

~

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 57

2.3 Conventional Propagator Matrix Analysis In the previous section we derived an explicit analytical form for the solution of the
governing differential equations within each layer. The coefficients of the solutions were then chosen inside each layer so as to guarantee agreement with the boundary conditions of the system.
In the conventional propagator matrix technique we apply exactly the same boundary condition in the bottom-most layer as we did in the last procedure . Having obtained an explicit form for the solution at this point, we propagate it using a technique from the theory of differential equations (see for example Braun 1983). In practice this technique requires significantly more analytical effort than the basic procedure outlined above but in return it avoids excessive matrix manipulations and provides increased numerical stability. We will discuss the theoretical foundations of this procedure in some detail as it also forms the basis for the wave propagation technique we derive in the next section.
Applying Fourier transforms to the quantities in (2.1.1 ), and using the identities given in (2.2.1) we may immediately see that:

crx;az1 =v1Z3 +Tµ13

carzx2; = V2 Z3 + Tµ23

(2.3.1)

crx; - c1Z3 _ A, +12µ (T33 -A{v1Z1 + v2z2))

Taking the derivative with respect to x3 of equations (2.3.1) and then substituting (2.2.4) into the resulting expressions yields:

µZ{ arl3 = - 1 I 4(A + µ) vT + (A + 2µ) vi) + µv1V2Z2 (3/4 + 2µ) + /4v1T33

I ar23 = - 1 µv1V2Z1(3A + 2µ) + µ2{( A + 2µ) vf + 4(/4 +µ)vi) + Av2T33

aT33 = -(v1T13 + V2T23)
X3

(2.3.2)

..4

1"'

1;

Defonnation ofa Stratified Semi-Infinite Elastic Body

Page 58

It should be noted here that we are once again assuming that there are no body

forces acting within the body. Equations (2.3.1) and (2.3.2) constitute a system of 6

simultaneous first order differential equations. The governing equation of this system

fl

may be written:

11(
I

aw~ ax3 = Aw

(2.3.3)

where the entries of A may be found from (2.3.1) and (2.3.2) and w is correspondingly

defined to be the vector whose entries are the modified transforms of the various

(i.e. components of stress and deformation

w = (z1, Z2, Z3, T13, T23, T33)T) .

I

I

L111

It should be noted that the coefficients in equations (2.3.1) and (2.3.2) (the entries

of the matrix A ) are dependent only on A, µ, v1, and v2 , and not on x3 except insofar as ).,

and µ are depth dependent.

r:

I'
I

(A+ Defining k1 so that k1 = 2µr I ' A may be written explicitly:

0

0

V1

0

0

V2

A = I

-k1AV1

- k1AV2

0

µ( vf vi) 4 k1(A + µ) +

k1q3A + 2µ) V1V2

0

q31 k 1 +2µ)v v1 2 µ(vf + 4 k 1(A + µ) v~) 0

0

0

0

1
µ

0 0

0

1
µ

0

0

0 k1

0

0 k1AV1

0

0 k1AV2

-VI

-V2

0

(2.3.4)

From the general theory of differential equations (again see Braun 1983) we see that

if we have a scalar function w that satisfies the equation:

dawx =f(x)w

(2 .3.5)

for some function f , then w has the form:
w = w ef~f (x) dx' 0

(2.3.6)

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 59

w w(x where 0 = 0). Assuming we are in a region of the halfspace where the Lame
parameters do not vary with depth (i.e. where the matrix A is constant with respect to x3 ) then using (2.3.6) we see that analogously, the solution to (2.3.3) is:

exp(D w(xJ =

A(x/)dx/) = eAfxi -x~Jw0

( 2 -.3. 7)

§2.3.1 Calculating the Exponential ofa Matrix We may extend the definition of the exponential of a scalar to the case of an arbitrary
m x m matrix A using the Taylor series expansion of the exponential function:

=~-, A

oo An

A2

e n=O n. =l+A+2- + ....

(2.3.8)

where I is once again the Kronecker identity matrix. This definition is valid for both
scalar values of A (corresponding to the case where m = I ) and the case where A is a
square matrix. In the latter instance however, equation (2.3.8) is not a particularly convenient expression to evaluate directly for m > 2 . Instead it is usually more appropriate to apply a standard technique from linear algebra (see for example Strang 1980).
Consider an m x m matrix A with m linearly independent eigenvectors,
h 1, b2, .... bm, each with corresponding eigenvalues a 1, ¼, .... am. We then define the
diagonalising matrix of A to be the matrix, D that satisfies the equation:

v- 1AD =A

(2.3.9)

where A= diag (a 1, ¼, .... am) is the matrix whose diagonal entries are the eigenvalues
of A and whose non-diagonal entries are zero. Pre-multiplying both sides of equation (2.3.9) by D it is clear that the columns of
D are the eigenvectors b 1, b2, .... b m • Post-multiplying both sides by v- 1 similarly shows that the rows of v- 1 are the left-eigenvectors of A , yT, yJ, .... y~ (the solutions to
yr yr ). the equation A = ai
Rearranging equation (2.3.9) and substituting back into equation (2.3.8) gives:

....

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 60

A

_ 1

(nAD-

1 )

(vAD-

1 )

(vAD-

1 )

(vAD-

1 )

(vAD-

1 )

e = I + DAD +

I"\

+

,.

+ ....

{I,

_, DA(Dv-')AD- 1 DA(DD- 1)A(nD- 1)AD- 1

r

= I + DAD +

I"\

+

,.

+ ....

=I +DAD-'

+D-A-2n-1

DA3D-I
+ 6

+ ....

= DeAD- 1

(2.3.10)

but eA = diag (ea1, ea2, •••• eam) , so that once D is known we may easily calculate eA.
In the case where A is an m x m matrix but has only k linearly independent eigenvectors, b 1, b2, •.•. bk , with corresponding eigenvalues, a 1, ai,, .... ak , and k < m , we may still evaluate eA using a generalisation of the diagonalisation technique outlined above.
It may be shown (see for example Strang 1980) that for the matrix A there exists a matrix D such that:

A 1 0 0 ··· 0

0 A2 0 ··· 0

v- 1AD

= A

= I

o

o
. . .

A...3

.•.··

o..

0 0 0 ··· A 1

(2.3.11)

where A is a block diagonal matrix whose non-diagonal entries are all zero and whose diagonal entries, Ai , are block matrices of the form:

Ai =I

a- 1 0 ... 0 l

I 0
0

a /- l · · · 0 a-
l

0
. . .

1

0 0 •·· O al-

(2.3.12)

..4

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 61

So that the diagonal entries of Ai are simply ai and the only other non-zero entries are
1's in the superior diagonal immediately adjacent to the main diagonal. A is called the Jordan form of A .
In the case where A has m distinct eigenvectors the A i are simply scalar entries, A is identical to the diagonal matrix in (2.3.9), and D may be constructed from the eigenvectors of A . The current case however is slightly more complicated.
From (2.3.9) we have that AD= DA, so that if A i is an ni x ni matrix then there will be a corresponding set of ni adjacent columns of D , d~, d~, .. .. d~. , such that:
I

Ad~= a1d~
= Ad~ a1d~ + d~ _1 j =2, 3, .... ni

(2.3.13)

It can immediately be seen that d~ =b1 , the eigenvector corresponding to the eigenvalue
d; ai, the other are called the generalised eigenvectors of A . This generalisation applies

also to the left-eigenvectors, although because of the change in the direction of

(a~JT a;+1. multiplication the ordering is reversed, so that

A= aid~i and (a;)T A= aid;+

Once we have constructed A and D we may calculate eA using (2.3 .10) which

still applies, though evaluating eA is not as simple in this case.

From substituting (2.3.11) back into (2.3.8) it is clear that the off-diagonal zero

entries in A suppress any cross multiplication of the A i terms so that:

eA1 0 0 ··· 0

0 eA2 0 ••• 0

I eA =

0

0 eA3 • • • 0
... ... ...

I

0 0 0 .. . eA1

(2 .3 . 14)

We now need to develop a technique for conveniently calculating the exponential of a non-diagonal matrix of the form of A i .
Given an m x m matrix A with m distinct eigenvalues, a 1, ~ am ' •••• (and thus m
linearly independent eigenvectors) there exists a matrix D such that (2.3.9) holds. Given
an arbitrary function f then f(A) may be evaluated using a Taylor series expansion for
J:

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 62

=D(t f(A) = f f,,A" = f fn(DAD-1)"

fnAn)v-1 = Df(A'\n-1

n=O n! n=O n!

n=O n!

JL'

(2.3.15)

where the fn are the Taylor series coefficients for f. Again applying the Taylor series

formulation for f and noting that all of the off-diagonal entries of A are zero, it is clear
that the i- th diagonal element of f(A) is f( ai) and that the off-diagonal elements are
(!( zero i.e. J(A) = diag a1 )J(a,), .... J(am)} .

We now define Mi to be the m x m matrix whose only non-zero entry is M~- = a .

ll

I

and note that:

Mi = ai fI (A - aJ)

]j:=:t;

i i

ai - aJ.

m
A=LMi
i= I

(2.3.16)

Substituting back into (2.3.15) and again using the fact the off-diagonal entries of A and Mi are zero yields:

J(A)

=D~,tM}r

1=D[,t1(M)]ir

1 =

~(vf(M;)v- 1)

fI (A - (A - f(a-)D a.J) =lt1v ai

ajl) -II = I"(~ 1, I

Jj::=t;iI ai - aJ.

rrm

} 1...-,.-}

jj:=:t,ii

a i -

aj

=,tf(a;) D(:-_ajl)

(2.3.17)

J::t;i I aj

which is the standard form for this type of expression (see for example Gantmacher 1960). The case where two of the eigenvalues are equal (i.e. ai = aj for some i and j)
may be considered by letting aj = ai + c in (2.3.17) and letting c approach zero. Doing
so yields:

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 63

))+ (D :--~:)+ (A-azl) f ~(a~a1 f(A)= n=I

f(ai)(r-

I "i: I ,]

I

f'(ai)

k "i: I ,]

I

ktJ(ak) fr k "i: I,]

ll-;kI ak- al

(2.3.18)

which is the same result as that given by Gantmacher (1960). We may obviously use (2.3 .18) to calculate eA directly except that taking the limit
as £ approaches zero becomes prohibitively complicated when the matrix A has more than two identical eigenvalues. In the case of the matrix A given in (2.3.4) however, (2.3.18) may be used in conjunction with the Jordan form approach to give an explicit form for eA . In this case we are interested in calculating the exponential of a 2 x 2 matrix, Q , of the form:

Z) Q = ( O0JZ OJZ

(2.3.19)

which may be substituted back into (2.3.18) to yield:
en = ( ewz zewz) 0 e(J)z

(2.3.20)

§2.3.2 Application to the Problem ofan Elastic Medium The characteristic equation of the matrix A given in (2.3 .4) is:
det (A - a/) =(af - a2)3 =0

(2.3.21)

where a is as defined for equation (2.2.8). Equation (2.3.21) has two solutions, both of order 3, a a 1,2,3 = and a4,s,6 =-a. These eigenvalues each have two corresponding
linearly independent eigenvectors of the form:

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 64

V2

-Vi

I d1,4 = /31

0 µv2a14,

-µv1a14,

0

d2,s = /32

V1
V2
-<Xi,s
2µv1 lli,s I
2µv2£Xis, -2µa2

(2.3.22)

where /31 and /32 are constant scaling factors to guarantee orthonormality of the system.
For· both i = l, 4 the equation (A - a/)x =di has no solution so that the
generalised eigenvectors of A (as defined in equation (2.3.13)) are the solutions to this equation for i = 2, 5 and may be written:

d3,6 = /32

k2V1
lX:3,6
k2V2
lX:3,6 k2 µvl µv2 µ~ ,6

(2.3.23)

where we have defined a new quantity k2 so that
l+3µ
k2 = 2(A+ µ)

(2.3.24)

These vectors di form the columns of the matrix D that gives the Jordan form of the matrix A given in equation (2.3.4). Solving the equations for the left-eigenvectors of A gives corresponding solutions:

µv2a14,

-µv1a14,

I Y1,4 = 'fJ 1,4

0

V2

-Vi
0

= Y3,6 /32

2µvl lX:,,6

2µv2lX:,,6

2µa2 V1

I

V2

ll:3,6

The generalised left-eigenvectors of A take the form:

(2.3.25a)

lF

t•·. '

Deformation ofa Stratified Semi-Infinite Elastic Body

....
Page 65

µvl

µv2

-µai ,s Y2,s = /32 k2V1

(2.3.25b)

¼ ,s

k2V2

¼,s

I

llli '

-k2

= = /3 Setting

/3___
1 1,4

2µ

ci2

a

1,4

,

and

(<2
JJ2

4µ (AA++2µµ)~u2-

in the above expressions gives:

* Yi ·dj = 0 if i j

Y1 ·d1 =y4 ·d4 =2µci2a1,4/31'/J1,4= 1

Y2 ·d2 = Y3 ·d3 = Ys ·ds = Y6 ·d6 = 4µ ( AA++2µµ) d2/322 = l

(2.3.26)

where it is usual to choose /31 = '/3 1 = - '/3 4 • With these relations holding it is clear that
v- 1 is the matrix whose rows are YT When combined these results yield a convenient form for the matrices D and v- 1
such that (2.3.11) holds for the matrix A given in (2.3.4). It remains only to calculate
the matrix eAz where z =x3 - xf and A , the Jordan form of A , can be calculated from
(2.3.11) and (2.3.12):

a OOOO0
0a 10 00
A=I00a000I
0 0 0 -a O O
O O O O -a 1 0 0 0 0 0 -a

(2.3.27)

Using (2.3.15) and (2.3.21) gives:

eaz 0 0 0 0 0

0 eaz zeaz 0 0 0

eAz =I 0
0

0 eaz 0 0 0 0 e-az 0

0
0 I

(2.3.28)

0 0 0 0 e-az ze-az

'
,n,II I

I 0 0 0 0 0 e-az

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 66

Combining equations (2.3.22), (2.3.23), (2.3.25) (2.3.26), and (2.3.28) we have
an explicit expression for eAz = D eAz v- 1 which we may substitute into (2.3.7). The
entries of the matrix eAz are given explicitly in section A.2. If we set v2 equal to zero throughout A.2 and recall that z is negative then our expression for eAz corresponds to
that given by Wolf (1985a). It should again be noted that A , D, A , and v- 1 are not
dependent on x3 except insofar as A and µ are depth dependent.

§2.3.3 Application to a Stratified Body Taking the same stratification as that given in section 2.2.2 we see that if x3 = d for
some depth d such that hn _1 < d < hn then:

dawi; I =Anw(d)
X3 =d

(2.3.29)

where the entries of An are obtained by substituting An and µn into (2.3.4). Using

111:

I I

(2.3.7) we see that (2.3.29) has solution:

~· _;

w(d) = eAn[d -do] w(do)

(2.3.30)

I

for some do such that hn - I <do< hn . In particular if we define Zn = hn- I - hn then:

w(hn - I)= eAn[hn -1 - hn] w(hn) = eAnZn w(hn)

(2.3.31)

11. \
where eAnzn may be calculated by substituting Zn for z in the equations in appendix A. Continuity of the stress-displacement vector across the boundaries between layers
may again be used recursively to yield:

(h ) -[rrN (h ) w(0)

- e e _ A1 z 1 A2z2

···

w eAN-

2 1 N-

I

N-I

-

eAnzn] w
n =I

N-I

(2.3.32)

§2.3.4 Boundary Conditions Once again we apply the condition of continuity of the stress-displacement vector
across the boundaries between layers, apply our load at the surface of the half-space where the horizontal components of stress are taken to be zero, and require that the various components of stress and displacement vanish as distance from the centre of the load increases.

Deformation of a Stratified Semi-Infinite Elastic Body

Page 67

This last boundary condition may be introduced by noting that from section 2.2.2 we have:

w(hN-1) =M~hN-1)YN

(2.3.33)

using the same definitions for Mn and Yn as were given in equation (2.2.1 Oa). Substituting back into (2.3.33) and applying the layer boundary conditions yields:

w(o) = [}] eA,,,]M,,(hN-1)YN

(2.3.34)

which may be re-written:

X1
X2
X3
0 0
¢( V1,V1)

0

0

G

0 DN

EN

FN

(2.3.35)

It I
We may again collect all the unknowns onto the one side of the equation as we did in section 2.2.3 and so reduce the problem to solving the equation:

0 0 0 0 0
¢( V1, V1)

-1 0 Q gl4 g15 g16

0 - 1 0 g24 g25 g26

--

0 0 - 1 g34 g35 g36 0 0 0 g44 g45 g46

0 0 0 g54 g55 g56

0 0 0 g64 g65 g66

Z1
Z2
I Z3
DN EN FN

(2.3.36)

=lu
where </J is the loading function discussed in section 2.2.3, and the entries of J may be calculated from (2.3.34), A.1, and A.2.
As in McConnell's analysis we calculate the matrix J and solve (2.3.36) over a grid of values for v1 and v2 , and then invert from the Fourier domain to get the solution in real space.

~

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 68

§2.3.5 Discussion Like the basic propagator matrix technique the advantages of this procedure are its
relative theoretical simplicity (particularly when applied to the two dimensional problem) and the ease with which it may be adapted to calculate the stress-deformation vector at depth. The slightly more complex theoretical development has however eliminated the need to calculate the inverse matrices used by the basic technique which should make this procedure somewhat more robust numerically.
Directly solving the equations of motion in the bottom-most layer however retains the problem of significant variation in the magnitude of the matrix entries at large depths.
Also of some concern is the fact that the entries of the matrices eAnznconsist largely of
cosh and sinh functions. At large depths a computer with finite accuracy will be unable to distinguish between these functions which will render the system prone to numerical error in the solving routines. These issues will be discussed in greater detail in section 2.5.

I:
I:
I
/11r
Irr:

.Ill'!
' i
I
.......

Deformation of a Stratified Semi-Infinite Elastic Body

Page 69

2.4 Wave Propagation Analysis

This technique was originally developed to model the propagation of waves through

an elastic half-space. It is designed to achieve increased numerical stability by completely

analytically decoupling exponentially increasing and decreasing components of

deformation and stress (the up-going and down-going components of the stress-

displacement field). It does so in such a fashion that only exponentially decreasing

components need be considered during the computational phase, greatly enhancing

stability but at the cost of significantly increased analytical complexity.

I

We start by introducing some new notation, for convenience the more important

1111;
quantities introduced in this section and the last will be presented in the table below.

Much of the new notation is devoted to the projection of standard quantities (the stress-

displacement vector, propagator matrix, etc.) into the eigenspace of A. This change of

basis plays an important part_in the mathematical development of the technique.

Quantity
A D
w(x3)
B(x3)
P(x3, xf)
v(x3) Q(x3, xf)
F(x3, xf)
R,T E'E'
n,m
%,D
g
tp
X
s

Description
Governing matrix of the stress-displacement system Matrix that gives the Jordan form of A
Stress-displacement vector at depth x3
Fundamental matrix solution of 2.3.4
Propagator matrix for w from depth xf to x3
Field vector - projection of w (x3) into eigenspace of A Field propagator - propagates v from depth xf to x3 Transform propagator - propagates from v(xf) to w (x3)
Reflection and Transmission matrices Sub-partitions of Q for a uniform layer Sub-partitions of D1 Up- and down-going components of v Stress-displacement vector for the source of deformation Discontinuity in w due to g Discontinuity in v due to g
Surface vector - projection of P to the surface

Equation 2.3.4 2.3.11 2.3.7 2.4.32 2.4.1 2.4.2
2.4.3 2.4.46
2.4.8 2.4.5 2.4.49 2.4.6 2.4.31 2.4.38 2.4.58 2.4.39

Table 2.1: List of important mathematical quantities in the development of the wave propagator technique

r

Deformation ofa Stratified Semi-Infinite Elastic Body

111111
Page 70

The original formulation of this procedure (Kennett 1981) was intended to model the propagation of waves through an elastic body. In his derivation Kennett keeps the acceleration term in equation (1.1.1), and introduces frequency explicitly into the form for his governing matrix. The development we will give will therefore correspond to the zero frequency case of Kennett. This does not affect the validity of the formalism which is still applicable to the case of static deformations.
We know from the analysis of the previous section that the solution to (2.3.3), the
stress-displacement vector, w(x3), satisfies the relation w(x3) = P(x3, x~) w(x~), where x~) from the general form of (2.3.6) we see that the propagator matrix P(x3, is given:

(E' P(x3, xJ) = exp A(x'Jftx'3)

(2.4.1)

where exp is the exponential function. From this definition it immediately follows that
P(~, x~) = P(xj, x]) P(xL x~).
Now if we choose a matrix D and a corresponding field vector v so that they
D(x satisfy the relation, w(x3) = 3) v(x3), then from (2.3.3) and (2.3.7) we have that:

d(Dv) =A(Dv)
X3

(2.4.2)

and

v(xJ) =[n(xJr 1P(xJ, X3)D(x3)] v(x3) = Q(xJ, X3) v(x3)

(2.4.3)

which may be taken as a definition of the field propagator matrix, Q . From this definition and the properties of the original propagator matrix, P , we have that:

Q(xj, x~) = Q(xj, xn Q(xL x~)

(2.4.4)

Given the form of Q in (2.4.3) it is usually convenient to choose D to be the matrix that gives the Jordan form of A .
It is important to note that this change of basis into the eigenspace of A is purely a mathem~tical formalism and does not correspond to a physical process. Performing this transformation does however make the mathematical representation of some physical effects more intuitive, as we shall see later in this section.

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 71

We follow Kennett by defining two classes of contribution to the stressdisplacement vector and the field vector, the up-going and down-going components. We define an up-going vector quantity to be any vector function whose magnitude increases as depth increases and a down-going vector quantity to be any vector function whose magnitude decreases as depth increases. These definitions coming from the fact that physically, waves decay as they travel further from from their source.
Letting z = x3 - x~ and defining matrices E(z) and E'(z) so that:

eaz O 0 E(z) = 0 eaz zeaz
0 0 eaz

eaz 0 0
E'(z) = I O eaz - zeaz
0 0 eaz

(2.4.5)

then from (2.4.3) and (2.3.28) it is clear that in a region of the halfspace with uniform

elastic properties we may choose D to be the matrix defined in (2.3.22), (2.3.23), and

i

1

(2.3.24) that gives the Jordan form of A. Then substituting the Jordan form from

11I
I

equation (2.3 .11) into the definition of the field propagator matrix yields:

eaz 0 0 0 0 0

0 eaz zeaz 0 0 0

v(x3) =(~) =

0 0

0 eaz 0 0 0 0 e-az 0

0 Iv(x~)
0

0 0 0 0 e-az ze-az

0 0 0 0 0 e-az

E(z) O j v(xf) - , 0 E'(-z)

= Q(x3, x~) v(x~)

(2.4.6)

Thus at any point inside a stratified halfspace the field vector v has both an up-going component, represented by the vector </Ju , and a down-going component, represented by the vector </Jo .

-, \

Deformation of a Stratified Semi-Infinite Elastic Body

Page 72

§2.4.1 Reflection and Transmission ofDeformation Fields The notation and terminology of this procedure were developed for a dynamic case
and it is not immediately obvious that the same concepts extend meaningfully to the static case (though the mathematical formalism is certainly still valid). The concepts of up- and down-going fields and their behaviour might be best illustrated by an example.
Consider a flat semi-infinite, elastic half-space with constant Lame parameters
Ao and ~ throughout. If we place a source of deformation at some depth xf within the
body then there will be a resulting stress-displacement field throughout the body
represented by the vertical displacement u~(x3) which will take its maximum value at ~,
as illustrated schematically in Figure 2.2a. In contrast, consider the case where we have a semi-infinite elastic body with
constant Lame parameters Ao and~ in the region x3 < x~ (with x~ > ~) and constant Lame parameters Ai and µ 1 in the region x3 > x~ (we will assume for this example that
the bottom layer is significantly weaker than the upper, with half the rigidity and bulk modulus). Given the same source of deformation, this second body will have a different
distribution of vertical displacement Uj (x3).
Inserting a layer below the level of the source of displacement will have an effect on the displacement field in the upper layer. In the case of an elastic wave this effect would

i--- -

-

-

-

-

-

A0 , µo

-----Ao,~

.x!; I .,...........,..,,

...............,,...... . .

. .,, . .. .....I X:,

-~

u3 (~ )
Figure 2.2a: Schematic comparison of the difference in displacement fields produced by insertion of a layer boundary . The source of deformation is the same for both bodies, the two displacement curves are compared in the graph on the right-hand side of the figure.
~

Deformation of a Stratified Semi-Infinite Elastic Body

Page 73

0.6 L -· - :_ - - -0 -

a,-.._

0.55 ~

- - 0 - -o _

'-"
·-i:::
.0...

0.5

ac:d 0.45

c.8

(1)
Q
-c:d (.)
·-t::

0.4 0.35

(1)
>

0.3 [

X

0
u 3

- - 0- -
- -0- -

- - 0 - - U1

0.25 L

3

0

10

20

30

40

50

Depth (km)

Figure 2.2b: The displacement curves u~ (4 ) and Uj (.x:,) due to a surface load at the origin (~ = 0 ),
uH plotted against depth, .x:, . 4 ) was calculated for a uniform elastic body of rigidity 3.6x 1010 Pa, and (x uH bulk modulus 6.5x 1010 Pa. Uj 3) was calculated for an half-space identical to that used for 4 )
except for a layer boundary at 50 km below which the rigidity and bulk modulus were both halved.

2

1.98

a,-.._
'-"

1.96

.g~

(1) (.)

1.94

C

(..1.).

·~a- 1.92

1.9

1.88
0

10

20

30

40

50

Depth (km)

Figure 2.2c: Difference between displacement fields , Li u3 (4 )= ul (.x:,)- u~ (4 ), plotted against depth.

~

Deformation of a Stratified Semi-Infinite Elastic Body

--
Page 74

0.25

I

I

I

I

.-_a-_-.,

0. 2 ~ - - - ¾ - - - - 7 E - ~ - - - ¾ - - ~ -

....0ri..:.o.:..:.

0.15 ~

-

§ c.8 0.1 ~

¥-- Component 1 -

-a(1)
..(r.o.).

0.05 ~

- - o - -Component 2

t::

(1)

>

0 r_ - - - - -0 - - - - -G - - - - -G- - - - - - -0- - - - -G- - - - - -

-0.05
0

I

I

10

20

30

40

50

Depth(km)

Figure 2.-2d: Up-going components of the difference in displacement fields ~u3 (.x:,) = u1(.x:,)- u~(.x:,) . The residual down-going component is of the order of 1o-3 m and is not shown in the figure .

l
'Iii

include internal reflections and reverberations between the surface and the layer boundary.

Analogously, we would therefore expect the difference between the two distributions in

the static case to include both an up-going component (corresponding to reflection off the

layer boundary) and a significantly smaller down-going component (corresponding to the

up-going component in tum being reflected off the surface layer).

it

The numerical results obtained for the two displacement curves, u~ and ul , are

I

shown in figure 2.2b. u~ was calculated for a uniform elastic body of bulk modulus

6.5 x 1010 Pa and rigidity 3.6 x 1010 Pa , while ul was calculated for a half-space identical

to the first except for a layer boundary at depth 50 km, below which both bulk modulus

and rigidity were halved. Both deformation fields are due to a square surface load of side

I
Im

100 km corresponding to an ice sheet of uniform thickness, 400 m, centred on the origin.

I

The difference between the two deformation fields , ~ u3 (x3) = ul (x3) - u~(x3) , is

11

given in figure 2.2c, which shows that the contribution due to the inserted layer boundary

I

throughout the uppermost layer is dominated by two separate up-going components,

l,i'

~
I

illustrated explicitly in figure 2.2d. The down-going component of the difference between
the two fields, not shown in the figure, was of the order of 1o-3 m.

1

11'
la:

In the dynamic case reflection off a layer boundary may be defined as the effect of

1:

that boundary on the stress-displacement vector in the region above it. This example

i

i

demonstrates that a completely analogous effect occurs even in the static case. It is

.....olll

11111

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 75

therefore mathematically and conceptually appropriate to retain both the notation and terminology of wave-propagation in the zero frequency case.
We now consider a region of a halfspace, x~ < x3 < xi , sandwiched between two
uniform layers, inside which the elastic parameters of the body are allowed to vary with depth. If we consider a stress-displacement field (represented by the field vector) being
propagated from x~ - to xi+ then from the definition of the field propagator matrix in equation (2.4.3), we have that v(xf-) = Q(x~ -, xi+) v(xi +). This may be re-written:

:o ) v(x~) =(

=[ 11

12

Q Q ] ( :U2 )

DO

Q21 Q22 D2

(2.4.7)

where the Qij are 3 x 3 matrices whose entries are taken from the corresponding parts of
the matrix Q, and the tlJu,Dk are 3 x 1 vectors representing the corresponding up-going
and down-going entries of the field vector v(~). We know from (2.4.6) that v(x~ -) and
v(xi +) may be written in this form inside a uniform region of the halfspace and from
continuity must therefore also have this form at the upper and lower boundaries of the region we are considering. We may now define the transmission and reflection matrices for the up and down-going components of the field vector.
Consider a source of deformation at depth xf < xf that produces only a down-going
component of deformation and assume that there are no other sources of deformation
within the body. Then the down-going component of the field vector at xi+ will be the
portion of the down-going deformation at xf- that has been transmitted from xf- to
xI + . We are not yet considering the effect of any layer boundaries below x3 =xi so that
there will be no up-going component of deformation at this level due to the source at xf ,
and the up-going component of the stress-displacement field at x~ - will be the reflected portion of the initial down-going displacement field. Defining reflection and transmission
matrices RZ2 and I'°o2 for down-going field components such that:

<PD2 = T)}<PDO

RZ <I>uo =

2
<I>vo

(2.4.8)

then substituting into (2.4.7) yields:

l l(; (:uo )= QII Q12

? )

DO

Q2I Q22

D_

(2.4.9)

~

•

--

1111111

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 76

from which it is clear that:

,.,..02
ln

-
-

Q

-212

RD02 = Q12(Q22)- l

(2.4.10)

Similarly, if we consider a source of deformation at depth ~ > xI that produces
only an up-going component of deformation and assume that there are no other sources of deformation within the body, then the up-going component of the field vector at x~ - will
be the portion of the up-going deformation at xI + that has been transmitted from xI + to x~ - . We are not yet considering the effect of any layer boundaries above depth x3 =x~ so that there will be no down-going component of deformation due to the source at xI ,
and the down-going component of the stress-displacement field at ~ + will be the
reflected portion of the initial up-going displacement field. If we define reflection and
transmission matrices Ri2 and T// for up-going field components such that:

t.Puo = T// t.Pu2

t.Pn2 = Ri2<Pu2

(2.4.11)

then substituting back into (2.4.7) yields:

'I

l(=~:) (";r) =[~:: ~::

1:

(2.4.12)

from which it is clear that:

rc/ = QI 1- Q12(Q22r 1Q21

Ru02 = (Q22)- l Q21

(2.4.13)

We may use the forms for the reflection and transmission matrices given in equations (2.4.10) and (2.4.13) to rewrite the field propagator matrix Q . Doing so yields one of the fundamental relations in the formalism of the wave propagation procedure, the partition equation:

Q(xo3 ,

x 2+)

-_ I rz2

-

RZ2(rz2)- IR02 U

Ro2(rr-02)- l
D ln

3

-{~2r1Ri2 (~2rl

(2.4.14)

which relates the reflection and transmission matrices to the field propagator matrix.

11:

~

Deformation of a Stratified Semi-Infinite Elastic Body

Page 77

§2.4.2 Composite Reflection and Transmission matrices

I

Consider the case of two adjacent regions of the halfspace wedged between two uniform layers as illustrated in Figure 2.3 below.

= X3 X~ X3 =xl x3 =Xj

Uniform Layer l Region l Region 2
Uniform Layer 2

Figure 2.3: Side on view of two adjacent regions of the half-space

j

wedged between two uniform layers.

1

l

i

Each region having its own elastic properties which may vary with depth, and

corresponding transmission and reflection matrices. We would like to calculate the overall

transmission and reflection matrices for the two regions considered as a whole.

The chain rule for the field propagator matrix Q is given in equation (2.4.4),

substituting the general form of the partition equation (equation (2.4.14)) back into this

expression yields:

r r ~/ - Rz2(n} IR~2 Rz2(~ 2 I
Q(x~, xi) -
-(~2r1R~2 (~2rl

r?,1 _ ROl(~1)- 1Ro1 ROI(r:_1)- 1 7,1 2 _ R12(7,12)- 1R1 2 R12(7,12)- I

U

D D

U

D D

U

D D

U

D D

-(rolr IRil

X

(~Ir !

-(r}}r 1R!}

(T~2r I

= Q(xt xnQ(xl,xf)

(2.4 .18)

By identifying the corresponding matrix entries we may immediately see that:

....61

Deformation ofa Stratified Semi-Infinite Elastic Body
T~2r (rg2r I =(TZlr 1(1 - Ri1R2)( I
r r r RZ2(ro2 =(rv1 -Rg1(rvr Ri1)R!i(Tb2 +Rg1(11iT1(rb2

Page 78

~ -Rg2 (I1;2f1R?; =(rv1 -Rg1(11iT1Ri1)(rb2 -R2(Tb2tRb2)
- Rz1 (rolr l(T~2r IR~2

(2.4.19)

(ro2r IRi2 (11iT1Ri1(rt2 -R2(Tb2tRb2) + (11iT1(rb2T1R2

from which we may obtain explicit forms for the combined reflection and transmission matrices:

R02 - ROI + r!,l R12 (1 - ROI R12)- IT'!..,l

D- D

U D

U D

D

Rou2

-_

R12 U

+ T12 Ro1(l -R12 Ro1)- 1T12

DU

DU

U

(1 - R T'!..,2 - Tl2 ROI 12)- 1T'!..,l

D- D

U D

D

(2.4.20)

(1 - 1,y-,u()2

--

r,,I U

Rl2 ROI)- ITl2

D U

U

Kennett discusses the physical interpretation and validity of these expressions when applied to elastic waves, they include the effect of all internal reverberations as well as direct transmission and reflection.
It should be noted that if there is a discontinuity in elastic parameters at any of the
regional boundaries, x3 =x~, x3 = xL or x3 =xl then instead of propagating our solution from x~ to xl via xl , we propagate our solution from x~ - to ~ + via xl + . It is not important which side of xi we choose to put the discontinuity as long as we are consistent
throughout. In effect, rather than using the normal propagator matrices, Q(x~, xl} and Q(xL xn 'throughout our analysis we use Q(x~ -, Xj +) and Q(xl +, Xj +) instead.

~

1 -

~

Deformation of a Stratified Semi-Infinite Elastic Body

Iii

11111

Page 79

§2.4.3 Application to a Stratified Body If we consider the case of a stratified body as discussed in section 2.2.2 then we
may use the forms for the composite reflection and transmission matrices given in equation (2.4.20) to construct reflection and transmission matrices for the halfspace as a whole.
In the case where both x3 and x~ lie inside the n-th layer of the halfspace we may choose Dn to be the matrix that yields the Jordan form of An (i.e. the matrix that satisfies
equation (2.3.11)). From section 2.3.3 we see that if An and µn do not vary with depth
then neither do the entries of Dn. Upon substituting Dn for D, equations (2.4.2) and (2.4.3) become:

tx; = [(vnr'Afln] v= [(Dnr'(DnAn(Dnr?} = Anv

(2.4.21)

v(x~) = [(vnr eAn[_,g -x,]vn]V (x3) = [(vnr 1(vn e-A,,z (vnr )vn]v(x3) = e-An'v(x3) (2.4.22)

The form for the eigenvalues in this case is completely independent of depth (from equation (2.3.21)), so the entries of the matrix An are independent of the layer number, n . This is not the case once pre-stress advection or internal buoyancy are included however, as we shall see in chapter 4, so it is useful to retain this notation.
Within each layer of the halfspace the forms for the field propagator matrix given in equations (2.4.6) and (2.4.22) hold. Noting that

v(x~) = Q(x~, x3) v(x3) = Q(xf, x3) Q(x3, x~) v(xf)

(2.4.23)

gives the relation:

Q(x3, x~) =[Q(xi, x3JT1

(2.4.24)

Applying this result to the field propagator matrix of a uniform layer (given in equation (2.4.6)) yields:

Q(xt X3) = I E(zr I 0
0 E'(-zr I

E(-z) 0 0 E'(z)

(2.4.25)

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 80

Which may also be obtained by substituting -z for z in (2.4.6). This result may be
compared with the the partition equation, (2.4.14), to give explicit forms for the reflection and transmission matrices of a uniform layer of thickness Zn , lying between
X3 =hn-l andx3 =hn:

T~ = E'(zn) r; = E(zn)

R~ =0

R; =0

(2.4.26)

using the same definition for Zn as was given in section 2.3.3 (i.e. Zn= hn- I - hn ) . The situation at the boundary between two layers at depth x3 = hn is illustrated in
Figure 2.4 below.

X3 = hn

Layern: ~~ =Anw =[DnA(Dnr]w
f'] Layer n +1: ~~ =An+ ,w = [vn+ ,A(Dn+ 1 w

Figure 2.4: Side-on view of the boundary between two layers and the equations that hold either side of the boundary.

From continuity of the stress-displacement vector across the boundaries between layers we have:

v(hn -) = [(Dn-1r v}(hn +) = Q(hn -, hn +) v(hn +)

(2.4.27)

We may then apply the partition equation, (2.4.14), to this result to obtain the reflection and transmission matrices for the boundary. We may also use equation (2.4.26) to obtain the reflection and transmission matrices for a uniform layer. Once both sets of matrices have been determined they may be combined by applying the form for composite reflection and transmission matrices given in equation (2.4.20). Repeated applications of these results allows us to evaluate the field propagator matrix, Q , throughout the entire half-space, starting from the bottom-most layer and working upwards to the surface.
To illustrate our recursive procedure consider Figure 2.5 below. Starting from hN+ we may apply equation (2.4.27) to calculate the reflection and transmission matrices for
the boundary x3 =hN. Then at any boundary layer hn we assume we have been able to

,rrf.
I'

~·

Defonnation ofa Stratified Semi-Infinite Elastic Body

Iii

111111111

Page 81

calculate the reflection and transmission matrices for the region [hn + 1-, hN+] (using
equation (2.4.27) will yield the reflection and transmission matrices for the bottom-most
layer boundary [hN-, hN+] so that our procedure has a starting point).
Substituting the forms for the transmission and reflection matrices given in equations (2.4.26) and (2.4.27) into the recursion relation (2.4.20) allows us to calculate the reflection and transmission matrices for the region [hn -, hn + 1-] • Applying the recursion relation (2.4.20) again then allows us to combine our results for both regions
and obtain reflection and transmission matrices for the region [hn-, hN+]. We may then
move on to the next layer in the body.

X3 =hn

hn -
hn+

X3 = hn+ 1

= X3 hN-1

Uniform Layer

'Ii I

Figure 2.4 Application of recursive procedure upwards from the bottom layer.

We may derive explicit expressions for the transmission and reflection matrices during each step in our recursive procedure. For ease of notation we make the following definitions:

Ru,D(X3) = Ru,D(x3, hN +) ru,D(x3) = Ru,D(x3 -, X3 +)

Tu,D(x3) = Tu,D(x3, hN +) fu ,D(x3) = Tu,D(x3 -, X3 +)

(2.4.28)

Using these definitions, the recursion relation (2.4.20), and the form for the transmission and reflection matrices of a uniform layer (given in equation (2.4.26)) yields the identities given in equation (2.4.29) below.

1,:

....tllll

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 82

RD(hn +) =T~hn +, hn+l -) RD(hn+l -) TD(hn +, hn+l -)

= E(zn+1) RD(hn+l -) E'(zn+1)
+) TD(hn =TD(hn+ 1-) TD(hn +, hn+l -) =TD(hn+ 1-) E'(zn+ 1)
E( +) T~hn =T~hn+, hn + 1-) Tu(hn+ 1-) = Zn+ 1) Tu(hn+ 1-)
+) Ru(hn = R u(hn+1 -)

(2.4.29)

Now using the recursion relation, (2.4.20), to combine this result with the reflection and transmission matrices of a boundary (given by equation (2.4.27)) yields:
Rihn-)= rD(hn) +tu(hn) RD(hn +) [1-ru(hn)RD(hn +ftD(hn)

Ru(hn -) = Ru(hn +) +TD(hn +) ru(hn)[l -RD(hn +) r)hnfTu(hn +)
+f TAhn -) = TD(hn +) [1-ru(hn)RD(hn tD(hn)

(2.4.30)

Tu(hn -) = tu(hn) [1- RD(hn +)ru(hnfTu(hn +)

It is worth noting again that we are effectively considering the zero frequency case and our waves are not actually travelling through the body but are static deformations. So there is in fact no transmission or reflection due to the internal layer boundaries in the normal sense, but as discussed earlier there is behaviour analogous to transmission and reflection. The transmission and reflection matrix notation also provides a convenient formalism for analysing the zero frequency problem and the similarity in the governing equations of the static and dynamic regimes means that the same relations hold in both settings.
~: l

_....61111

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 83

§2.4.4 Introducing the Source ofDeformation We wish to consider the case of a body that has been deformed from its equilibrium
state by some force applied within it, in particular the case where this force takes the form of a vertical stress applied at the surface of the body. It is standard in this procedure to consider such a surface load as the limiting case of an internal source of deformation restricted to a small region that is allowed to approach the surface. We start therefore, by considering a general stress applied somewhere within the body, and developing a technique for calculating the value of the stress-displacement vector due to this stress at any other point within the body below the level of the source, and then allow first the source of deformation and then the receiver to approach the surface.
We consider then a stratified half-space as discussed in section 2.2.2 with a zero stress condition applying at the surface and some source of deformation, represented by
the stress-displacement vector r, applied within the body. In each layer of the halfspace we define a fundamental solution matrix Bn whose columns are linearly independent
solutions of the governing equation, (2.3.3), (where we have substituted An for A), so
that for some constant coefficient vector cn , we have that w(x3) = Bn(x3) cn • The existence
of such a set of solutions is guaranteed by the theory of differential equations (see for
example Braun 1983). We also define a corresponding field matrix Vn such that at any
depth within the body Bn(x3) = Dn vn(x3) (analogous to the field vector v defined earlier).
r From (2.1.3) we see that the forcing term must be introduced into the field
equations. Restricting our attention to the n -th layer of a stratified halfspace equation (2.3.3) becomes:

ox3w(x3) = An w(x3) + ,{x3)

(2.4.31)

a 4 . where represents the linear differential operator

½

0~

We once again follow the analysis of Kennett (1981) and note from the definitions

of the propagator matrix, P , defined in equation (2.4.1 ), and the fundamental solution

matrix, Bn, defined above, that since the columns of Bn are solutions of (2.3.3),

Bn(x3) =P(x3, x~) Bn(x~), so that where both x3 and x~ lie inside the n -th layer of the
stratification we have:

P(x3, x~) = Bn(x3)(Bn(x~))- I

(2.4.32)

A.

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 84

which may be combined with the governing equation, (2.3.3), to show that:

dx,[P(x3,

xf r']

= dx

,[Bn(xf

)Bn(xJ

1 ]

= Bn(xf

)dx,[ Bn(x3f1]

=

-Bn(xJ)Bn(xJ

1

dx,[Bn(x3)]Bn(x3r

1 _

(2.4.33)

=-[Bn(xf)Bn(x3f1]AnBn(x3)Bn(x3f =-P(x3, xgrAn

where

we

have

used

the

general

identity

O=ox(B B-

1 )

= B(oxB-

1 )

+

(oxB)B-

1

•

Premultiplying (2.4.31) by P(x3,x~r 1 and rearranging yields:

r P(x3, x~r l ,{x3) = P(x~, X3) ,{x3) = P(x3, x~ l [ ax3w(x3)]-P(x3, x~r I Anw(x3)

r = dx,[P(x3, xf 1w(x3)]

(2.4.34)

which may be integrated and then rearranged to yield:

w(x3)

=

P(x3, x~)w(x~)

+

3
J(x~x

P(x3~)

,{

s)

d

(

(2.4.35)

We wish to consider the case where our source of deformation is confined to the
plane x3 = xf , so that r has the form:

,{x3) = Yi ~X3 - xf) + 11 o(x3 -xf)

(2.4.36)

where 8 in this case is the Dirac delta function and 8 its derivative with respect to x3 • We will assume that the source lies within the n-th layer of the stratification (i.e. hn _1 < .xf < hn) . The form of equation (2.4.36) covers a large class of sources of
deformation, and in particular the case of a force applied only across the boundary x3 = 0 . Substituting this term back into equation (2.4.35) gives the integral the form:

J(x"~' P(x3, S) 1{S) d( = H(x3- xff) P(x3, xf)[Yi + riAn] = H(x3 - xn P(x3, xI) P

(2.4.37)

....

Defonnation ofa Stratified Semi-Infinite Elastic Body

Page 85

where H is the Heaviside step function. Equation (2.4.37) therefore defines a discontinuity in the stress-displacement vector w :

w(xf +)-w(xf-) = lJF

(2.4.38)

From equation (2.4.36) we see that for this class of body forces the governing
equation (2.4.31) has the same form as equation (2.3.3) for all x3 =t= .xf. All of the
relations we have derived so far are therefore still valid in these regions, we need simply
to include the effect of the stress-displacement discontinuity lJF at x3 = xf .

§2.4.5 Boundary Conditions for a Buried Source ofDeformation

At the surface, we can see immediately from the definition of the propagator matrix

P that the surface vector, S , due to the stress-displacement discontinuity P at depth

X3

=

s
X3

.
IS

gi.ven

b y:

S = ( ~;) = P(0, xf) 'I'

(2.4.39)

where Sw is the 3 x 1 vector whose entries are the various components of displacement and ST is the 3 x 1 vector whose entries are the various components of stress. That S may be written in this form is a direct consequence of the definition of the stressdisplacement vector w .
The total stress displacement vector at the surface must satisfy the zero stress boundary condition and may therefore be written:

w(o) = (~o)

(2.4.40)

The stress-displacement vector at a depth just below the source of deformation may be related to the stress-displacement field in the bottom-most layer using the propagator matrix P:

w(~ +) = P(~, hN-1)w(hN-1)

(2.4.41)

We may then use (2.4.38) to calculate the stress-displacement vector at a depth just above the source of deformation:

~

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 86

w(zj-) =P(zj, hN-1)w(hN-1) - P

(2.4.42)

which may in tum be propagated to the surface using the propagator matrix P:

w(O)=P(o,x;)w(~-)=P(o, hN_ 1)w(hN_ 1) -P(o,x;) P =P(o, hN-1)w(hN-1) -S

(2.4.43)

We now consider the boundary condition in the bottom-most layer where we require that displacement must be a down-going quantity (i.e. that it must decrease as depth increases). From the form of the field propagator matrix for a uniform layer given in equation (2.4.6), it is clear that within any such layer the solution of the governing equation (2.3.3) may be written as a linear combination of three down-going and three upgoing vector functions, we may therefore take these up-going and down-going vectors to be the basis functions for our solution space. In the bottom-most layer of the stratification we may take the fundamental solution matrix, Bn, to be the matrix whose columns are the up-going and down-going basis vectors of our solution space. Without loss of generality we may take the first three columns of Bn to be the up-going basis vectors and the last three columns to be the down-going basis vectors. Our solution vector will then be a linear combination of the last three columns of Bn and may be written:

w(hN-1)=BN(hN-1)( ~N)

(2.4.44)

Where Cn is a 3 x 1 constant vector whose entries will be determined by the boundary conditions applied at the surface of the body.
We are still free however to specify the initial values for our solution space basis
vectors (i.e. we are free to choose a value for the field matrix vn(x~) at some point, x~ , within the layer). In this case we choose vn(hn _1) =I, which, from the definition of Vn,
yields:

BN(hN-1 +) = DN(hN-1 +)

(2.4.45)

Substituting this result into the form for the solution vector in the bottom-most layer given in equation (2.4.44), and the resulting expression into equation (2.4.43) yields a form for the stress-displacement vector at the surface:

_.....,,.

Deformation ofa Stratified Semi-Infinite Elastic Body
)-s s w(o) =[P(o, hN_ ,) vN]( iN = F(O, hN-1 +)( iN) ~

Page 87 (2.4.46)

which may be taken as a definition of the transform propagator matrix, F .
From the definition of the field propagator matrix Q (equation (2.4.3)) we have:

r Q(a+, hN_I +) =n(a+ n( 1P(o, hN_I) hN_I +)

(2.4.47)

which may be combined with the definition of the transform propagator matrix to give:

F(o, hN-1 +)=n(o+)Q(o+, hN-1 +)=D1Q(o+, hN-1 +)

(2.4.48)

If we partition the matrix D1 into 3 x 3 sub-matrices:
D1 =(mu mD) nu nD

(2.4.49)

and adopt the following notation for the reflection and transmission matrices of the halfspace as a whole:

Ri~D= Ru,D(0 +, hN- 1 +)

ru~D = Tu,D(0 +, hN-1 +)

(2.4.50)

then upon substituting the partition equation and these last two results into equation
(2.4.48) we may re-write the transform propagator in terms of the sub-partitions of D1
and the reflection and transmission matrices of the half-space as a whole:

F12) +) F; 1
F(0, hN - I = ( F2, F22

=(mu mo) T/;N _RiN(~Nr IR~N RiN(~Nr l

nu no

-{~NrIRtN

(~Nrl

(2.4.51)

where the F . are 3 x 3 sub-matrices of F. I)

~

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 88

Now substituting this result, the definition of the surface vector, (2.4.39), and the surface boundary condition, (2.4.40), into the form given in equation (2.4.46) for the stress-displacement vector at the surface yields:

(w0o)

=

(F11
F21

F12) ( 0 ) -(Sw)
F22 CN ST

(2.4.52)

from which we immediately see that: CN = F22-I ST

(2.4.53)

which may in tum be used to show that:

Wo

=

(F12F22-

1 )

ST-

Sw

(2.4.54)

We may then use (2.4.51) to derive explicit formulae for the matrices F12 and F;2 :

r:r F12 = (mD + muR~)( 1

F22 = (nD + nuR~)( T~tr 1

(2.4.55)

which may be substituted into (2.4.54) to yield:

Wo=(mD +muR~)(nD +nuR~rI ST-SW
-RRD = (mD +muRDON) (I ~ ON)- I (nD)- I ST-Sw

(2.4.56)

where R~ = - nD-1 nu . We may use this expression to calculate the surface displacement in the case we are
considering (where we may take Sw = 0 and the only non-zero entry in the vector ST is the vertical component of stress which is given by the loading function </J ). But we will first develop a technique for calculating the stress-displacement vector at any point within the body.

~

r

11

Deformation of a Stratified Semi-Infinite Elastic Body

-

.......

Page 89

§2.4.6 Buried Sources and Receivers
At any depth x3 within the body we may write the field vector v(x3) in terms of its
up-going and down-going components:

v(x3) = I vu(x3)
vD(x3)

(2.4.57)

Now at the level of the source there is a discontinuity in the field vector v , corresponding to the discontinuity in the stress-displacement vector w . From the definition of the field vector we see that this discontinuity, I, is given by:

[v(xl)]: =~xl) =(~) =D(zjr 1 'P

(2.4.58)

then from this result and the definition of the field vector in equation (2.4.3) we see that:

v(xf-) = v(xf +)- ~xf) = Q(xf +, hN-1)v(hN-1 +)-~xi)

(2.4.59)

For notational convenience we define:
Rf;= RD(xf +, hN-1 +)

Ris = RJo, x;-)

(2.4.60)

with similar notation for the other transmission and reflection matrices for these regions. We take the receiver (the point at which we wish to calculate the stress-displacement
vector) to be at depth xf inside the body.
Consider the case where the receiver lies below the the source of deformation (i.e.
where xf < xf < hN- I) then since there are no sources below the receiver the up-going
component of the field vector at this point is simply the total reflected portion of the downgoing component. This may be written:

vu(xf) = RD(xf, hN- 1)vD(xf)

(2 .4.61)

Similarly, if our receiver lies above the source of deformation we have that:
vD(xf) = R~0, x:)vu(xf)

(2.4.62)

__,,_

Deformation ofa Stratified Semi-Infinite Elastic Body

Page 90

From the regularity condition in the bottom-most layer we know that there is no up-
+) going component of the field vector in this region ( i.e. vu(hN- 1 = 0 ) . Using the
notation of (2.4.60) we may apply this last expression, equation (2.4.61), and the
partition equation, (2.4.14) to equation (2.4.59) to yield:

vu(xf -) R~5 vu(xf-)

TisN U

_

RsN D

(risDN)-

1RsN U

RsDN(TsDN)- 1

-(TiNrIRf:

(T;NrI ( vD(h:_ 1 +))-( ~) (2.4.63)

+) We may eliminate vD(hN_1 from this expression to show that:
vu(xf-) = RiN[RZ5vu(~ -) + Eo]- Eu

(2.4.64)

So that from (2.4.62) and (2.4.64) we may derive expressions for both components of the field vector just above the source of deformation:

Vu(X3S

- )=

l[-

R DSR N

uOS]-

l

(

R

SN 0

Lr SD - L~u)

vD(xf-) = Rfsvu(xf-)

(2.4.65)

Using (2.4.58), (2.4.59), (2.4.61), (2.4.65), and the fact that for any two square
matrices A and B for which the relevant inverses exist we have the relations
A(I-Ar 1 =(l-Ar 1 -I and A(I-BAr 1 =(l-ABr 1A we may similarly derive
expressions for the up-going and down-going components of the field vector just below the source of deformation. Doing so yields:

= v0 (x3s +)

[I -

RuOS RDSN]-

l

L( r,5 0

-

RuOS L...D..,S)

v~xf +) = R;t vD(Xj +)

(2.4.66)

Given that we are only interested in displacements due to surface loads we need only consider the case where the receiver lies below the source of deformation .
In this case we may use the definition of the field propagator matrix and its inverse in equations (2.4.3) and (2.4.24) to show that:

~

Deformation ofa Stratified Semi-Infinite Elastic Body
r v(x:) = Q(x:, x; +)v(xf +) = Q(xf +,xf l v(xf +)

From (2.4.14) we then see that:

(Q(xf +, xf}f - (T:rl

- (TuSR)- IRDSR

RtR( TtRr I

ySR - RSR(TSR)- IRSR

D

U U

D

Page 91 (2.4.67)
(2.4.68)

which may be substituted back into (2.4.67) to yield:
vD(x:) = [Ti;R - RtR(TtRr (Ri;R - R:)]vhf +)
r vu(xf) = - (TZR l (R;R - R;N) VD(xf +)

(2.4.69)

Using the recursion relation in equation (2.4.20) we may write RiN in terms of

s TuSR ,

y SDR ,

RSR D ,

RSuR

and

RRN D

.

ubsti.tut1.ng the resu1tm. g expressi.on 1+or the reflect1.on matn.x

into (2.4.69) and applying the matrix identities we used to derive (2.4.66) yields:

R:r vo(x:) = [1 -RtR Ti;RvD(xf +)

(2 .4.70)

We may then extend (2.4.49) to D(xf), and use equation (2.4.70), the definition of
the field vector, v , and the relations between its up- and down-going components given
in equations (2.4.61) and (2.4.66) to derive the displacement vector at xf:

[ l[ ]- i- wR=mDR+mRuRDRN I-RuSR RDRN 1 TDSR [1-RuOS RDSN 1( ~-RuOS Eu)

(2.4.71)

which allows us to calculate the deformation at any point inside the body.

~

r

Deformation ofa Stratified Semi-Infinite Elastic Body

-

....

Page 92

§2.4. 7 Su,face Loads and Receivers If we consider a point x3 that lies above the source then from the general form of the
surface stress-displacement vector given in equation (2.4.46), and the surface boundary condition given in equation (2.4.52) we have that:

(:o) F11(x3) F12(x3) =F(O, x3)v(x3) =

v~x3)

F21(x3) F22(x3) vo(x3)

(2.4.72)

We may now use (2.4.65) and (2.4.72) to show:
v0(x3) = Ru(O, x3)vu(x3) = (F22(x3)tF21(x3)v~x3)

(2.4.73)

Now as x3 approaches zero and we enter the upper-most layer of the stratification, we
D(x have that F(o, x3) = D 3) = D1 , which may be combined with the partitioning of 1
(equation (2.4.49)), the form for the stress-displacement vector at the surface given in
equation (2.4.56), and the above result, to show that:

Rjo, 0+) =(F22(0 +}fF2i(0+) =-no- 1llu =R

(2.4.74)

We may immediately apply this result to find the displacement vector just below the surface. Using (2.4.74) and noting that there are no layer boundaries and negligible attenuation of the stress-displacement vector between the source and the receiver since
RtR riR xf-~ = 0 (so that = 0 and = I ) we have that Ri5 = R. This in tum yields:

Wo+ = [mD +muR~N][1-RR~Nr(~-REu) = [m0 + muR~] [n0 + nuR~r 1(n0 ~ + nuru)

(2.4.75)

Noting that P(o) =(Sw sr)T =D1.E(o) we see that Sr= (no~+nurv) so that
(2.4.75) may be rewritten:

Wo+ = [mo +muR~Nl [no +nuR~rlsT

(2.4.76)

which we could have obtained from (2.4.56) by setting Sw = 0

.....d

Deformation of a Stratified Semi-Infinite Elastic Body

Page 93

This total surface deformation is all we are in fact interested in, the case we are considering does not require us to allow non-zero values of the discontinuity in displacement, Sw , in our form for the source of deformation. We instead specify the vertical stress at the surface through the stress vector ST and calculate the total resulting surface deformation.
Although we could have used (2.4.56) to yield the form for the surface deformation, the formalism of sections 2.4.6 and 2.4.7 is necessary if we want to calculate deformation or stress throughout the body. While we could propagate from a surface vector into the body, doing so without the proper methodology re-introduces exponentially growing terms and undermines the central strength of the procedure. Maintaining complete decoupling of exponentially increasing and decreasing terms is strictly necessary.
As in the previous two procedures we use this technique to obtain the deformation over a grid of values in the transform domain and then numerically transform to get the stress-displacement field in the spatial domain.
' I
,111
§2.4.8 Discussion This procedure is significantly more theoretically complicated than either of the
procedures discussed previously. Despite this it is still reasonably simple to implement numerically though it is less conveniently adapted to determining the stress-displacement vector at depth.
The principle advantage of this technique is that the reflection matrix R~N 1s calculated recursively using (2.4.27), (2.4.29), and (2.4.30) which contain no
exponentially growing terms (E' and E containing positive exponentials of Zn which are defined to be negative), as compared with the propagator matrix techniques developed in sections 2.2 and 2.3. This substantially reduces the variation in magnitude of the matrix entries, making the procedure correspondingly more robust numerically.
Br The inverse matrices we are required to calculate are all of the form [A + I
where B may be a matrix with exponentially decaying entries but A has no depthdependence at all (A is usually either a sub-matrix of some Dn or the identity matrix as in (2.4.75) or (2.4.71)) so that the sum of the matrices approaches A as the entries of B get smaller. This has the effect of keeping variation in the magnitude of the matrix entries to an acceptable level, so that they are not particularly prone to computational inaccuracies during numerical matrix manipulation. A second consequence of this property is that the inverted matrices do not contain any exponentially increasing terms.

~

r

Deformation ofa Stratified Semi-Infinite Elastic Body

-....
Page 94

2.5 Numerical Stability of Propagator Matrix Techniques §2.5.1 Su,face Deformation due to a Square Load: Analytical Techniques
While any two numerical techniques may be consistent, absolute accuracy can only be determined by direct comparison with analytical results. In small-scale engineering problems it is standard (see for example Timoshenko & Goodier 1970, Love 1927) to neglect internal buoyancy and pre-stress advection. This simplification often makes analytical solutions possible and in the particular case of a unit load over a square of side a on the surface of a uniform elastic half-space, the deformation at the centre of the load is given by (Timoshenko & Goodier 1970):

u3(0)= ~ ln(l +12)a (A+lµ) z0.56a (A+lµ)

µ(A+µ)

µ(A+µ

(2.5.1)

where A and µ are the Lame parameters of the body. A similar form may also be
obtained for the average deformation over the area of the load:

ur"" 0.475 a (A+ 2µ)
µ(A+µ)

(2.5.2)

Both forms agree very well with the analysis by Love (1929) in which he derives the form for the deformation due to a rectangular surface load on a flat semi-infinite halfspace. It is worth considering Love's formulation in a little detail as it illustrates some of the difficulties of the engineering approximation.
Given an arbitrary point in the body x = ( x1, ~' x3)T and a point on the surface of
(xi, the body x' = -½, 0 )T , we define the distance function between the two:

1/z

r(x, x') = ( (x1 ·-

xi)2 + (~ -

-½)2 + xi )

(2.5.3)

We then define the logarithmic and Newtonian potentials:
fi r(x) = In (x3 + r(x,x'))dx; dx;
fl V(x)= r(x,x't1 dx;dx;
1,,

(2.5.4)
(2.5.5)

r··,,
I• ...d