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(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 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, n(r) Yn( 0, 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: 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+) = 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 , + 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, -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, . 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µ 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 \ (V1, V2) FN II: where (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

(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 - ~) 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: uo = 2 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