## The Modern Global Theory Of The Glacial Isostatic Adjustment Process

The theory of the GIA process is embodied in an integral "sea level equation," the solution of which describes the time-dependent separation between the surface of the solid earth and the gravitational equipotential surface that determines the equilibrium level of the sea. This equipotential surface is, of course, precisely the geoid of classical geodesy. The basic ingredients and initial form of this sea level equation (SLE) were developed in Peltier (1974, 1976), Peltier and Andrews (1976), and Farrell and Clark (1976). First solutions were published by Clark et al. (1978) and Peltier et al. (1978). This rather primitive initial theoretical structure has since been very significantly refined both in regard to the mathematical methods employed to solve the SLE and in regard to the range of physical effects included as influences upon the solution. The most general form of the SLE which governs the time dependence of relative sea level subsequent to the onset of any variation of land ice cover is as follows (see Peltier, 1998a for a most up-to-date review and detailed discussion), with S(6, A, t) the variation of relative sea level at latitude 6 and longitude A as a function of time t:

In Eq. (1) the function C(0, A, t) is the so-called "ocean function," which is, by definition, equal to unity where there is ocean and zero elsewhere. This equation is an integral equation because the surface mass load per unit area, L, is composed of both ocean-water and land-ice components. It takes the explicit form

in which 1(8, A, t) is the space-time history of ice-thickness variations and S(d, A, t) is the space-time variation of ocean-water thickness variations. The latter is, of course, simply relative sea level itself. Due to the composite property of L expressed in Eq. (2), the unknown observable function S appears both on the left-hand side of (1) and under the triple convolution operator on the right-hand side of (1). The theory embodied in (1) is thus of integral equation form and may be employed to determine S given an assumed known history of glaciation and deglaciation 1(9, A, t), together with suitable impulse response functions (Green functions) and Gj. These Green functions embody all of the information concerning the radial viscoelastic structure of the earth that is required to determine the impact of this structure upon relative sea level history. In terms of the surface load and tidal "Love numbers," (hL, kL) and (hT, k[ ), respectively, these Green functions have the following mathematical forms, which consist of infinite series expansions in Legendre polynomials Pf(cos 6) in which f is spherical harmonic degree:

G^ (6, A, t) = — 2 (1 + # (0 - h\ (r)) P( (cos 6), (3a)

Gl (6, A, 0 = - 2 (1 + kl (t) - hj (r)) Pt (cos 6). (3b)

Sf=0

The Love numbers in turn have Dirichlet series expansions of the form (Peltier 1976)

In these Dirichet series the hf'E and kY~F are precisely the elastic surface load Love numbers of Farrell (1972) which measure the magnitude of the response to application of a load that varies as a delta function in time. The quantity S(t) is the Dirac delta function, the if are a set of normal mode "poles" in the plane of the complex Laplace transform variable 'V representing the inverse relaxation times of modes of "viscous gravitational relaxation" (Peltier, 1976), and the rj and qj are the residues at these poles that may be determined in the usual way using Cauchy's residue theorem (Peltier, 1985) or the collocation technique described previously in Peltier (1974). The tidal Love number counterparts of {h\, k\), namely (hj, kj), have expansions of precisely the same form as (4a) and (4b) but they are determined subject to the boundary conditions that are appropriate when the responses arise through an imposed variation in the gravitational potential field rather than through an imposed surface mass load (see Farrell, 1972, for details).

Solutions to the sea level equation (1) are constructed using an iterative methodology that is based upon the recognition that the influence of time variations in the ocean function C(6, A, t) and the influence of the feedback onto sea level of the changing rotational state of the planet due to the variations of surface mass load described in (1) by the convolution of with Gj are both second-order effects. The complete algorithm that I have developed to incorporate these effects, as discussed in Peltier (1994, 1998a, 1999), begins with the construction of a solution to (1) which neglects both of these effects and which is derived using the semispectral algorithm described in Mitrovica and Peltier (1991)—an algorithm that delivers a solution S(0, A, t) in terms of the time-dependent coefficients in an expansion of the solution on a basis of spherical harmonics. Solutions are now typically computed at very high spatial resolution by truncating the spherical harmonic expansions to degree and order 512. Given this first approximation to the full solution, we have complete knowledge of the variations of surface mass load associated with the glacial cycle, the component 1(6, A, t) having been an assumed known input to (1) that may be varied to improve the fit to the observations and the component S(6, A, t) being the output from its solution. The surface load L(6, A, t) in (2) is thus entirely determined.

In the second step of this iterative process we next compute the rotational response to this history of surface loading by solving the Euler equation,

in which Jij is the moment of inertia tensor of the planet, oj, are the components of its angular velocity vector, and eijk is the Levi-Cevita alternating tensor. Assuming a biaxial form for the undeformed shaped of the planet, accurate solutions to (5) may be constructed by employing the standard perturbation expansion

 CO, = CI (8ij + m,) ^ = Iij, i ^ j Jn = A + In (6) J22 = A + I22 J33 = C + 733,

in which (A, A, C) are the principal moments of inertia, Q is the basic state angular velocity of the earth, and and m, are (assumed small) fluctuations away from the basic state values. To first order in these fluctuations (Munk and MacDonald, 1960), substitution of (6) into (5) leads to the following decoupled system for polar motion and earth rotation, respectively, as

Pr V

in which the overdot indicates a time derivative, the so-called excitation functions are and ^3, <rr = fl (C - A)!A is the Chandler wobble frequency of the rigid earth, m = wii + 1 w2, ^ = + i i = and the are

Now the functions I,j and thus the excitation functions ^ are entirely determined by the surface loading history L(6, A, t) in Eq. (2), which is known from the first step in the iterative procedure [the details of these relationships are discussed in Peltier (1982) and Wu and Peltier (1984)]. Because the changes in the moments of inertia induced by surface loading also depend upon the redistribution of mass in the interior of the planet that is induced by the surface loading process, and because the redistribution of "internal mass" also depends upon the (assumed radial) viscoelastic structure, the excitation functions also depend upon this structure. Peltier (1982) and Wu and Peltier (1984) show in detail how to employ Laplace transform methods to construct complete solutions to Eqs. (7a) and (7b) for the time dependence of the perturbations to the angular velocity vector of the planet û}i(t), i = 1,3. No useful purpose would be served by repeating here the details of the algebraically complex analyses that are required to construct this solution for the rotational response to déglaciation.

In the third step of the iterative procedure we next construct the forcing function in (1) required to compute the impact upon sea level due to this changing rotational state. Because the solution of Peltier (1982) and Wu and Peltier (1984) used to determine the w,(i) is based upon the assumption that wild <§ 1, this same assumption must be made in computing Dahlen (1976), in the context of his analysis of a different problem, has shown that the appropriate general solution for correct to first order in perturbation theory, has the spherical harmonic expansion

in which the coefficients are related to the w,(i) by

Given R we may now complete the third step in the iterative procedure by solving (1) once more for S(0, A, f), this time incorporating the influence of the changing rotational state of the planet, the first of the two second-order effects that the iterative procedure is intended to capture.

Now these first three steps in the iterative procedure have been based upon the assumption that the ocean function C(8, A, t) was constant and equal to its present-day form. The fourth step in the iterative procedure involves relaxing this assumption and this may be accomplished by following the procedure first described in Peltier (1994), which exploits the fact that the sea level equation (1), being itself a construct of first-order perturbation theory, delivers a solution for relative sea level S(6, A, t) that is "relative" to an unspecified and therefore arbitrary datum. This arbitrariness may be exploited to construct what I have previously referred to as a "topographically self-consistent" solution for RSL history. The procedure is simply as follows: First, we define a time-independent field T'(0, A) such that when we add it to the solution of the SLE computed for the present day (time t = ip), we obtain the known present-day topography of the planet with respect to sea level, as

in which Tp is the present-day topography wrt sea level, determined, say, by the ET0P05 data set. Using (11) to define T'(0, A) we may then determine the topography at any time in the past by computing

T(0, A, 0 = S(6, A, t) + [Tp(6, A) - S(9, A, ip)]. (12)

The true paleotopography, including the contribution from the ice sheets of thickness 1(0, A, t), is then

At any point in space and at any instant of time for which PT < 0 we therefore have ocean, whereas "land" is where PT > 0 (the only exception to this general rule concerns a number of low-lying inland seas (such as Caspian Sea) whose surfaces clearly do not lie on the same equipotential surface as does the global ocean). Subject to this caveat we may define a first approximation to the actual time-dependent ocean function as

In reconstructing the topogarphy PT it is also important to recognize the "implicit" component of the ice load (Peltier 1998c) that is activated in the solution of (1) because of the perturbation theory-based nature of this equation. Since the fact that the redefinition of land to sea that occurs (for example) when the initially ice-covered Hudson Bay, Barents and Kara Seas, and the Gulf of Bothnia become connected to the ocean is a nonperturbative effect, the "explicit" component of the ice load that is varied in (1) to enable the model to fit the data must be considered to be an "effective" load. To compute the actual load and thus the true paleotopography, one must add to the effective load the additional "implicit" component. This completes step 4 of the iterative procedure. The next step in the process involves executing the first four steps in the process a second time, ending with a second guess for the time-dependent ocean function C2(0, A, t) and, of course, with a new solution for relative sea level history S(6, A, t) that incorporates the full influence of rotational feedback.

Before discussing the impact of the global process of glacial isostatic adjustment on modern tide gauge-based estimates of secular sea level change, in the next section I will first discuss the global form of the solutions which the above-discussed theory delivers for the present-day rate of relative sea level rise and the extent to which it is able to provide satisfactory fits to a globally distributed set of relative sea level history observations. It will thereafter prove interesting not only to discuss the sea level signature that a specific tide gauge might be expected to measure but also to consider the signature of the time rate of change of absolute sea level (geoid height) that would be observed by an artificial Earth-orbiting satellite, either one equipped with an altimeter to measure distance from the satellite to the sea surface (e.g., TOPEX/POSEI-DON; see Chapter 6 for details of the way in which this system has already been employed to address the issue of global sea level change) or one that may be employed to infer this signal through measurement of the time dependence of the space-dependent strength of the gravitational field itself (expressed in terms of the time dependence of geoid height). The latter measurement is one that will be made at very high resolution in the near future using the GRACE satellite, and the CHAMP precursor, as discussed in greater detail in Peltier (1999).

0 0