Geophys. J. Int. (1994) 116,321-348 Lamb’s problem of gravitational viscoelasticity: the isochemical, incompressible planet Detlef Wolf Institute of Planetology, University of Miinster, Wilhelm-Klemm-Strape 10, 0-48149 Miinster, Germany Accepted 1993 June 28. Received 1993 June 23; in original form 1992 July 25 SUMMARY We consider a spherical, isochemical, incompressible, non-rotating fluid planet and study infinitesimal, quasi-static, gravitational-viscoelastic perturbations, induced by surface loads, of a hydrostatic initial state. The analytic solution to the incremental field equations and interface conditions governing the problem is derived using a formulation in terms of the isopotential incremental pressure measuring the increment of the hydrostatic initial pressure with respect t o a particular level surface of the gravitational potential. This admits the decoupling of the incremental equilibrium equation from the incremental potential equation. As result, two mutually independent (4 X 4) and (2 X 2) first-order ordinary differential systems in terms of the mechanical and gravitational quantities, respectively, a re obtained, whose integration is algebraically easier than that of the conventional (6 X 6) differential system. In support of various types of application, we provide transfer functions, impulse-response functions and Green’s functions for the full range of incremental field quantities of interest in studies of planetary deformations. The functional forms in the different solution domains involve explicit expressions for the Legendre degrees n = 0, n = 1 and n 2 2, apply to any location in the interior o r exterior of the planet and are valid for any type of generalized Maxwell viscoelasticity and for arbitrary surface loads. Key words: generalized Maxwell viscoelasticity, gravitational viscoelastodynamics, Green’s functions, impulse-response functions, surface loading, transfer functions. 1 INTRODUCTION The problem of the elastostatic deformation of a spherical body was first investigated by Lam6 (1854), who considered a spherical shell subjected to given volume forces and prescribed conditions on its inner and outer surfaces. LamC formulated the field equations in spherical coordinates and derived the solution for the displacement in terms of surface harmonics. Lame’s problem was independently solved by Thomson (1864). In contrast to LamC, he employed Cartesian coordinates and expanded the solution into solid harmonics. Applications of Lam& problem to deformation studies of planetary bodies were initially connected with the problem of correctly accounting for gravitation. The modifications introduced by gravitation were discussed by Love (1908), who pointed out two basic effects. One of them is related to the presence of initial stress in planetary interiors and requires a modification of the ordinary momentum equation valid in the absence of initial stress. Love (1911, pp. 89-93) implemented the necessary adjustments to the theory and derived the incremental momentum equation for a hydrostatic initial state. The other effect only arises if the planet is taken as compressible. In that case, the incremental gravitational force associated with perturbations of the initial density introduces a tendency towards instability. Normally, this tendency is, however, compensated by the opposing force resulting from the compressibility of the material. The stability of planets was studied in detail by Love (1911, pp. 89-104, A simplification of the investigations by Lame, Thomson and Love as far as applications to planets are concerned is the assumption of homogeneous distributions of density, bulk modulus and shear modulus in the initial state. This constraint was removed by Herglotz (1905) and Hoskins (1910, 1920), who gave analytic solutions for elastostatic and gravitational- elastostatic deformations of a sphere, due to tidal volume forces, for simple types of variation of density and elasticity with 1 1 1-125). 321 322 D. Wolf radial distance. Later, Takeuchi (1950) extended these studies in order to allow for more realistic radial variations. This generalization required the use of numerical integration techniques. Takeuchi’s approach was also employed in studies of gravitational-elastostatic perturbations of planetary bodies due to arbitrary surface loads. The latter problem was first investigated by Slichter & Caputo (1960), Caputo (1961, 1962) and Longman (1962, 1963), who calculated the Green’s functions for displacement and incremental gravity. Longman’s theory and numerical results were reviewed and extended by Farrell (1972). More recently, a number of studies were concerned with gravitational-elastostatic perturbations of a compressible planet with fluid core. The solution to this problem proved not to be straightforward and occupied the investigators involved for several years. The competing approaches were reviewed by Longman (1975) and reconciled by Dahlen & Fels (1978); in both works, comprehensive bibliographies of the relevant publications may be found. The development of the theory governing quasi-static, gravitational-viscoelastic perturbations of initially hydrostatic planets has been intimately related to the study of the earth’s response to glacial surface loads. The basic theoretical work was completed by Peltier (1974) and Wu & Peltier (1982); an alternative approach was taken by Sabadini, Yuen & Boschi (1982) and Spada et al. (1992). Essential to both approaches is the formal reduction of the problem to the corresponding elastostatic problem, which is then solved by recasting the incremental field equations into a (6 X 6) first-order ordinary differential system with respect to the radial coordinate, leading to six fundamental solutions. The virtue of the theory presented in these publications is that perturbations of realistic spherically symmetric planets can be calculated. On the other hand, explicit solution functions for elementary models are indispensable for a physical interpretation of the perturbations (e.g. Wolf 1991b) or for tests of the accuracy of numerical solution methods (e.g. Gasperini & Sabadini 1989; Wu 1992). So far, only a limited number of explicit solutions have been obtained on the assumption of incompressible Maxwell viscoelasticity. Thus, Sabadini et al. (1982) stated the six fundamental solutions for a homogeneous spherical layer. Considering the particular case of a homogeneous sphere, Wu & Peltier (1982) derived the special solution for displacement due to surface loading. Dragoni, Yuen & Boschi (1983) gave the special solution for displacement induced by volume forces in a sphere consisting of a homogeneous elastic lithosphere overlying a homogeneous viscoelastic mantle. Wolf (1984) and Amelung & Wolf (1994) studied selected two-layer models and derived special solutions for surface loading. For the same type of forcing, Wu (1990) analysed gravitational-viscoelastic perturbations of a two-layer sphere with arbitrary contrasts of density, shear modulus and viscosity across the interface. In this study, infinitesimal, quasi-static, gravitational-viscoelastic perturbations, due to surface loads, of a spherical, isochemical, incompressible, non-rotating, fluid planet initially in hydrostatic equilibrium are reconsidered. The distinctive features of our analysis are the following. We show that the incremental field equations can be recast into a form in which the equation for the (mechanical) momentum is decoupled from the equation for the (gravitational) potential. The coupling between the mechanical and gravitational aspects of the problem is then restricted to density discontinuities and expressed by incremental interface conditions. Instrumental to the decoupling of the incremental field equations is the use of a field quantity referred to as isopotenrial incremental pressure measuring the increment of the hydrostatic initial pressure with respect to a (perturbed) level surface of the gravitational potential (Section 2). Using the appropriate ansatz for the decoupled incremental field equations, we then establish two mutually independent (4 X 4) and (2 X 2) first-order ordinary differential systems for the mechanical and gravitational aspects of the problem, respectively. The deduction of the general solutions to these systems is algebraically simpler than the deduction of the general solution to the conventional ( 6 x 6 ) differential system; the special solution is obtained in the usual way by means of the incremental interface conditions (Section 3). For the main portion of our study, we are concerned with the derivation of special solution functions. In contrast to previous studies, a comprehensive catalogue of formulae covering all field quantities of interest in studies of planetary deformations is provided. Furthermore, transfer functions, impulse-response functions and Green’s functions for the incremental fields in the appropriate solution domains are collected. The solution functions given involve explicit expressions for the Legendre degrees n = 0, n = 1 and n 2 2 and apply to any location in the interior or exterior of planetary bodies and for arbitrary surface loads. Of theoretical and practical interest is the consideration of generalized Maxwell viscoelasticity, which includes a stability analysis of the solution for this type of viscoelasticity (Section 4). We conclude our study with an assessment of the results obtained and a brief outlook on possible consequences (Section 5 ) . 2 FIELD EQUATIONS A N D INTERFACE CONDITIONS In this section, the basic equations governing infinitesimal, quasi-static, gravitational-viscoelastic perturbations, due to surface loads, of a spherical, isochemical, incompressible, non-rotating, fluid planet initially in hydrostatic equilibrium are compiled. Section 2.1 presents the Cartesian-tensor forms of the field equations and interface conditions for the initial fields (Section 2.1.1) and the incremental fields and their Laplace transforms (Section 2.1.2). Section 2.2 defines the geometry of the problem Lamk’s problem 323 to be solved (Section 2.2.1) and gives the scalar forms in spherical coordinates of the equations for the initial fields (Section 2.2.2) and the Laplace-transformed incremental fields (Section 2.2.3). Key points of this section are (i) the uniform use of the Lagrangian kinematic formulation in the internal and external domains and (ii) the decoupling of the incremental field equations for the mechanical quantities from those for the gravitational quantities. Complete derivations from first principles of the equations in this section can be found elsewhere (Wolf 1993, pp. 14-22, 27-31, 48-50). 2.1 Tensor equations We consider Cartesian-tensor fields in indicia1 notation and imply for them the usual summation and differentiation conventions. Throughout this study, we employ the Lagrangian formulation, Aj... = L j . . . ( X , t ) , where the field value refers to the position, X, + u,(X, t ) , at the current time, t , of a material point whose position, Xi, at the initial time, t = 0, is taken as the spatial argument. The spatial domains of definition for the field equations and interface conditions are 2- U %+ and a%, respectively, with E- the internal domain, %+ the external domain and 8% the interface; 8 = 2- U %+ U 8% is the Euclidian space domain. We also regard the total field value, i, ..., as perturbation of the initial field value, f$!. , such that ij... = f j;!. + f :,?. applies, with fg!. the material incremental value. Accordingly, the temporal argument is the initial time, t = 0, for the initial fields, the current time, r E 5, for the total and the incremental fields and the inverse Laplace time, s E 9, for the Laplace-transformed incremental fields, with 9 the time domain [0, 2) and Spthe (complex) inverse Laplace-time domain. For all Xi E 2- U %+ and t E .T, the field values are taken as continuously differentiable with respect to the arguments as many times as required; jump discontinuities are admitted for Xi E 8%. We begin on the assumption that the fluid completely fills 8 and is isochemical in %- and %+, respectively. More details on the mathematical and notational concepts underlying this study can be found in Appendix A. 2.1.1 For a fluid with the properties specified above, the initial field equations take the forms Equations f o r the initial fields - p y + p4‘y’ = 0, gjo’ = 4(!) . I * 4‘;) = -4aGp, (3) where G is Newton’s gravitational constant, g, the gravity (gravitational force per unit mass), p the (mechanical) pressure, p the volume-mass density and 6, the (gravitational) potential. With p prescribed in 2- and %+, respectively, (1)-(3) constitute a system of partial differential equations for gjo), p(O) and 4(”). The associated initial interface conditions are 2.1.2 Equations f o r the incremental fields On the assumption of infinitesimal, quasi-static, gravitational-viscoelastic perturbations of the specified fluid, the material form of the incremental field equations is given by where tij is the (Cauchy) stress, t‘ the excitation time, ui the material displacement, 6, the Kronecker symbol and p(t - t‘) the shear-relaxation function. With p(t - t’) and p prescribed parameters in 2?- and %+, respectively, and p(O) and +(’) given as special solution to the equations for the initial fields, (7)-( 11) constitute a system of partial differential equations for 324 D. Wolf gi(S),p(*), fl;), ui and 4(6). Supposing a material sheet on the interface, the material form of the associated incremental interface conditions can be written as with nio) the outward unit normal on 82, y = -n!O)gjo) the magnitude of g!') on 8% and u the (incremental) interface-mass density. Note that the definition of y implies that gjo) and nj") are anti-parallel (e.g. Batchelor 1967, pp. 14-20). The derivation of the solution simplifies if the incremental field equations and interface conditions are put into their isopotential-local form. The definitions of material, isopotential and local incremental fields, f j:?., f:!. and fp?., respectively, are considered in Appendix A2. In particular, we use (309) to express p(') and r!;) in terms of p(') and rg) and (310) to express gjs) and I$(') in terms of gjA) and 4('). Taking the Laplace transforms of the resulting equations (e.g. LePage 1980, pp. 285-328) and denoting the Laplace transform of f ( t ) by f(s), we arrive at the following system of partial differential equations and associated interface conditions for &A), $'), Ty), Ei and 6 (I): We note that the incremental equilibrium equation, (18), no longer includes terms accounting for effects due to hydrostatic initial stress and gravitational perturbations and thus formally agrees with the corresponding equation valid in the absence of gravitation. However, such effects now appear in the traction interface condition, (22), which explicitly involves the initial pressure gradient, p+':), and the local incremental potential, &(a). Next, we eliminate ElA) and f f ) from the above equations. For this purpose, we introduce the rotation, Gi, defined by G . = l e . . I 2 q k U k , j , - (25) P"y + 2sp€ijkGk,j =0. with eijk the Levi-Civita symbol. Using (16), (17), (25) and the identity eijkekCm = 6i,Sjm - 6im6jC, eq. (18) takes the form (26) Eqs (16), (20), (25) and (26) constitute an alternative system of partial differential equations in terms of $'), iii, & ( A ) and 8,. The incremental interface condition for follows upon substitution of (17) into (22): [nj0)p(') - spnjo'(ni,j + E ~ , ~ ) + pnjO)($(A) + 4::)~~)lT = -ynjo)el. (27) Eqs (21), (23), (24) and (27) constitute the incremental interface conditions associated with the alternative system of differential equations given above. 2.2 Scalar equations in spherical coordinates 2.2.1 Geometrical considerations We proceed on the assumption that the fluid is initially confined to E. In this case, p(O) = 0, p( t - t ' ) = 0 and p = 0 for X i E %+ and it follows from the incremental field equations and interface conditions that iii = continuous for X i E 8%, Ei = indeterminate for X , E %+ and the remaining mechanical quantities vanish for Xi E 2+. must be parallel planes, co-axial cylinders or concentric spheres (e.g. Batchelor 1967, pp. 14-20). Here, we consider spherically It can be shown that, for a hydrostatic initial state in a non-rotating fluid, the level surfaces of p(O) and LamC's problem 325 symmetric level surfaces and take their common centre as the origin, 0, of a Cartesian coordinate system, OX,X,X3. The spherical coordinates, r , 8 and A , are related to the Cartesian coordinates, X,, X, and X,, by x2 x3 A =tan-' - where r E (0, m) is the radial distance, 8 E (0, z) the colatitude and A E [0,2z) the longitude of the observation point. For brevity, we refer to the fluid initially occupying E- as planet and to the material sheet initially occupying dEas load. With u the radius of the planet, we then have 2- = {Xi I r E (0, a)}, 2+ = {Xi 1 r E (u , y-)} and a%= {Xi I r = a}. Henceforth. we append to symbols of tensors the label subscripts r , 8 and A to denote their appropriate scalar components in spherical coordinates; the summation convention no longer applies. 2.2.2 Equations for the initial$elds On account of the spherical symmetry, the relevant components of the initial field equations and interface conditions, (1)-(6), reduce to -4nr2Gp, r < a r > a ' r = a. [p(O']- = 0 [4"'']' = 0 [4:;']' = 0 These equations are to be supplemented by appropriate conditions ensuring that the initial fields remain bounded as r + 0 and r 4 ~ . 2.2.3 Equations for the incremental fields Since the solution functions in the (r, 8, A , t ) domain are to be expressed in terms of Green's functions representing the contributions from point loads (Section 4.3), it is sufficient to restrict the following analysis to uxisyrnnzetric perturbations. For convenience, we let the X, axis coincide with the symmetry axis. Then, the relevant components of (16), (20), (25) and (26) reduce to I Er, , - (riie),r + 2 6 , = 0 (40) (41) sin 8(rz&$?),r + (sin e i f;)),@ = 0, r # a. Eqs (37)-(40) are four partial differential equations of first order for the mechanical quantities pea), Lf,., Ee and 6,; they are decoupled from (41), which is a second-order partial differential equation in terms of the gravitational quantity & ( A ) . The solutions to these equations must satisfy the appropriate incremental interface conditions. In view of the supposed symmetries, we find for the relevant components of (21), (23), (24) and (27) the expressions (42) I F ' a ) - 2spa,, +,@a' + 4!,"'E,)]- = y b r =a. t (43) (44) 326 D. Wolf These equations are to be supplemented by appropriate conditions ensuring the boundedness of the incremental fields as r + 0 and r + m. Before proceeding to the integration of the incremental equations, we note that a related decoupling method was briefly commented on by Richards & Hager (1984). In contrast to the present study, their scheme is limited to viscous-gravitational perturbations. Hence, the coupled form of the incremental field equations and interface conditions basic to their study involves the local incremental stress, tj?’, and therefore differs from (7)-(15) above, which are written in terms of the material incremental stress, t y ) (cf Wolf 1991a for details on this difference). 3 INTEGRATION OF THE EQUATIONS In this section, the scalar equations in spherical coordinates compiled above are solved. In Section 3.1 the special solution to the initial field equations and interface conditions is given. In Section 3.2, fundamental solutions with respect to 8 satisfying the incremental field equations and interface conditions are sought in terms of Legendre polynomials. Based on this assumption, we deduce the special solution with respect to r in three steps. First, we establish the general solution to the (4 X 4) differential system governing the mechanical quantities (Section 3.2.1). After that, we derive the general solution to the (2 X 2) differential system for the gravitational quantities (Section 3.2.2). Finally, we determine the integration coefficients using the incremental interface conditions (Section 3.2.3). 3.1 Solution for the initial fields The special solution to (31)-(33) which satisfies (34)-(36) and remains bounded as r+ 0 and r+ = is well known (e.g. Ramsey 1981, pp. 45-51). Upon introducing the non-dimensional radial distance, R = r / a , the following formulae are obtained: p c O ) = $ a y p ( l - R2), R < 1, (46) where the additive constant in the potential function has been chosen such that 1imr-= 4 ( O ) = 0. Since y = 4nGap/3, the initial state is completely determined if any two of the parameters a, y and p are given. 3.2 Solution for the incremental fields: ( r , n , s ) domain In the following, we seek solutions to (37)-(41) subject to (42)-(45) in terms of Legendre polynomials of the first kind, P,(cos 8) , where n E {0,1,. . .} is the Legendre degree and P,(cos 8 ) satisfies Legendre’s equation (e.g. Lebedev 1972, pp. 44-5 1). a r , 8, S) = - f ihn(r , s)l,(~)p,,,(cos e), (52) where f , ( s ) is the non-dimensional Legendre coefficient of a(@, s) (Section 3.2.3) and p,(r, s) the normalized Legendre coefficient of f ( r , 8, s). Note that P,(r, s) is assumed to be a Laplace transform, which will be confirmed below (Section 4.3.2). To proceed further, we distinguish the Legendre degrees n = 0 and n 2 1. For brevity, the arguments of the functions will usually be suppressed. Degree n = 0. With Po,, = 0, we may put (53) - I u*o = Qh0 = 0, Lami's problem 327 whence (37)-(40) reduce to ,.jj(a) - 0 ] r < a . 0,r - + 20, = o General solutions are 0 -A(2)R-2 pf) = A(1) ] R < 1 , m- with A(') and A(') arbitrary integration coefficients. Degrees n 2 1. Upon substitution of (49)-(52) and use of Legendre's equation, (37)-(40) take the forms r 1 ' (A) - 3ayj R < 1 Degree n = 1. Using (62), (81), (103) and (104) and substituting the appropriate equations for the eigenvectors and integration coefficients, we obtain I U,, = -3a (126) (127) p"',"'= 0 (128) &$A) = 0, R ZO. (129) tr,, = -3a ] R < 1 , Degrees n 2 2. A s for Legendre degree n = 1, we find YP or,, = -a[n(n + 2)R"-' - (n' - l)R"+'] kns$ + Y P a YP 8,, = - [n(n + 2)R"-' - (n - l)(n + 3 ) R n f 1 ] n knS$ + Y P (2n + 3)(n2 - 1) 2n2 + 4n + 3 kns$ k,sF + y p F z ) = 2a y p R" R<1, With formulae for the four basic incremental field components, ii,., ii,,p"(') and $(A), for n = 0, n = 1 and n 2 2, respectively, we can now give formulae for any other incremental field quantity of interest. In order to relate isopotential increments to local or material increments, we need the solution for the isopotential height, 5, which is given in the following section. 4.1.2 Isopotential height With di the isopotential displacement (Appendix Al) and $(a) = 0 by definition of the isopotential increment (Appendix A 2 ) , (309) and (310) lead to 4(F)d, = -$ (A) . Since n:') is anti-parallel to c$(;), can be replaced by nIo)n,(O)4::)Z, = -$('). An elementary consideration shows that, on the assumption of infinitesimal perturbations, njO)d, equals the height of the current level surface, 4 = 4', passing through X i + di with respect to the associated initial level surface, 4(O) = +', passing through Xi, as measured in the direction of niO). With - = njO)d, the isopotential height and gy' = +(:I, we thus obtain or, in spherical coordinates, Using (121) and K = t7,fnPn, 332 D. Wolf we get - 1 - ff = --@(A) n n , r f a . do) Degree n = 0. Substituting (48) and (125), the preceding equation gives 3aR-', R < 1 "= {3aR, R > 1 ' We note that, with g!"' vanishing for R = 0 and R * 33, Go becomes singular at these points. Degree n = 1. In view of (48) and (129), it follows that R, = 0. Degrees n 2 2. Upon substitution of (48) and (133), we obtain (137) 4.1.3 Strain In Cartesian-tensor notation, the strain, Eij , is defined by e".. = &. . + E. .), i , j j . z Xi E E-. In spherical coordinates, the non-vanishing components of Fij are given by Lamk's problem 333 Degree n = 0. In view of (122) and (123), it follows that I - I E , = E,,, = EFJO = E(2) - E(1) = g.2) = 0 880- A A O A h 0 > R < l . Degree n = 1. Substituting (126) and (127), we arrive at 1 E,, = 0 E,,, = 3R-1 Degrees n 2 2. Using (130) and (131), we obtain the following formulae: Y P Err,, = -(n - l ) [n(n + 2)R"-' - (n + 1)'R"I knsF + Y P Er8,, = (n - l)(n + 2)(Rn-'- R") B$&, = (n + 2)[n2Rn-' - (n' - l ) R n ] YP knsP + Y P Y P knsF + Y P 1 YP n knsF + YP EfJ,, = - - [n(n + 2)Rn-' - (n - l)(n + 3)R"I Ei?,, = -[n(n + 2)R"-' - (n' - l )Rn] I YP knsp + YP 1 Y P n knsF + Y P a Eyin = - [n(n + 2)Rn-2 - (n - l ) ( n + 3)Rn] R < 1 4.1.4 Rotation In spherical coordinates, the non-vanishing component of (25) is Using (118), (119) and 'A = -TiAn!?nPn,8, this takes the form - 1 QAn = (Urn + rUen,r + UTen), r urn Degree n = 0. Upon substitution of (48), (122), (138) and (224), the preceding equations become 0, R < 1, (237) &;) = Degree n = 1. In view of (48), (126), (139) and (226), the following expressions result: 3Y, R < 1 , (239) G:;) = 0, R Z 1 . (240) = Degrees n 2 2. Considering (48), (130), (140) and (227), we obtain ~ 4.1.7 Relations to other studies Provided the underlying field equations and interface conditions are identical, the special solution functions derived above must be identical with those obtained from the conventional (6 X 6) differential system. Special solution functions for models similar to ours have been published by Wu & Peltier (1982), Spada et al. (1992) and Amelung & Wolf (1994). Here, we inspect the formulae in the (r, n, s) domain given by the first authors, which apply to a model identical to ours. Considering the radial displacement as an example, we write (130) in the form Ern = onPn with fin = C,n/[2(2n + 3)]rnt1 + C,r"-'. Upon expressing C, and C, in terms of the model parameters and with the appropriate notational changes applied, some algebra establishes the identity of (130) with the corresponding formulae in Wu & Peltier (1982) [ c t in particular their eqs (5), (29) and (32)]. Of some interest is also the case where the wavelength of planetary deformations is sufficiently short compared to the radius of the planet that the sphericity becomes irrelevant and gravitational perturbations may be neglected. Then, the half-space theory developed elsewhere (e.g. Wolf 1991b, 1993, pp. 33-45) is profitably employed and yields the desired results more easily than the spherical theory considered here. The accuracy of the half-space approximation has been tested computationally for a number of earth models (Wolf 1984; Amelung & Wolf 1994); here, we show how the solution functions for the sphere can be formally reduced to those for the half-space. Since the problem is primarily of theoretical interest, we restrict our analysis to the radial displacement. Using (93), (118), (130), R = 1 + x / a , a, = Ex and the definition of k,, it is expressible as (1 +;)(1+ y-'- (1 -$)(I +;)"" 1 n (2 +:)[ (2 +; + g ) ;sF + Y P Y@,P,. ax = - 1 4 3 n (243) If we put n /a = k, this can be recast into where (1 +i)(l +g - (1 -f)(l + ; k x ) 1 r, = ' (2 +:I[ (2 +; 4 3 + 2 ) k ~ i i + n Supposing now that x /a +. 0 and n + x such that k = finite, some manipulation yields 1 - k x lim r, = n - r 2ksp + y p ' 1 " lim (1 + ; k x ) = ekr. n+s5 (247) According to Hobson (1931, p. 299), En can always be chosen such that, correct to the order l/&, the relation EnPx = s sin [(n + $)@ + n/4] applies in some neighbourhood of 8. Selecting 0 = n/2 and putting ne = ky, this simplifies to @,P, = s cos ( k y ) . (249) In view of (244) and (247)-(249), we obtain for x /a +. 0 and l / n --;r 0 the limit " (1 - kx)ek" cos (ky). (250) - u, = - 2ks& + yp Eqs (249) and (250) are fully consistent with the special solution functions in Wolf (1991b, 1993, pp. 33-45) [ c j in particular eqs (4.40), (4.41) and (4.61) in the second reference] and, accordingly, give the vertical displacement in a half-space subject to a fixed and homogeneous gravity field and a sinusoidal load of wavelength 27r/k and amplitude s. 4.1.8 Transfer functions Inspection of the solution functions listed above shows that the ordinary Legendre coefficients, f n ( r , s) = Fn(r, s ) f , ( s ) , of the field quantities analysed can be decomposed according to Sn(r, s) = E z ( r ) T n ( s ) t n ( s ) , (251) which is the general form of the solution functions in the (r, n, s) domain. Function E ( r ) = F,(r, s ) / T , ( s ) specifies the radial dependence of fn ( r , s) and can be directly obtained from the individual solution functions. Function Tn(s) is referred to as transfer function and found to be of either of two types: c 1, n = l As in the foregoing equations, the arguments of the functions considered will henceforth be displayed for clarity. Lamk’s problem 339 4.2 Generalized Maxwell viscoelasticity For the inversion of ?LA)@) for n 2 2 and T;*)(s) for n 2 1, it is necessary to specify jZ (s). We start from the general formula which expresses p(t - t’) in terms of its spectrum, ~ ( c u ‘ ) . We can approximate the latter to any degree of accuracy required by F(a ’) = x.,“=l p,(9)8(a‘ - a(9)), where a(q) > 0 is the qth inverse elemental Maxwell time and p(q) > 0 the qth elemental elastic shear modulus, both prescribed for q E (1, 2, . . . , Q}; as usual, 8(a’ - a(q) ) denotes the (shifted) Dirac delta function. Using this approximation, (254) reduces to which is the shear-relaxation function for generalized Maxwell viscoelasticity (e.g. Christensen 1982, pp. 16-20; Muller 1986; Wang 1986). Defining p, = limr-r,+o p(t - t’), we obtain in particular pe = xf=l p(9), which is the elastic shear modulus. The Laplace transform of (255) with respect to t - t’ is 4.3 Functions in the (r, n, t ) domain 4.3. I Impulse-response functions We proceed with the transformation of the solution functions specified in Sections 4.1.1-4.1.6 from the (r, n, s) to the (r, n, t ) domain. This requires inverse Laplace transformation of (251)-(253). Details on the inversion of the regular and singular functions entering into these equations can be found elsewhere (e.g. LePage 1980, pp. 285-328; Wolf 1993, pp. 80-82). Inverse Laplace transformation of (251) gives as general form of the solution functions in the ( r , n, t ) domain. Function T,(t - t ’ ) is the impulse-response function associated with the transfer function T,(s), which is of type TiA)(s) or TiB)(s). With p(s) specified and the general functional form of f,(r, t ) established, (252) and (253) can now be inverted. We consider the Legendre degrees n = 0, n = 1 and n 2 2 individually. Degree n = 0. The shifted inverse Laplace transform of (253) is T p ( t - t’) = s(t - t’). Degree n = 1. The shifted inverse Laplace transform of (252) is T$A’(t - t’) = s(t - t’) , Upon substitution of (256) into (253), it follows that whence the shifted inverse Laplace transform takes the form 340 D. Wolf Degrees n 2 2. Substituting (256), we obtain from (252) the equation After some algebraic manipulation, this can be recast into where The shifted inverse Laplace transform of (263) can formally be written as T y ( t - t ’) = y p [ 6 ( t - t’) + Wn(t - t ’ )] . I knpe + YP Comparing (252) and (253), we get T y ) ( s ) = 1 - T P ) ( s ) , whose shifted inverse Laplace transform is T y ( t - t ’ ) = y p [ 525 6 ( t - t ’ ) - Wn(t - t ’ ) ] . k n P e + Y p Y P ~ ~~ 4.3.2 Stability analysis It remains to establish the functional form of Wn(t - t’). Inspecting (264)-(266), we note that Fn(s) can be rewritten as the quotient of two polynomials in s (without common roots) of degrees L = Q - 1 in the enumerator and M = Q in the denominator. Hence, the inverse Laplace transform of m ( s ) exist and, according to the complex inversion formula and the residue theorem, can be specified upon determination of the roots of Vn(s). To prove that all roots are simple and negative, we assume that the M poles of Vn(s) have been ordered such that 0 > -a(’) > -a(2) > . . . > a(M) . Considering the interval $ ( l ) = ( -a(’) , 0) first, we note that lim3+--a(,)+0 V,(s) = --cc and V,(O) = 1, whence one root must lie in 9(’). The remaining roots are found by considering the interval @-) = ( -a(m), -a(m-’)), where rn E {2,3, . . . , M}. Since lims--m(m)+o Vn(s) = --co and l i m s ~ ~ a ~ , , - l ~ ~ o Vn(s) = 2, one root must also lie in each of &*), 9(3), . . . , 9(M). However, Vn(s) can have either M roots if all are simple, or less than M roots if at least one is multiple. Taking this into account, it follows that there is exactly one root in each of 9(’), 9(”), . . . , 9a(M). Having established that vn(s ) has M simple and negative roots, the functional form of Wn(t - t ’ ) can now be given. Denoting the pole in 9‘”) by -Pim), evaluation by means of the complex inversion formula and the residue theorem yields Lamk’s problem 341 The mth term of the sum in (270) is called rnth relaxation mode, with ULm)/VLm) the modal amplitude, pirn) the inverse modal relaxation time and M the total number of relaxation modes. Eqs (270)-(272) completely determine the impulse-response functions TL*)(t - t’) and TLB)(t - t‘) for n 2 2 and, thus, the solution functions in the (r, n, t ) domain. 4.3.3 Maxwell and Burgers viscoelasticity General methods of obtaining closed-form expressions of the roots, -pLm), of V,(s) exist only for M 5 4 . In all other cases, numerical methods must normally be applied. In practice, such methods are even used for M = 3 or M = 4. Here, we evaluate W,(t - 1’) exactly for M = 1 and M = 2. Since M = Q, this is equivalent to Q = 1 and Q = 2, which correspond to Maxwell viscoelasticity and Burgers viscoelasticity, respectively. Case M = 1. Eq. (266) takes the simplified form with the root whereas (270)-(272) lead to The solution for Maxwell viscoelasticity was discussed by Wu & Peltier (1982) and, more recently, by Amelung & Wolf (1993). Its characteristic feature is the single exponential decay mode, with the modal amplitude and modal relaxation time being simple functions of the Legendre degree and the parameters of the planet. Case M = 2. Eq. (266) reduces to which has the following roots: The solution is now characterized by the superposition of two exponential decay modes. Compared to Maxwell viscoelasticity, the complexity of the functional forms for the modal amplitudes and relaxation times is greatly enhanced. A detailed discussion of the response characteristics of elementary models with Burgers viscoelasticity can be found in Riimpker (1990). 4.4 Functions in the (r , 8, A, c) domain The final step is the transformation of the solution function (257) from the (r, n, t ) to the ( r , 8, A , t ) domain. In the following, we will distinguish axisymmetric and non-symmetric loads and calculate the respective Green’s functions. 4.4.1 Axisymmetric Green’s functions We first consider the Green’s function for axisymmetric loads whose distribution is given by c(O’, t ‘ ) , where 8‘ is the colatitude of the excitation point. On the assumption that a(@’, t’) is twice continuously differentiable with respect to 8’ in (0, n) and that 342 D. Wolf J,”[v(8’, t‘)I2 sin 8’ do’ is finite, the distribution can be expanded into a convergent Legendre series (e.g. Lebedev 1972, pp. 53-60): u,(t’) = (n + i) v(8’, t’) sin B’P,(cos 8 ’ ) d8’. c We recall that only the term v,(t‘) corresponds to a net mass; accordingly, no = 0 applies if v(8’, t’) specifies an accreted load, and no = 1 if it specifies a redistributed load. In view of the linearity of the problem, the solution of the load prescribed by (279) can be expressed as which, upon use of (93) and (257), takes the form 1 = P,(cos 8) -cot ep,,,(cos e) 2 F,O [ -P,,,(cos 8 ) ] LT,(t - t‘)a,(t’) dt’. n=no2n + 1 Substituting (280) and changing the sequence of summation and integrations, this becomes where f(””)(r, 8, 8’ , t - t ‘ ) denotes the axisymmetric Green’s function in the (r, 8, t ) domain. 4.4.2 Non-symmetric Green’s functions It is now straightforward to deduce the Green’s function for non-symmetric loads described by the distribution u(8 ’ , A’, t‘), where 8’ and A ’ are the colatitude and longitude of the excitation point, respectively. For this, we take into account that f(””)(r, 8, 8’ , t - t ’ ) is the normalized contribution to f ( r , 8, t ) lrom an annular load at colatitude 8‘. The contribution to f ( r , 8, t ) from a point load on the symmetry axis thus follows from Noting that, for a non-symmetric load, f ( r , 8, A, t ) can be obtained by superposing the contributions from the appropriate distribution of point loads, the generalizations of (283) and (284) are j 7 2 R f f ( r , 8, A , t ) = 1 ~ f “ ’ ) ( r , 6 , t - t ’ ) r ( 8 ’ , A‘, t’) sin 8‘ do’ dA‘ dt’, 1 “ 0 0 -cot 6P,,.(cos 6) f @ ) ( r , 6, t - t ’ ) = - c F,(r) 2ap where f @ ) ( r , 6, t - t’) denotes the non-symmetric Green’s function in the (r, 8, A, t ) domain and 6 the angle between the observation and excitation points, with cos 6 = cos (8 - 8’) cos (A - A‘). Since <(r) = F,,(r, s)/T,(s) is implied by the solution functions listed in Sections 4.1.1-4.1.6, T,(t - t ’ ) of the types T:*)(t - t ’ ) or TLB)(t - t’) given in Section 4.3.1 and u(8 ’ , A’, t’) prescribed, the solution to Lame’s problem of gravitational viscoelasticity is completely specified. Lamt’s problem 343 5 CONCLUDING REMARKS The main results of our study are the following. (1) We have given a complete and rigorous solution describing infinitesimal, quasi-static, gravitational-viscoelastic perturbations, induced by surface loads, of a spherical, isochemical, incompressible, non-rotating, fluid planet initially in hydrostatic equilibrium. The kinematic formulation is uniformly Lagrangian in the internal and external domains of the planet. (2) The solution method adopted differs from the methods conventionally employed for the type of problem studied. Its main characteristic is that the incremental field equations are recast into two mutually decoupled (4 X 4) and (2 X 2) first-order ordinary differential systems in terms of the mechanical and gravitational quantities of the problem, respectively, the coupling being restricted to the incremental interface conditions. A useful property of the decoupled differential systems is that the complexity of the algebraic manipulations necessary to solve them is significantly less than for the (6 X 6) differential system commonly used. It is suggested that this simplification is also of some consequence when considering perturbations of layered planets. (3) The decoupling of the incremental field equations is contingent upon the use of the isopotential incremental pressure measuring the increment of the hydrostatic pressure with respect to a particular isopotential surface. The rigorous definition of isopotential increments and the Lagrangian expressions relating them to the material and local increments conventionally employed in the mechanics of continua subject to initial stress are given in Appendix A. (4) The resulting solution functions are specified for the ( r , n, s), (r , n, t ) , (r , 0, t ) and (r , 8, A , r ) domains. They involve explicit expressions for the Legendre degrees n = 0, n = 1 and n 2 2, are valid at any location in the interior or exterior of the planet and comprise all incremental field quantities of interest in the dynamics of planetary bodies. Of significance is that the solution functions apply to arbitrary types of generalized Maxwell viscoelasticity. The inverse relaxation times characterizing the particular type are given as the poles of the quotient of two polynomials in terms of the Laplace frequency, s. Since all poles are simple and negative, the planet is always stable and its impulse response involves a series of exponential decay modes. Some interesting consequences of our study are the following. (1) The (4 X 4) differential system obtained formally agrees with the system governing the corresponding non-gravitating problem. Available solutions to this simpler problem can therefore be generalized a posteriori in order to include gravitation. Since the numerical modelling of viscoelastic perturbations of planets has so far been based on techniques developed for non-gravitating continua, our solution method opens a way of accounting for initial stress and gravitational perturbations when using these techniques. (2) Allowance for perturbations due to non-gravitational volume forces can be made by an additional term, p(o ) f i s ) , on the left-hand side of (9). Since j$’) = 0, we then have instead of (18) the expression T K ) + p(0)fjA) = 0 and therefore again formal agreement between the isopotential-local form of the incremental field equations and the corresponding equations valid in the absence of initial stress and gravitation. (3) In the case of perturbations of a compressible planet, we have ui,i # 0 and p = p‘” + p‘”. As a consequence, (16) no longer applies, T f j - g~o’(p‘o’iij),j = 0 replaces (18) and p‘” is to be substituted for p in (16)-(24). Hence, no additional mechanical-gravitational coupling is introduced by compressibility, but the formal reduction to the non-gravitating case is not achieved. A special type of compressibility where ( p ‘ ” ) ~ ~ ) , ~ vanishes has recently been studied by Li & Yuen (1987) and Wolf (1991 a). ACKNOWLEDGMENTS Part of this research was carried out while the author was a visiting scientist at the Geophysics Division, Geological Survey of Canada, Ottawa; the generous support provided by A. Lambert and C. Thomson is gratefully acknowledged. Comments by R. Sabadini and two anonymous reviewers on a first version of this paper have been very helpful. REFERENCES Amelung, F. & Wolf, D., 1994. Viscoelastic perturbations of the earth: significance of the incremental gravitational force in models of glacial Batchelor, G. K., 1967. An fnrroducrion ro Fluid Dynamics, Cambridge University Press, Cambridge. Caputo, M., 1961. Deformation of a layered earth by an axially symmetric surface mass distribution, J . geophys. Res., 66, 1479-1483. Caputo, M., 1962. Tables for the deformation of an earth model by surface mass distributions, J. geophys. Rex, 67, 1611-1616. Christensen, R. M., 1982. Theory of Viscoelasticity, 2nd edn, Academic Press, New York. isostasy, Geophys. J . Int., in press. 344 D. Wolf Dahlen, F. A., 1974. On the static deformation of an earth model with a fluid core, Geophys. J . R. astr. SOC., 36, 461-485. Dahlen, F. A. & Fels, S. B., 1978. A physical explanation of the static core paradox, Geophys. J. R. astr. Soc., 55, 317-331. Dragoni, M., Yuen, D. A. & Boschi, E., 1983. Global postseismic deformation in a stratified viscoelastic earth: effects on Chandler wobble Farrell, W. E., 1972. Deformation of the earth by surface loads, Rev. Geophys. Space Phys., 10,761-797. Gasperini, P. & Sabadini, R., 1989. Lateral heterogeneities in mantle viscosity and post-glacial rebound, Geophys. J . Int., 98,413-428. Grafarend, E. W., 1982. Six lectures on geodesy and global geodynamics, Mitt. Geodat. Znst. Tech. Uniu. Graz, 41, 531-685. Herglotz, G., 1905. Uber die Elastizitat der Erde bei Berucksichtigung ihrer variablen Dichte, Z. Marh. Phys., 52, 2i5-299. Hobson, E. W., 1931. The Theory of Spherical and Ellipsoidal Harmonics, Cambridge University Press, Cambridge. Hoskins, L. M., 1910. The strain of a non-gravitating sphere of variable density, Trans. Am. math. SOC., 11,494-504. Hoskins, L. M., 1920. The strain of a gravitating sphere of variable density and elasticity, Trans. Am. math. SOC., 21, 1-43. Lamt, G., 1854. Memoire sur I’tquilibre d’elasticitt des enveloppes sphkriques, Liouville’s J . Marh. Pures Appl., 19, 51-87. Lebedev, N. N., 1972. Special Functions and their Applications, Dover, New York. LePage, W. R., 1980. Complex Variables and the Laplace Transform for Engineers, Dover, New York. Li, G. & Yuen, D. A,, 1987. Viscous relaxation of a compressible spherical shell, Geophys. Res. Lett., 14, 1227-1230. Longman, I. M., 1962. A Green’s function for determining the deformation of the earth under surface mass loads I: theory, J . geophys. Res., Longman, I. M., 1963. A Green’s function for determining the deformation of the earth under surface mass loads 11: computations and Longman, I. M., 1975. On the stability of a spherical gravitating compressible liquid planet without spin, Geophys. J . R. astr. SOC., 42, 621-635. Love, A. E. H., 1908. The gravitational stability of the earth, Phil. Trans. R. SOC. Lond. Ser. A, 207, 171-241. Love, A. E. H., 1911. Some Problems of Geodynamics, Cambridge University Press, Cambridge. Miiller, G., 1986. Generalized Maxwell bodies and estimates of mantle viscosity, Geophys. J . R. asrr. Soc., 87, 1113-1141. Peltier, W. R., 1974. The impulse response of a Maxwell earth, Rev. Geophys. Space Phys., 12,649-669. Ramsey, A. S . , 1981. Newronian Attraction, Cambridge University Press, Cambridge. Richards, M. A. & Hager, B. H., 1984. Geoid anomalies in a dynamic earth, J . geophys. Res., 89, 5987-6002. Riimpker, G., 1990. Viskoelastische Perturbationen eines Burgers-Halbraums: Theorie und Anwendung, Diploma thesis, Westfalische Sabadini, R., Yuen, D. A. & Boschi, E., 1982. Polar wandering and the forced responses of a rotating, multilayered, viscoelastic planet, J. Slichter, L. B. & Caputo, M., 1960. Deformation of an earth model by surface pressures, J. geophys. Res., 65,4151-4156. Spada, G., Sabadini, R., Yuen, D. A. & Ricard, Y., 1992. Effects on post-glacial rebound from the hard rheology in the transition zone, Takeuchi, H., 1950. On the earth tide of the compressible earth of variable density and elasticity, Trans. Am. geophys. Un., 31, 651-689. Thomson, W., 1864. Dynamical problems regarding elastic spheroidal shells and spheroids of incompressible liquid, Phil. Trans. R. SOC. Lond., Wang, R., 1986. Das viskoelastische Verhalten der Erde auf langfristige Gezeitenterme, Diploma thesis, Christian-Albrechts-Universitat, Kiel. Wolf, D., 1984. The relaxation of spherical and flat Maxwell earth models and effects due to the presence of the lithosphere, J. Geophys., 56, Wolf, D., 1991 a. Viscoelastodynamics of a stratified, compressible planet: incremental field equations and short- and long-time asymptotes, Wolf, D., 1991b. Boussinesq’s problem of viscoelasticity, Terra Nova, 3, 401-407. Wolf, D., 1993. Gravitational viscoelastodynamics for a hydrostatic planet, Habilitation thesis, Westfalische Wilhelms-Universitat, Munster. Wu, P., 1990. Deformation of internal boundaries in a viscoelastic earth and topographic coupling between the mantle and the core, Geophys. Wu, P., 1992. Deformation of an incompressible viscoelastic flat earth with power-law creep: a finite element approach, Geophys. J . Int., 108, Wu, P. & Peltier, W. R., 1982. Viscous gravitational relaxation, Geophys. J . R. asrr. SOC., 70,435-485. excitation, J . geophys. Res., 88, 2240-2250. 67,845-850. numerical results, J. geophys. Res., 68, 485-496. Wilhelms-Universitat, Miinster. geophys. Res., 87,2885-2903. Geophys. J. Int., 109,683-700. 153, 583-616. 24-33. Geophys. J. Znt., 104,401-417. J. lnt., 101, 213-231. 35-51. A P P E N D I X A: MATHEMATICAL A N D N O T A T I O N A L CONCEPTS We modify and extend the concepts developed in Dahlen (1974), Grafarend (1982) and Wolf (1991a); a complete exposition of the following analysis can be found in Wolf (1993, pp. 4-10). A1 Kinematic formulations Suppose Cartesian-tensor fields and employ for them the usual indicia1 notation and summation convention. Denote by 8 the unbounded 3-D Euclidian space domain, by 5 the time domain [ O , x ) and consider the mapping ri = ri(X, r ) , where X i E 8 and t E 9. We stipulate that ri = r j ( X , t ) is a one-to-one mapping of 8 onto itself with the property X i = r f ( X , 0) and continuous differentiability with respect to the arguments as many times as required. We call t current time, ri current position, t = 0 initial time and Xi initial position. Assume now that % is completely filled by a gravitating fluid. A particular mapping satisfying our assumptions then is rlL) = rjyx, t ) = X i + U i ( X , t ) , xi E 8, t E 3, Lamk's problem 345 which is the kinematic formulation used for material points (particles). The mapping identifies each material point by its initial position, Xi, and relates to the point its current position, r!'), in terms of the material displacement, uir from its initial position. In addition, we define isopotential points, which move in the direction of the gradient of the gravitational potential currently existing at their respective positions such that the potential experienced by each of them remains constant during their motion. A second mapping satisfying our assumptions thus is riN) = r j N ) ( X , t ) = X, + d i ( x , t ) , xi E 8, t E 9, (289) which is the kinematic formulation used for isopotential points. The mapping identifies each isopotential point by its initial position, Xi, and relates to the point its current position, rj", in terms of the isopotential displacement, di, from its initial position. The inverse mappings to (288) and (289) are, respectively, In contrast to (288) and (289), eqs (290) and (291), refer to local points (places) identified by their position, ri. In particular, (290) relates to each local point the initial position, Xj"), of the material point currently at ri by means of the material displacement, Ui. Similarly, (291) relates to each local point the initial position, XIN', of the isopotential point currently at ri by means of the isopotential displacement, D,. In view of the general assumptions to be satisfied by the mappings r j L ) ( X , t ) and r j N ) ( X , t ) , (290) and (291) define one-to-one mappings that are continuously differentiable with respect to the arguments as many times as required. We proceed by specifying the domains of definition of the mappings (288)-(291) more closely. Beginning with (288), we decompose 8 into two open subdomains: the simply connected internal domain, gPL), and the complementary external domain, @:). With a@') the 2-D interface between the two domains, it then follows that 8 = @-L) U @:) U a@'). We now define B.L")(t) = { r$L) (X , t ) I xi E @:), r E T}, M P ) ( t ) = {r;')(X, t ) I xi E a@'), t E q. (292) (293) Considering the physical interpretation of (288), %F)(t) and 89?("(t) are the current domains of those material points initially occupying @,"I and a@'", respectively. In this study, we suppose that 9 y L ) ( t ) and 9?YL)(t) are domains of continuity for the parameters of the fluid and that a%?(t) is an interface of discontinuity for these parameters. We therefore take B.F)(f) = 9LN)(t), d%(')(t) = a9?("(t) and define @T)(t) = {XjN)(r, t ) I ri E %iN)(f ) , t E 9, aPN)( t ) = {XjN)(r, t ) I ri E ag("(t), t E 91. (294) (295) In view of the physical interpretation of (291), @T)(t) and a@"(r) are the initial domains of those isopotential points currently occupying %?zc,")(t) and 3B.(N)(t), respectively. Since 9?kL)(t) = %LN)(t) and d9?Z'L)(t) = 3 9 ( N ) ( t ) , no distinction is required and the symbols 9?+(t) and a9?(r) are used henceforth. Next, we give formulations equivalent to (288)-(291) for arbitrary field quantities. Since we wish to allow for the possibility that the values of such fields or their gradients are discontinuous on d % ( t ) , all material, isopotential and local points currently on this interface are excluded. We thus disregard material points for which X i E a@'), isopotential points for which Xi E a@"(t) and local points for which ri E a9?(t). With this, the generalizations to (288) and (289) are f$!. = f F ! . ( X , t ) . xi E a"'" U @p, t E 9, (296) f$N! =fkN!(X, t ) , xi E @-N)(t) u Pp(r), t E 9. (297) The quantity f:,!?. in (296) is the current value of an arbitrary field at the material point whose initial position is Xi. Similarly, fkN! in (297) is the current value of that field at the isopotential point whose initial position is X i . Eq. (296) is commonly referred to as Lagrangian formulation of the field. Eq. (297) is non-conventional and here referred to as Newtonian formulation. The generalization to (290) and (291) is I$... = F, ,... (r,r), r, E %(t ) LJ 9?+(t), t E 9. (298) This equation relates to each local point identified by its position, r,, the current value, cj..., of an arbitrary field at this point; it is commonly called Eulerian formulation of the field. The mappings defined in (296)-(298) are assumed to be single-valued and 346 D. Wolf continuously differentiable with respect to the arguments as many times as required. As in the preceding equations, we use lower case symbols for the Lagrangian and Newtonian formulations of fields and upper case symbols for the Eulerian formulation. The Lagrangian and Newtonian formulations are distinguished by the label superscripts L and N. A2 Perturbation equations We assume that the current value of an arbitrary field represents a perturbation of its initial value. Allowing for discontinuities of the field values on d%(t), the Newtonian and Eulerian formulations of the perturbation equation are then straightforward only for isopotential and local points that are initially in %-(O) and currently in %-(t) or that are initially in %+(O) and currently in %+(t). We call such points strictZy internal or external. For conciseness, we define $:)(ti = $c)(o) n @F)(t), (299) = %*(o) n %+(t). (300) On account of these equations, the necessary and sufficient conditions for strictly internal or external isopotential points and for strictly internal or external local points, respectively, are therefore Xi E @-”(t) U g?)(t) and ri E % ( t ) U i%+(t). Using this, the Lagrangian, Newtonian and Eulerian forms of the perturbation equation can be written as follows: fkL?.(X, t ) =f&!.(X, 0) + SfkL?.(X, t ) , f$N!(x, t ) =f$N’I!(x, 0) + afLN!(x, t ) , q j . . . ( r , t ) = cj,,.(r, 0) + A<,.,,(r, t ) , xi E @-=) u @!I, xi E @-”(t) u %!F)(t), rj E i%(t) U a+(?), t E q (301) (302) (303) t E z t E .T. We refer to the left-hand sides of the equations as total fields, to the first terms on the right-hand sides as initial fields and to the second terms on the right-hand sides as incremental fields. In particular, Sf $?!.(X, t ) is called material increment, df;N!(X, t ) isopotential increment and Aei ...( r, t ) local increment. In some neighbourhood of d%(t ) , isopotential and local points are initially in %-(O) and currently in %+(t) or vice versa. Since the field values are not necessarily continuous on a%(t), such hybrid points require special consideration. In order that this be avoided, we need the Lagrangian forms of (302) and (303). Using the abbreviation fF?.,,(X, t ) for the gradient of fkL?.(X, t ) with respect to X , and assuming infinitesimal perturbations, we obtain upon some algebraic manipulations (Wolf 1993, pp. 7-9) fF?.(X, t ) =fk?.(X, 0) + dfkL?.(X, t ) +fbL!,,(X, O)[u,(X, t ) - dlf-’(X, t ) ] , xi E 2FL’ u @?’, t E z f f? . (X, t )=f~~. (X,O)+Aff? , (X, t )+f~) ,k(X,O)u, (X, t ) , Xi E@-~)U@>), t E 3 For notational convenience, we adopt several simplifications: (i) the arguments Xi, ri and t are suppressed; (ii) the argument t = 0 is indicated by the label superscript 0 appended to the function symbols; (iii) the material, isopotential and local increments are indicated by the label superscripts 6, a and A appended to the function symbols. With these conventions, the three alternative forms of the Lagrangian perturbation equation, (301), (304) and (305), reduce to (306) xi E SL) u S:), t E r, (307) (308) fp =fp f f y f f?. = f k.0) + fF.A’ + f y i u , fg?. = f&.o’+f~.?”+ fEP,!(u, -d j f - ) ) whence The second terms on the right-hand sides of (309) and (310) are called aduectiue increments. They account for those parts of the increments resulting from the component of the motion of material or isopotential points parallel to the gradient of the initial field. In the present study, only the Lagrangian formulation is employed, allowing us to suppress L. A3 Interface .conditions We consider the behaviour of field values on a%. In order to formulate a condition expressing this behaviour, we locally assign to 8% (the Lagrangian form of) the unit normal directed outward into %+. Denoting this vector by nj and assuming E > O , Lamk's problem 347 we define The interface condition for fi,.. can then be written as [Lj ...I l=f; ...) X , € d % t € % where f$. .. is the increase of fi,... in the direction of n,. APPENDIX B: LIST OF IMPORTANT SYMBOLS B1 Latin symbols integration coefficient of fundamental solution to (4 x 4) system radius of planet integration coefficient of fundamental solution to (2 x 2) system isopotential displacement 2.71828.. . strain normalized Legendre coefficient off r-dependent part of Fn scalar field or function Laplace transform of f non-gravitational force per unit mass ordinary Legendre coefficient off axisymmetric Green's function for f non-symmetric Green's function for f Cartesian tensor field of arbitrary rank gradient of j& with respect to Xk initial value of J j . . . local increment off$?. material increment off 1;9!. isopotential increment of f:!. Newton's gravitational constant gravity (gravitational force per unit mass) isopotential height sequential number of fundamental solution to (4 X 4) system Legendre wave number B2 Greek symbols inverse spectral time inverse elemental Maxwell time inverse modal relaxation time magnitude of gjo) on da" Dirac delta function Kronecker symbol partial-derivative operator with respect to t Levi-Civita symbol non-dimensional Legendre coefficient of CT colatitude of observation point colatitude of excitation point angle between observation and excitation points eigenvalue of (4 X 4) system c (313) sequential number of fundamental solution to (2 X 2) system total number of relaxation modes sequential number of relaxation mode Legendre degree outward unit normal on 8% origin of coordinate system Legendre polynomial of the first kind (mechanical) pressure total number of Maxwell elements sequential number of Maxwell. element non-dimensional radial distance of observation point radial distance of observation point current position of material point inverse Laplace time deviatoric incremental stress impulse-response function transfer function current time excitation time (Cauchy) stress material displacement initial position of material point solution vector of (4 X 4) system eigenvector of (4 X 4) system solution vector of (2 X 2) system eigenvector of (2 X 2) system longitude of observation point longitude of excitation point eigenvalue of (2 x 2) system shear-relaxation function shear-relaxation spectrum elastic shear modulus elemental elastic shear modulus 3.14159.. . volume-mass density (incremental) interface-mass density gravitational potential rotation 348 D. Wolf B3 Calligraphic symbols % Euclidian space domain 9- internal domain of ri 9+ external domain of ri 9’ domain of s 9 domain of t 2- internal domain of X , 2+ external domain of X, 8% 8% interface between 9- and %+ interface between %- and %+