Abstract

We present a generalized formalism for computing gravitationally self-consistent sea level changes driven by the combined effects of dynamic topography, geoid perturbations due to mantle convection, ice mass fluctuations and sediment redistribution on a deforming Earth. Our mathematical treatment conserves mass of the surface (ice plus ocean) load and the solid Earth. Moreover, it takes precise account of shoreline migration and the associated ocean loading. The new formalism avoids a variety of approximations adopted in previous models of sea level change driven by dynamic topography, including the assumption that a spatially fixed isostatic amplification of ‘air-loaded’ dynamic topography accurately accounts for ocean loading effects. While our approach is valid for Earth models of arbitrary complexity, we present numerical results for a set of simple cases in which a pattern of dynamic topography is imposed, the response to surface mass loading assumes that Earth structure varies only with depth and that isostatic equilibrium is maintained at all times. These calculations, involving fluid Love number theory, indicate that the largest errors in previous predictions of sea level change driven by dynamic topography occur in regions of shoreline migration, and thus in the vicinity of most geological markers of ancient sea level. We conclude that a gravitationally self-consistent treatment of long-term sea level change is necessary in any effort to use such geological markers to estimate ancient ice volumes.

1 INTRODUCTION

The earliest discussions of plate tectonics focused on horizontal motions at the Earth's surface, reflecting a preoccupation with the concepts of continental drift and seafloor spreading. However, the driving force for these motions—thermochemical convection within the Earth's mantle—is also responsible for vertical deflections of the Earth's surface (McKenzie 1977) that have come to be known as dynamic topography (e.g. Hager et al.1985). The earliest simulations of time-varying dynamic topography connected the process to epeirogenic motions of continents and long-term regional- and global-scale sea level changes, and typically focused on subduction-controlled dynamics (e.g. Mitrovica et al.1989; Gurnis 1990, 1992, 1993; Holt & Stern 1994; Coakley & Gurnis 1995; Mitrovica et al.1996; Lithgow-Bertelloni & Gurnis 1997; Pysklywec & Mitrovica 1998; Gurnis et al.1998; Pysklywec & Mitrovica 1999). However, several studies also demonstrated a clear link between upwelling mantle flow and large-scale fluctuations in the topography of continents (e.g. Lithgow-Bertelloni & Silver 1998; Gurnis et al.2000; Forte et al.2010).

Predictions of present-day, global-scale dynamic topography have been derived from viscous flow modelling constrained by seismic inferences of mantle heterogeneity (e.g. Hager et al.1985; Ricard & Wuming 1991; Forte et al.1993; Wen & Anderson 1997; Steinberger 2007). As the resolution of seismic imaging has progressively improved, viscous flow modelling of time-varying dynamic topography has shown concomitant advances in resolution, allowing applications to the detailed, regional-scale geological record of continents (e.g. Daradich et al.2003; Moucha et al.2009; Rowley et al.2013; Liu et al.2014). These studies have had fundamental implications for our understanding of long-term sea level change. For example, they have highlighted the fact that dynamic topography is present globally, not only at plate boundaries, and that it can significantly contaminate records of long-term sea level change at all locations, even passive continental margins (e.g. Moucha et al.2008; Müller et al.2008; Conrad & Husson 2009; Spasojevic & Gurnis 2012; Rowley et al.2013). Moreover, these analyses have shown that the globally averaged change in the shape of ocean basins is governed by the integrated effect of dynamic uplift and subsidence throughout these basins, not just changes in elevation linked to fluctuations in the rate of plate creation.

While the connection between dynamic topography and regional and global sea level changes has been established, a rigorous mathematical and numerical formalism for accurately mapping viscous flow predictions of time-dependent dynamic topography into geographically variable sea level change is lacking. Previous efforts to compute this mapping have several shortcomings:

  1. A change in dynamic topography will redistribute water, and the loading (or unloading) associated with this redistribution will, in turn, drive displacements of the crust and hence cause sea level changes. This loading effect is generally treated in an approximate sense through isostatic amplification of the ‘air-loaded’ dynamic topography. In this regard, some studies (e.g. Moucha et al.2008; Rowley et al.2013) assume water cover everywhere and apply a global amplification to the air-loaded dynamic topography obtained from their viscous flow modelling. In contrast, Müller et al. (2008) applied no amplification to the dynamic topography change for their predictions along the U.S. East Coast. A hybrid approach only amplifies the dynamic topography at sites with ocean cover (Steinberger 2007; Flament et al.2013), however, it and the other approaches neglect water loading within any region of shoreline migration. This assumption is problematic given that the error introduced into the prediction will be largest in the vicinity of an ancient sea level marker.

  2. The above studies also ignore the possibility of flexural effects in the response to water loading, that is, they assume local hydrostatic equilibrium in the crustal response. Flexure has a direct effect on the computed sea level change via its impact on crustal elevation (e.g. the crustal displacement will extend outside the region of loading). Moreover, it has indirect effects on sea level change because any departures from hydrostatic equilibrium will drive gravitational and rotational effects that feedback into sea level.

  3. Studies of dynamic topography are often only concerned with the change in crustal elevation due to mantle stresses. Sea level, however, is the difference between the elevation of both the gravitational equipotential defining the sea surface and the crust, and therefore changes in the gravitational field that arise from mantle convection will perturb sea level, an effect only rarely taken into account (Spasojevic & Gurnis 2012; Conrad & Husson 2009).

  4. A further shortcoming involves the assumption of a local application of isostasy (e.g. Müller et al.2008; Conrad & Husson 2009; Spasojevic & Gurnis 2012). Any loading-induced uplift (subsidence) of the crust must be matched by subsidence (uplift) elsewhere to conserve mass of the solid Earth. Local treatments of isostasy (whether they incorporate flexure or not) ignore this global-scale compensation.

  5. Finally, predictions of long-term sea level change should also conserve the mass of the ocean (or ice plus ocean), but this constraint is not always applied in such calculations (e.g. Flament et al.2013).

In this article, we present a gravitationally self-consistent mathematical treatment of global-scale sea level change driven by time-varying dynamic topography, convection-driven geoid perturbations, ice mass changes and sediment redistribution. Our approach is based on a formalism developed in studies of ice age sea level change (Farrell & Clark 1976; Mitrovica & Milne 2003; Kendall et al.2005; Dalca et al.2013) and it conserves mass of both the surface (ice plus ocean) load and the solid Earth while taking precise account of the migration of shorelines. In its most general form, and under the assumption that the dynamic topography, convection-driven geoid perturbations, ice mass and sediment histories are known, the method is valid for viscoelastic Earth models of arbitrary complexity; one only needs to have a numerical method for computing the response of the Earth model to an evolving surface mass load. Our approach may be used to consider dynamic topography in isolation; however, we have included the possibility of contemporaneous ice mass changes and sediment fluxes to provide a framework that avoids errors that may be incurred by assuming that the total sea level change can be computed by treating each of the driving processes in isolation. Finally, we also outline the extension necessary when the timescale being considered is sufficiently long that plate motions must be taken into account.

In addition, we present successive, special cases that assume: (1) isostatic equilibrium is maintained throughout the evolution of the system and (2) the Earth model varies with depth alone. The first assumption simplifies the numerical calculations by ignoring time-dependent viscous adjustments and this allows one to evolve the system backwards in time starting from a known (present-day) topography. Taken together, the two assumptions allow us to make use of fluid Love number theory to compute the response of the Earth model to the surface mass loading. In this case, we adopt profiles of density and elastic moduli from the PREM seismic model (Dziewonski & Anderson 1981) and the thickness of the elastic lithosphere is a free parameter in the modelling. To illustrate the new, gravitationally self-consistent treatment of sea level change driven by dynamic topography, we present a series of numerical predictions based on these special cases.

2 CALCULATING LONG-TERM SEA LEVEL CHANGE

2.1 Calculating dynamic topography from mantle convection

A variety of approaches have been adopted to compute dynamic topography in mantle convection simulations. The most rigorous but computationally expensive approach is to allow the surface to deform freely. These free surface calculations allow the upper boundary of the computational domain to respond dynamically to evolving stresses in the mantle (Gurnis et al.1996; Zhong et al.1996; Kaus et al.2010; Kramer et al.2012). An alternate, fixed surface approach introduces a top layer in the computational domain with a comparatively low viscosity and low density, which serves as a proxy for air and yields quasi-free surface deformation. This approach is often referred to as the ‘sticky-air’ method (Schmeling et al.2008; Crameri et al.2012). The third and most common approach models mantle convection within a fixed domain and computes the radial stresses at the top of the domain. The dynamic topography (DT) is then calculated as the compensation height that balances these radial stresses (McKenzie 1977; Hager et al.1985; Zhong et al.1993; Flament et al.2013):
(1)
where σrr is the radial stress at the upper surface of the mantle flow model, Δρ is the density contrast between the crust and the overlying material (e.g. air, water or sediment) and g is the gravitational acceleration. This approach can be improved by including self-gravitation associated with the disturbed topography (Zhong et al.2008) or flexural effects of an elastic lithosphere in response to the internal deformation (Golle et al.2012). The three distinct approaches yield comparable predictions for the dynamic topography (Crameri et al.2012); the last approach is characterized by a slight shift backwards in time, compared to the others, because it assumes an instantaneous isostatically adjusted response of the surface to internal loading.

As is clear from the basic principles of isostasy that underlie eq. (1), the amplitude of predicted dynamic topography arising from viscous stresses associated with mantle flow increases as the density difference between the crust and overlying material decreases. Thus, the choice of an appropriate density contrast is a key issue in comparing predictions of dynamic topography with the geological record. Regardless of the method employed for modelling dynamic topography, one has to choose an appropriate density contrast between the crust and the overlying material. As described in Section 1, this is a shortcoming of existing models of dynamic topography and it can be overcome using the framework discussed below.

2.2 Sea level changes in the presence of dynamic topography

The derivation presented below is based on the ice age sea level theory developed by Mitrovica & Milne (2003), Kendall et al. (2005) and Dalca et al. (2013), following the canonical treatment of Farrell & Clark (1976). Farrell & Clark (1976) derived the so-called ‘sea level equation’ governing gravitationally self-consistent sea level changes on a viscoelastic Earth driven by global ice mass variations. Their derivation assumed that the position of shorelines was fixed in time, which is equivalent to assuming that all shorelines are characterized by steep vertical cliffs such that any local changes in sea level produce no transgression or regression (i.e. no shoreline migration). Their calculations were based on spherically symmetric, Maxwell viscoelastic and non-rotating Earth models, and in this case they adopted viscoelastic Love number theory (Peltier 1974; Wu & Peltier 1982) to compute time-dependent perturbations in the elevation of the solid surface and gravitational field. Mitrovica & Milne (2003) extended the sea level equation to accurately incorporate shoreline migration and Kendall et al. (2005) presented an algorithm for solving the equation in the case of 1-D and 3-D Earth models of arbitrary complexity. The Mitrovica & Milne (2003) treatment allows for the feedback of perturbations in Earth rotation into sea level, and in the case of a 1-D, rotating Earth model, Kendall et al. (2005) derived the necessary equations making use of viscoelastic Love number theory. The Dalca et al. (2013) study extended these treatments to incorporate sediment redistribution into the gravitationally self-consistent sea level calculation. Their form of the theory, which defined the topography as including ice and sediment, is advantageous to our treatment of dynamic topography, and so it forms the basis of our derivation.

We adopt the following symbolism:

R - Earth's bedrock elevation;

H - Sediment thickness;

I - Thickness of land-based and grounded marine-based ice;

G - Height of the equipotential surface that defines the sea surface.

All these quantities have a dependence on position, expressed in terms of colatitude θ and east longitude φ, and time t. We will henceforth suppress these explicit dependencies to simplify the expressions. In this regard, we abbreviate a generic function χ(θ, φ, tj) by χ, or χj when it is necessary to make the time dependence explicit.

Sea level is defined as the difference between the heights of two bounding surfaces: the sea surface equipotential and the solid surface (including the height of ice and/or sediments), where the latter is given by R + H + I:
(2)
In this case, topography is simply defined as the negative of sea level:
(3)
Moreover, the ocean height is simply the projection of sea level onto the location of oceans:
(4)
where C is the ocean function defined by:
(5)
The total surface mass load L is comprised of water (oceans), ice and sediments:
(6)
where, ρW, ρI and ρH are the densities of water, ice and sediment, respectively.
Geological markers record local changes in sea level, as defined in eq. (2), and thus we are interested in how changes in grounded ice, sediments, water and dynamic topography, perturb the elevation of the bedrock and sea surface height from some initial state. Any time-dependent field χ at some time tj can be expressed as a sum of the initial value at the onset of loading, at time t0, and a perturbation from this state:
(7)
This form may be applied to any of the quantities R, H, I, G, SL, T, S, C and L. As an example, the changing bedrock elevation may be written as:
(8)
where the term R0 includes the background signal from dynamic topography. The perturbation of the bedrock elevation can furthermore be decomposed into a time-dependent perturbation in air-loaded dynamic topography plus a perturbation due to the response of the bedrock to the changing surface mass (ice plus ocean plus sediment) load. If we denote these terms by ΔDTA and ΔRL, respectively, we have:
(9)

It is important to emphasize that it is air-loaded dynamic topography that should appear in this equation. The amplification from air-loaded to water-loaded dynamic topography at sites under water, as indicated by eq. (1), is naturally incorporated into eq. (9) (as well as the broader sea level calculation) via the term |$\Delta R^L_j$|⁠. That is, an increase in the air-loaded dynamic topography (i.e. an uplift) at sites within the oceans will displace water and this unloading will act to amplify the uplift, as required. Similarly a dynamic subsidence computed in the air-loaded case at such sites will increase the local water load and amplify the subsidence. To be consistent, the air-loaded dynamic topography change and the load-induced perturbation to the bedrock elevation (the two terms on the right-hand-side of eq. (9)) should be computed using the same Earth structure model.

Using the above expressions, the change in sea level is given by:
(10)
The change in ocean height S cannot simply be written as a combination of perturbations in G, R, I and H because it involves the projection onto the ocean function given by eq. (4), which is itself time-dependent. Using this equation, together with the general form (7) and the simple relationship between sea level and topography (3) yields:
(11)

Eqs (3), (5), (10) and (11) represent the generalized form of the sea level equation valid for any timescale of forcing and Earth models of arbitrary complexity. Physically, eq. (11) indicates that the ocean height change at time tj is equal to the change in sea level projected onto the ocean function, both at time tj, plus a correction term that accounts for shoreline migration (Mitrovica & Milne 2003). The latter involves a projection of the initial topography onto a field (ΔCj) that is only non-zero in regions experiencing shoreline migration between times t0 and tj (i.e. ΔCj is only non-zero when C0Cj).

It is also useful to have an expression for the change in the surface mass load from its initial state. Using eqs (6) and (7), this change is given by
(12)
The change in the water load driven by dynamic topography is embedded in both fields on the right-hand-side of eq. (4). Sea level, SL, is dependent on dynamic topography via eq. (2) (and eqs 8 and 9), and this dependence also alters the ocean function C via eq. (5). The latter demonstrates that ocean mass redistribution driven by dynamic topography will, in general, lead to shoreline migration, even at shorelines which experience no local change in dynamic topography. The sea level equation incorporates this physics while maintaining gravitational self-consistency.
The sea surface equipotential G is perturbed by both the surface mass loading and by mantle convection. We could have incorporated this second contribution into our definition for dynamic topography; however, models of mantle flow-induced dynamic topography typically define dynamic topography as an absolute measure of bedrock height (or as being relative to the centre of mass of the Earth), not as a height relative to the geoid (or sea surface equipotential). The perturbation in the height of the sea-surface may be expressed as:
(13)
where the first term on the right-hand-side denotes the perturbation due to mantle convection and the second term is the sea surface equipotential due to surface mass loading. The former is the perturbation to the geoid that mantle flow modellers often compute in tandem with dynamic topography.

Eqs (10)–(12) highlight the integral nature of the sea level eq. (11). In particular, from eqs (9), (10), (11) and (13), the change in ocean height, ΔSj, depends on surface load-induced perturbations to the height of the bedrock (⁠|$\Delta R^L_j$|⁠) and sea surface potential (⁠|$\Delta G^L_j$|⁠). However, both |$\Delta R^L_j$| and |$\Delta G^L_j$| are, in turn, dependent on the ocean height change since the latter forms a component of the total surface mass load given by eq. (12).

It is traditional in gravitationally self-consistent ice age sea level theory to decompose the right-hand side of eq. (13) into geographic variable and globally uniform terms (Farrell & Clark 1976). To this end, we can write:
(14)
where |$\Delta \mathcal {G}^{mc}$| and |$\Delta \mathcal {G}^L$| are the geographically variable components of ΔGmc and ΔGL, respectively. While the sea surface must remain an equipotential surface, the specific equipotential that defines the sea surface need not remain constant over time (Farrell & Clark 1976; Dahlen 1976). Indeed, changes in the mass of the ocean in response to variations in global ice volumes, as well as changes in the shape of the ocean basin associated with load-induced deformation or the redistribution of sediments, will perturb the equipotential that defines the sea surface. In eq. (14), this time-dependent, globally uniform perturbation is denoted by ΔΦj, and dividing this value by the surface gravitational acceleration, g, converts the perturbation into a geographically uniform perturbation in the height of the initial sea surface equipotential.
Using eq. (14) into eq. (10) allows us to decompose the sea level change as follows:
(15)
where
(16)
and thus the sea level equation (11) can be rewritten as
(17)
where |$\Delta \mathcal {SL}_j$| is given by eq. (16) and Cj by eq. (5).
The term ΔΦj/g can be solved for by invoking conservation of mass of the ocean plus ice reservoirs, which requires that
(18)
where Ω spans the full surface area of the Earth. Integrating eq. (17) over the Earth's surface, using eq. (18) and solving for ΔΦj/g yields:
(19)
In this expression,
(20)

The sea level eq. (17), together with the set of eqs (3), (5), (16), (19) and (20), governs gravitationally self-consistent sea level changes in response to mantle-convection-induced deformation of the Earth's solid surface (dynamic topography) and gravitational field, ice mass variations and sediment redistribution.

The time-dependent fields ΔDTA, |$\Delta \mathcal {G}^{mc}$|⁠, ΔI, and ΔH are provided on input, and the full solution yields the fields ΔRL, |$\Delta \mathcal {G}^L$| and ΔC, the (time-dependent) scalar ΔΦ and, both the global sea level change ΔSL (via eq. 15) and the ocean height change ΔS. The global sea level change and the ocean height change are related through the sea level eq. (17), and the fact that the former is in general non-zero over land has an important physical interpretation, as elucidated by Dahlen (1976). In particular, at some specific site (θ, φ) on land at time t, the value ΔSL(θ, φ, t) represents the change in the local depth of water that would be measured if one were to dig a deep, infinitesimally thin canal connecting the site to the ocean.

As we noted above, the sea level equation is an integral equation since the quantity being solved for, ocean height change, is a component of the surface mass loading. In practice, this complication is generally addressed numerically by solving the sea level equation (17) iteratively, whereby at each time step an initial guess of the ocean height change is input into the right-hand-side of the equation, which is then solved to yield an updated estimate of ΔSj, and the process is repeated until convergence. One then moves to the next time step and the process is repeated. A reasonable first guess to the solution is to assume that the load-induced perturbation in the elevation of the bedrock and sea-surface are equal to zero. In this case:
(21)
and
(22)
with
(23)
and
(24)
In these equations, the superscript i denotes the iteration counter.
The fields Ij and ΔIj represent grounded ice height. In solving the sea level equation, a check is made at each time step to ensure that the ice thickness in marine sectors is sufficient to ensure that the ice is grounded, and if it is not, the ice model is revised to remove any floating ice (Kendall et al.2005). If we denote the input model ice history as |$I^{{\rm model}}_j$|⁠, then this check is based on the following criterion:
(25)

The sea level equation (17) requires that the initial topography T0 is specified. In traditional ice age applications, T0 is computed through a second iterative procedure. In particular, one begins with a first guess for T0 and this is used to compute sea level changes from the beginning of the simulation to the present day. The difference between the computed present-day topography and the observed topography is then used to update the estimate of T0 and the process is repeated until convergence, that is, until the computed present-day topography matches the observed topography to within a specified tolerance. This ‘outer’ iteration, which loops over the full time span of the simulation, should, in general, be adopted in the present application, however, we note an important exception below.

The set of equations we have derived here are valid for arbitrarily complex Earth models. Regardless of the level of complexity, the main requirement is a method for computing perturbations in the elevation of the Earth's bedrock surface and sea surface, ΔRL and |$\Delta \mathcal {G}^L$| respectively, driven by surface mass loading. For models with 3-D mantle structure, ΔRL and |$\Delta \mathcal {G}^L$| can be computed using a variety of numerical techniques developed for problems in glacial isostatic adjustment (e.g. Wu & van der Wal 2003; Zhong et al.2003; Latychev et al.2005).

2.2.1 Some special cases

There are a variety of special cases one might consider in computing sea level changes driven by dynamic topography. For example, if one adopts a spherically symmetric Earth model, then the equations can be cast in the spectral domain making use of viscoelastic Love number theory (Peltier 1974; Wu & Peltier 1982). The relevant equations are given in Dalca et al. (2013).

For long timescale applications that are generally of interest when considering changes in dynamic topography, one might additionally assume that the Earth system reaches full isostatic equilibrium across each time step in which the dynamic topography is updated. That is, one may assume that these time steps are long enough that all viscous stresses in response to the surface loading have relaxed and isostatic equilibrium has been reached. Departures from isostatic equilibrium may arise, for example, if the viscoelastic lithosphere has a relaxation time comparable to the timescale of the convective forcing. Furthermore, internal density discontinuities contribute buoyancy in response to the convective flow for timescales of order several million years (Piromallo et al.1997). Isostatic adjustment in response to a surface mass load is also impacted by these internal discontinuities, and thus reaching complete isostatic equilibrium would require a similarly long timescale. However, the departure from the equilibrium state would be relatively small for shorter timescales since the viscous modes associated with the internal discontinuities only contribute at the percent level to the load-induced surface displacement.

The assumption of isostatic equilibrium further simplifies the algorithm by making the ‘outer’ iteration (to find T0) unnecessary. In this case, since time-dependent viscous effects are ignored, the simulation is reversible in time and T0 can be interpreted as the present-day topography, which is known. The simulation then steps backwards in time through the history of ice, sediment, dynamic topography and gravity fields. The assumption of isostatic equilibrium also allows one to evolve the system directly from t = 0 to tj without calculating intermediate time steps.

In the examples discussed below, we assume isostatic equilibrium and a 1-D Earth model with a purely elastic lithosphere of specified thickness. In this case, the viscoelastic Love number theory appropriate for a 1-D Earth can be simplified to use the long-term, fluid limit of the Love numbers (Wu & Peltier 1982). The relevant equations for this case are provided in the Appendix. The Earth's viscoelastic lithosphere is characterized by a suite of relaxation times. Thus, the assumption that the lithosphere retains any elastic strength is only valid if the timescale of convective flow is much shorter than the longest of these relaxation times (in this case, the lithospheric thickness we are adopting in the numerical model should be interpreted as an effective elastic thickness associated with the unrelaxed lithospheric modes). Furthermore, the assumption of isostatic equilibrium requires that all other viscous relaxation is complete (i.e. that the timescale of flow is much longer than the other relaxation times associated with the viscoelastic lithosphere). Finally, the adoption of a 1-D Earth model ignores important complexity associated with lateral variations in lithospheric strength, including the presence of weak zones at plate boundaries.

While we make the assumption that the adoption of fluid Love number theory is valid in the results presented below, we emphasize that a time-forward and iterative procedure is necessary in calculations in which viscous effects cannot be ignored, as is the case for Pleistocene ice age calculations where changes in ice and ocean loading occur on timescales such that glacial isostatic adjustment must be modelled as a viscoelastic phenomenon.

2.3 Schematic illustration of sea level physics in the presence of dynamic topography

Fig. 1 provides a schematic illustration of the physics of sea level change near a shoreline in the presence of the combined effects of a local perturbation in dynamic topography and a lowering of the sea-surface height (neglecting for this illustration the flexural effects of the elastic lithosphere). The latter may arise from a variety of processes, including, for example, the growth of ice cover at some location distant from the shoreline. Figs 1(a)–(c) decompose the total sea level change into a two-step evolution, the first associated with the perturbation in dynamic topography and the second from the lowering of the sea surface height. If the perturbation in dynamic topography and sea surface height occurs in tandem (i.e. in one step), then the sea level change is given by the evolution from Fig. 1(a) to (directly) Fig. 1(d).

Schematic illustration of the physics of sea level change near a shoreline subject to dynamic topography and a change in sea surface height (neglecting lithospheric flexure).
Figure 1.

Schematic illustration of the physics of sea level change near a shoreline subject to dynamic topography and a change in sea surface height (neglecting lithospheric flexure).

In the first time step of the two-step evolution (t0 to t1 or frame a to b), a relatively constant perturbation to the air-loaded dynamic topography occurs producing an uplift of ΔDTA,1 (Fig. 1b). In the region landwards of the original shoreline, there is no amplification in the dynamic topography since the region begins and ends under air cover. As one moves from left to right between the shoreline position at t0 and t1, progressively more water is displaced by the dynamic topography and the unloading leads to a crustal uplift (or amplification from the air-loaded dynamic topography) that increases linearly (⁠|$\Delta R_1^L$|⁠) up to the new shoreline because at this point the crustal uplift has raised the bedrock to the sea-surface. To the right of this shoreline, the water unloading associated with the air-loaded dynamic topography is constant since the amount of water that is displaced is the same everywhere right of the shoreline at t1, and hence the response to this changing load is the same.

In the second time step (t1 to t2 or frame b to c), the sea surface height drops to the level G2, and the unloading produces crustal uplift at all locations to the right of the shoreline position at t1 (Fig. 1c). The uplift of the crust increases linearly as one moves from this shoreline to the location where the newly uplifted crust intersects the new sea surface height (⁠|$\Delta R_2^L$|⁠), and this intersection defines the shoreline at the end of the second time step. To the right of the shoreline at t2 the crustal uplift due to unloading is constant. The final position of the solid surface is given by the line R2.

The two step processes merges into one step in the evolution from Fig. 1(a) directly to Fig. 1(d). The drop in the sea surface from G0 to G1 and the contemporaneous perturbation in air-loaded dynamic topography, DTA,1, leads to a crustal uplift from R0 to R1 and a migration of the shoreline from t0 in Fig. 1(a) to its position at t1 in Fig. 1(d). To the left of the original shoreline the crust remains subaerial throughout the evolution and there is no ocean unloading effect. To the right of the new shoreline position, the crust remains under ocean cover and the unloading-induced uplift, |$\Delta R_1^L$|⁠, is constant. Between the old and new shoreline this crustal uplift increases monotonically towards the latter.

This simple example illustrates the equivalence between the water unloading and the amplification of air-loaded dynamic topography as discussed in the context of eq. (1). We emphasize, however, that as mentioned above, these schematics miss an important component of the sea level response, namely elastic flexure. If one includes an elastic crust/lithosphere layer in the physical response, as we do in the calculations described below, then the crustal uplift will extend laterally beyond the region of unloading, and this will impact both the amplitude of the sea level change and the predicted migration of the shoreline.

2.4 Including plate motion

We have thus far ignored the signal associated with the horizontal motion of tectonic plates, which would alter the geometry of the surface mass load and move any given geographic site through the field of mantle convection induced dynamic topography and geoid fluctuations. On timescales of 103–105 yr, this effect can generally be neglected. However, on longer timescales, the impact of plate motions needs to be considered.

Let us assume, in the following, that the Earth remains in a state of isostatic equilibrium in response to any change in the surface mass load associated with plate motions. In this case, the impact of changes in ice and sediment loads, dynamic topography and the gravity field on sea level can, for a given site, be most simply expressed in a reference frame fixed to the site (i.e. a Lagrangian frame of reference) rather than in a fixed mantle reference frame (e.g. hot spot reference frame).

As described in eq. (7), a perturbation in a general time-dependent quantity χ from the onset of loading, at time t0, to a time tj is simply the difference between the initial value of the quantity and its value at time j. If, across this time interval, a specific point on a given plate travels from (θ0, φ0) to (θj, φj) as a consequence of plate motion, then the response of the Earth at this site only depends on the net difference in the load and dynamic topography that the site experiences across this time interval, and not on the full time history from t0 to tj:
(26)

If we consider, for example, the air-loaded dynamic topography, then ΔDTA, j expressed in eq. (26) will change if the dynamic topography due to mantle convection evolves, or if the plate moves through a gradient in dynamic topography. The change in the parameter χ can be positioned in a coordinate system fixed to either the initial geography (i.e. Δχ = Δχ(θ0, φ0)) or the final geography (i.e. Δχ = Δχ(θj, φj)). In either case, there will be regions for which Δχ is undefined; areas of plate destruction (e.g. subduction zones) for the case Δχ(θ0, φ0), or areas of plate creation (e.g. mid-ocean ridges) for Δχ(θj, φj).

Within the sea level formalism, changes in parameters such as ice, sediments, the gravity field and dynamic topography must be defined globally to guarantee mass conservation. Therefore, any undefined regions have to be interpolated while preserving the conservation of water plus ice mass and solid Earth mass (or volume) between t0 and tj. This is a relatively easy procedure for changes in ice volume because ice sheets are generally located away from plate boundaries. For all other quantities, an appropriate interpolation scheme must be chosen.

3 APPLICATIONS

In the following section, we apply the gravitationally self-consistent sea level methodology derived above to a series of case studies that assume a spherically symmetric (i.e. depth varying) Earth that maintains isostatic equilibrium at all times (i.e. the case considered in the Appendix). Our goal is to highlight the accuracy of the new approach relative to previous, approximate treatments of sea level changes driven by dynamic topography. We use the density and elastic structure of the seismic model PREM (Dziewonski & Anderson 1981) in all calculations. Furthermore, while our standard results will adopt an elastic lithosphere of thickness 90 km, we will also consider a thickness of 130 km to investigate the sensitivity of the results to this parameter. While the calculations further assume a globally constant crustal density of 3300 kg/m3, the method developed in this paper does allow to include spatially varying crustal densities. Finally, to emphasize the underlying physics of sea level changes, all dynamic topography fields discussed below are simple, hypothetical geometries that are imposed rather than computed from a realistic simulation of mantle convection.

3.1 Example 1: isostatic amplification of dynamic topography

To illustrate the results of the new sea level methodology, we begin with a simple scenario in which an axisymmetric air-loaded dynamic topography perturbation, centred near Norfolk, West Virginia, is applied along the U.S. East Coast. The air-loaded perturbation is characterized by a peak uplift of 160 m. Figs 2(a)–(c) show the sea level change that one would predict by simply scaling the input air-loaded dynamic topography using different density contrast models, as discussed in Section 2.1. Specifically, we consider three cases in which we apply: (1) no amplification (air cover); (2) an amplification assuming a water cover; and (3) no amplification in zones of present-day land cover and water cover amplification in present-day ocean regions. Fig. 1(d) is the sea level change predicted using the full, gravitationally self-consistent treatment described in the last section. In each frame of the figure we superimpose the predicted shoreline prior to the perturbation in dynamic topography. Note that the shoreline after the perturbation is applied corresponds to the modern coastline and the shoreline prior to the perturbation is the ancient coastline. In practice, we calculate this by going backwards in time (as discussed in Section 2.2), starting at the known present-day topography and calculating the unknown ancient topography by applying the negative change in dynamic topography and in the case of Fig. 2(d), solving the sea level equation.

Change in sea level near the U.S. East Coast in response to a synthetic perturbation in dynamic topography. Frame (d) is the sea level change based on the gravitationally self-consistent sea level methodology described in the text, while frames (a)–(c) are based on the following assumptions: (1) air load everywhere, (2) water load everywhere and (3) air load over continents and water load over oceans. The white lines are the predicted shorelines prior to the change in dynamic topography.
Figure 2.

Change in sea level near the U.S. East Coast in response to a synthetic perturbation in dynamic topography. Frame (d) is the sea level change based on the gravitationally self-consistent sea level methodology described in the text, while frames (a)–(c) are based on the following assumptions: (1) air load everywhere, (2) water load everywhere and (3) air load over continents and water load over oceans. The white lines are the predicted shorelines prior to the change in dynamic topography.

To explore the error in the three approximate solutions, Figs 3(a)–(c) show differences between the prediction in Fig. 2(d) and each of the fields in Figs 2(a)–(c). The second column on Fig. 3 shows the same fields plotted in the first column, except that we zoom in on the U.S. East Coast. In all frames of Fig. 3, the white line is reproduced from the ancient shoreline position predicted from the gravitationally self-consistent calculation (Fig. 2d). The reason for this choice is that the frames of Fig. 3 represent errors in the prediction incurred by using an approximate treatment of sea level change, and our intent is to demonstrate how these errors will map onto the elevation of the “true”location of the ancient shoreline.

(a–c) Difference between the predicted sea level change calculated using a fixed amplification of an imposed, air-loaded dynamic topography (Figs 2a–c) and the sea level change calculated using our gravitationally self-consistent sea level methodology (Fig. 2d). (d–f) Close up of the frames (a–c), respectively.
Figure 3.

(a–c) Difference between the predicted sea level change calculated using a fixed amplification of an imposed, air-loaded dynamic topography (Figs 2a–c) and the sea level change calculated using our gravitationally self-consistent sea level methodology (Fig. 2d). (d–f) Close up of the frames (a–c), respectively.

As illustrated in the schematic of Fig. 1(b), the gravitationally self-consistent simulation accurately accounts for water mass unloading driven by dynamic topography. This unloading tends to increase as one moves from the ancient to modern shoreline and it remains relatively constant over the present-day ocean. The unloading associated with the simulation in Fig. 2(d) is more complex than in Fig. 1(b) because the initial bedrock topography (R0) is more variable, and because the former incorporates elastic flexure.

In regard to the approximate calculations of sea level change, Fig. 2(a) does not incorporate any water unloading and therefore the air-loaded dynamic topography is not amplified in this case. The error in the calculation of the sea level change (Figs 3a and d) is largely limited to regions oceanwards of the ancient shoreline (i.e. |$\Delta R_1^L$| in Fig. 1b is not captured). This is not strictly true because elastic flexure extends the effect of the water unloading landwards of the ancient shoreline in the accurate sea level calculation (Fig. 2d). In contrast to this case, Fig. 2(b) assumes that water unloading occurs everywhere. This leads to a significant error landwards of the ancient shoreline and an error that decreases as one moves from the ancient to modern shoreline (Figs 3b and e). The error in Figs 3(b) and (d) remains non-zero over a localized region oceanwards of the modern shoreline because flexural effects associated with water unloading are not accounted for. Finally, the calculation that assumes a water-amplified dynamic topography oceanwards of the modern shoreline, and no amplification landwards of this shoreline (Fig. 2c), neglects the effects of ocean unloading, and associated flexure, within the region of shoreline migration. As a consequence, the error is largest in the zone between the two shorelines, and it tapers to zero on either side of this region (Figs 3c and f). The geometries of the errors associated with the approximate treatments is dependent on whether the region is subject to a regression or transgression, but these plots nevertheless provide a sense of the scale and amplitude of these errors.

Fig. 4(a) shows the sea level change predicted for sites along the ancient shoreline in the gravitationally self-consistent calculation (blue; see also Fig. 2d). For the purpose of comparison, we also show predictions along the same profile for the three approximate treatments of sea level change in Figs 2(a)–(c). (The results for the case where no isostatic amplification is applied, Fig. 2(a), and the hybrid case treated in Fig. 2(c), are indistinguishable on Fig. 4(a) and therefore simply plot on top of each other.) All calculations predict a significant drop in sea level across the profile, as expected given the dynamic uplift of the region in this simulation (see Fig. 2). However, the air-loaded and hybrid calculations predict a sea level fall ∼15 per cent less than the gravitationally self-consistent calculation, while the water load approximation yields an amplitude ∼25 per cent greater than the calculation based on our new methodology.

Change in sea level calculated using an air load approximation (Fig. 2a), a water load approximation (Fig. 2b), a hybrid approximation (Fig. 2c) and the full gravitationally self-consistent methodology (Fig. 2d) with a lithospheric thickness of 90 km (solid blue) and 130 km (dashed blue). (a) Prediction along the location of the ancient shoreline (see the inset) as calculated with the full methodology. The air load approximation and hybrid approximation are identical over the modern continents and therefore these two approximations yield the same values. (b) Prediction perpendicular to the shoreline (see the inset). The hybrid approximation changes from the air load approximation to the water load approximation at the position of the modern shoreline. A thicker lithosphere leads to more pronounced flexure that extends further into the continent and oceans.
Figure 4.

Change in sea level calculated using an air load approximation (Fig. 2a), a water load approximation (Fig. 2b), a hybrid approximation (Fig. 2c) and the full gravitationally self-consistent methodology (Fig. 2d) with a lithospheric thickness of 90 km (solid blue) and 130 km (dashed blue). (a) Prediction along the location of the ancient shoreline (see the inset) as calculated with the full methodology. The air load approximation and hybrid approximation are identical over the modern continents and therefore these two approximations yield the same values. (b) Prediction perpendicular to the shoreline (see the inset). The hybrid approximation changes from the air load approximation to the water load approximation at the position of the modern shoreline. A thicker lithosphere leads to more pronounced flexure that extends further into the continent and oceans.

The sign of the error associated with the air load and hybrid calculations can be understood with reference to Fig. 1(b), which illustrates the cross-sectional profile expected from a gravitationally self-consistent solution of sea level change. In contrast to this case, the approximate solutions neglect the water unloading in the region between the ancient and modern shoreline (and the flexural effects of this unloading) because they assume air loading for all sites landwards of the modern shoreline (corresponding to the shoreline at t1). Thus, within the region of shoreline migration the approximate solution does not accurately capture the uplift of the crust and thus it predicts a sea level fall of smaller amplitude than the gravitationally self-consistent solution in the vicinity of this region. In contrast, the water load approximation overestimates the change in load (and hence crustal uplift) in all regions landwards of the modern shoreline and therefore predicts a sea level fall of greater amplitude than the self-consistent calculation.

Fig. 4(b) shows the computed sea level change perpendicular to the shoreline. The water load calculation overestimates the sea level change over continents while the air load calculation underestimates the change over oceans. The self-consistent calculation is characterized by a smooth transition between the two regions across the zone of shoreline migration and slightly beyond. The modern shoreline is marked by the transition of the hybrid approximation from the air to water load prediction. Assuming a thicker lithosphere (dashed blue line) results in broader flexural effects and a wider zone of transition.

3.2 Example 2: can independent signals be summed?

The gravitationally self-consistent sea level methodology described above can be used to consider the impact of dynamic topography in isolation (as in the example of the last section), however, it also admits the possibility of contemporaneous dynamic topography fluctuations, ice mass changes and sediment fluxes. To date, studies have computed these effects independently and summed them up (e.g. Müller et al.2008; Rowley et al.2013). This assumes the processes are independent of one another, an assumption that does not necessarily hold. To consider the potential error introduced in the latter approach, we return to, and extend, the numerical example treated in the last section. In particular, we apply the same change in air-loaded dynamic topography as shown in Fig. 2(a) but in addition we assume a contemporaneous growth of an ice sheet in the far field of the U.S. East Coast that leads to a fall in sea level over the oceans of approximately 30m. These tandem effects are analogous to the example illustrated schematically in Fig. 1.

To begin, we perform a single, gravitationally self-consistent calculation that considers both processes, a change in dynamic topography and an increase in ice volume, together. The change in sea level computed for this case using our new sea level methodology is given in Fig. 5(a). Next, we perform two independent calculations. The first involves the perturbation in dynamic topography and no far-field change in ice volume, and the second includes the change in ice volume in the absence of a signal from dynamic topography. Both calculations are also based on our gravitationally self-consistent sea level methodology. We then sum the sea level changes obtained from these two calculations. Fig. 5(b) shows the difference between the latter ‘superposition’ procedure and the full calculation shown in Fig. 5(a).

(a) Change in sea level near the U.S. East Coast in response to a perturbation in both dynamic topography and ice volume in the far field of the region, as computed using the gravitationally self-consistent sea level methodology described in the text. The white line shows the predicted location of the ancient shoreline. (b) Difference between the calculation shown in frame (a) and an approximate approach in which the sea level change due to both perturbations is computed separately using the same sea level method and the total response is calculated by summing the two effects.
Figure 5.

(a) Change in sea level near the U.S. East Coast in response to a perturbation in both dynamic topography and ice volume in the far field of the region, as computed using the gravitationally self-consistent sea level methodology described in the text. The white line shows the predicted location of the ancient shoreline. (b) Difference between the calculation shown in frame (a) and an approximate approach in which the sea level change due to both perturbations is computed separately using the same sea level method and the total response is calculated by summing the two effects.

The error incurred in the approximate solution is greatest in areas of shoreline migration, but it extends beyond this region as a consequence of flexural effects. We can understand this trend by once again considering Fig. 1. The full solution accurately accounts for the total migration of the shoreline as a consequence of the tandem effects of a dynamic topography change and far-field drawdown of the sea surface (i.e. Figs 1a–d). However, the approximate solution yields a different initial shoreline for each process. In the example of Fig. 5, both processes act to move the shoreline oceanwards, and thus the approximate solution errs by not accounting for the full lateral extent of the shoreline migration (i.e. the migration of the shoreline at t0 in Fig. 1a to the position at t1 in Fig. 1d) and thus underestimates the water loading. The sea level change computed using the approximate treatment will predict a sea level fall of lower magnitude than the full solution, and this leads to the sign of the error in Fig. 5(b).

The amplitude of the error in the sea level prediction along the predicted location of the ancient shoreline is about 10 per cent of the sea level fall over the oceans due to far field ice growth (30 m). This is the level of error one would incur if the geological record were corrected for dynamic topography (using a gravitationally self-consistent calculation) and the residual used to infer the ice volume change. The error would of course be much larger if one had used the approximate treatments of sea level change due to dynamic topography (as illustrated in Fig. 4). Finally, differences in the location of the ancient shoreline, as predicted using the gravitationally self-consistent treatment of sea level change and the approximation in which the effects are computed independently, differ negligibly (compare the white lines in both frames of Fig. 5). We note, however, that this agreement largely reflects the relatively high topographic slope in the region of the ancient shoreline. A shallow slope would, of course, result in a greater discrepancy in the predicted shoreline migration, and this would in turn influence the amplitude of the error incurred in the approximate solution.

4 SUMMARY

In this study we have outlined a gravitationally self-consistent formalism for including dynamic topography and perturbations in the gravity field due to mantle convection in calculations of sea level change. Our formalism is based on ice age sea level theory (Mitrovica & Milne 2003; Kendall et al.2005; Dalca et al.2013) and indeed includes the possibility of ice mass changes and sediment redistribution. It further accurately accounts for shoreline migration, conservation of the mass of ice plus water and, of course, the mass of the solid Earth. (The latter is enforced by ensuring that the globally averaged radial deformation of the crust in response to the total surface mass load is zero.)

Our formulation is valid for arbitrary Earth models but we also consider the special case of a 1-D (depth varying) Earth model that maintains isostatic equilibrium throughout the simulation. In this case, the relevant response equations are based on fluid Love number theory. Numerical examples based on this special case have allowed us to explore the level of error incurred in previous predictions of sea level change that (1) assume a simplified amplification of air-loaded dynamic topography to account for the influence of water cover and (2) consider the impact of ice mass changes and dynamic topography on sea level in isolation. These demonstrate that the error associated with the simplified treatment of the amplification of dynamic topography can vary from ∼−15 per cent to ∼25 per cent of the predicted sea level change. Superimposing ice mass and dynamic topography effects introduces an additional error of 10 per cent in the ice volume estimate. Rowley et al. (2013) have argued for a dynamic topography change in water-loaded dynamic topography of up to 50 m along the U.S. East Coast since the mid-Pliocene climate optimum. In this case, the two sources of error we have identified might lead to an error of ∼10 m in any effort to use geological data from this region to infer sea level equivalent changes in ice mass flux over the last 3 Myr. This level of error is significant and comparable to the total eustatic sea level change associated with the collapse of the Greenland and West Antarctic Ice Sheets.

As interest in ice age sea level change and ice sheet stability broadens from the usual focus on Holocene and modern sea level to encompass records of early Pleistocene and Pliocene age, traditional glacial isostatic adjustment analyses must account for the sea level signal associated with mantle dynamics. The methodology presented here is an attempt to incorporate the latter signal, and in particular dynamic topography and associated geoid fluctuations, into the gravitationally self-consistent theory of ice age sea level change that dates to the pioneering work of Farrell & Clark (1976).

We thank two anonymous reviewers for their constructive comments on this manuscript, as well as Alessandro Forte, Rob Moucha and Dave Rowley for many discussions related to the ideas in this paper. We acknowledge funding from Harvard University and the Harvard University Center for the Environment.

REFERENCES

Chan
N.-H.
Mitrovica
J.
Matsuyama
I.
Creveling
J.
Stanley
S.
2011
The rotational stability of a convecting earth: assessing inferences of rapid TPW in the late cretaceous
Geophys. J. Int.
187
3
1319
1333

Coakley
B.
Gurnis
M.
1995
Far-field tilting of Laurentia during the Ordovician continent
J. geophys. Res.
100
6313
6327

Conrad
C.P.
Husson
L.
2009
Influence of dynamic topography on sea level and its rate of change
Lithosphere
1
2
110
120

Crameri
F.
et al.
2012
A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the ‘sticky air’ method
Geophys. J. Int.
189
1
38
54

Dahlen
F.A.
1976
The passive influence of the oceans upon the rotation of the Earth
Geophys. J. R. astr. Soc.
46
363
406

Dalca
A.V.
Ferrier
K.L.
Mitrovica
J.X.
Perron
J.T.
Milne
G.A.
Creveling
J.R.
2013
On postglacial sea level-III. Incorporating sediment redistribution
Geophys. J. Int.
194
1
45
60

Daradich
A.
Mitrovica
J.X.
Pysklywec
R.N.
Willett
S.D.
Forte
A.M.
2003
Mantle flow, dynamic topography, and rift-flank uplift of Arabia
Geology
31
10
901
904

Dziewonski
A.
Anderson
D.
1981
Preliminary reference Earth model (PREM)
Phys. Earth planet. Inter.
25
297
356

Farrell
W.E.
Clark
J.A.
1976
On postglacial sea level
Geophys. J. R. astr. Soc.
46
3
647
667

Flament
N.
Gurnis
M.
Müller
R.D.
2013
A review of observations and models of dynamic topography
Lithosphere
5
2
189
210

Forte
A.M.
Peltier
W.R.
Dziewonski
A.M.
Woodward
R.L.
1993
Dynamic surface topography: a new interpretation based upon mantle flow models derived from seismic tomography
Geophys. Res. Lett.
20
3
225
228

Forte
A.M.
Quéré
S.
Moucha
R.
Simmons
N.A.
Grand
S.P.
Mitrovica
J.X.
Rowley
D.B.
2010
Joint seismic–geodynamic-mineral physical modelling of African geodynamics: a reconciliation of deep-mantle convection with surface geophysical constraints
Earth planet. Sci. Lett.
295
3–4
329
341

Golle
O.
Dumoulin
C.
Choblet
G.
Čadek
O.
2012
Topography and geoid induced by a convecting mantle beneath an elastic lithosphere
Geophys. J. Int.
189
1
55
72

Gurnis
M.
1990
Bounds on global dynamic topography from Phanerozoic flooding of continental platforms
Nature
344
754
756

Gurnis
M.
1992
Rapid continental subsidence following the initiation and evolution of subduction
Science
255
1556
1558

Gurnis
M.
1993
Phanerozoic marine inundation of continents driven by dynamic topography above subducting slabs
Nature
364
589
593

Gurnis
M.
Eloy
C.
Zhong
S.
1996
Free-surface formulation of mantle convection-II. Implication for subduction-zone observables
Geophys. J. Int.
127
719
721

Gurnis
M.
Müller
R.D.
Moresi
L.
1998
Cretaceous vertical motion of Australia and the Australian–Antarctic discordance
Science
279
5356
1499
1504

Gurnis
M.
Mitrovica
J.X.
Ritsema
J.
van Heijst
H.-J.
2000
Constraining mantle density structure using geological evidence of surface uplift rates: the case of the African superplume
Geochem. Geophys. Geosyst.
1
7
doi:10.1029/1999GC000035

Hager
B.H.
Clayton
R.W.
Richards
M.A.
Comer
R.P.
Dziewonski
A.M.
1985
Lower mantle heterogeneity, dynamic topography and the geoid
Nature
313
541
545

Han
D.
Wahr
J.
Cohen
S. C.
VanííEk
P.
1989
Post-glacial rebound analysis for a rotating earth
Slow Deformation and Transmission of Stress in the Earth
American Geophysical Union
1
6

Holt
W.E.
Stern
T.A.
1994
Subduction, platform subsidence, and foreland thrust loading: the late tertiary development of Taranaki Basin, New Zealand
Tectonics
13
5
1068
1092

Kaus
B.J.
Mühlhaus
H.
May
D.A.
2010
A stabilization algorithm for geodynamic numerical simulations with a free surface
Phys. Earth planet. Inter.
181
1-2
12
20

Kendall
R.A.
Mitrovica
J.X.
Milne
G.A.
2005
On post-glacial sea level-II. Numerical formulation and comparative results on spherically symmetric models
Geophys. J. Int.
161
3
679
706

Kramer
S.C.
Wilson
C.R.
Davies
D.R.
2012
An implicit free surface algorithm for geodynamical simulations
Phys. Earth planet. Inter.
194–195
25
37

Latychev
K.
Mitrovica
J.X.
Tromp
J.
Tamisiea
M.E.
Komatitsch
D.
Christara
C.C.
2005
Glacial isostatic adjustment on 3-D Earth models: a finite-volume formulation
Geophys. J. Int.
161
2
421
444

Lithgow-Bertelloni
C.
Gurnis
M.
1997
Cenozoic subsidence and uplift of continents from time-varying dynamic topography
Geology
25
8
735
738

Lithgow-Bertelloni
C.
Silver
P.G.
1998
Dynamic topography, plate driving forces and the African superswell
Nature
395
345
348

Liu
S.
Nummedal
D.
Gurnis
M.
2014
Dynamic versus flexural controls of Late Cretaceous Western Interior Basin, USA
Earth planet. Sci. Lett.
389
221
229

McKenzie
D.
1977
Surface deformation, gravity anomalies and convection
Geophys. J. R. astr. Soc.
48
211
238

Milne
G.A.
Mitrovica
J.X.
1996
Postglacial sea-level change on a rotating Earth: first results from a gravitationally self-consistent sea-level equation
Geophys. J. Int.
126
3
F13
F20

Milne
G.A.
Mitrovica
J.X.
1998
Postglacial sea-level change on a rotating Earth
Geophys. J. Int.
133
1
1
19

Mitrovica
J.X.
Milne
G.A.
2003
On post-glacial sea level: I. General theory
Geophys. J. Int.
154
2
253
267

Mitrovica
J.X.
Peltier
W.R.
1991
On postglacial geoid subsidence over the equatorial oceans
J. geophys. Res.
96
91
20 053
20 071

Mitrovica
J.X.
Beaumont
C.
Jarvis
G.T.
1989
Tilting of continental interiors by the dynamical effects of subduction
Tectonics
8
5
1079
1094

Mitrovica
J.X.
Pysklywec
R.N.
Beaumont
C.
Rutty
A.
1996
The Devonian to Permian sedimentation of the Russian platform: an example of subduction controlled long-wavelength tilting of continents
J. Geodyn.
22
112
79
96

Mitrovica
J.X.
Wahr
J.
Matsuyama
I.
Paulson
A.
2005
The rotational stability of an ice-age earth
Geophys. J. Int.
161
2
491
506

Moucha
R.
Forte
A.M.
Mitrovica
J.X.
Rowley
D.B.
Quéré
S.
Simmons
N.A.
Grand
S.P.
2008
Dynamic topography and long-term sea-level variations: there is no such thing as a stable continental platform
Earth planet. Sci. Lett.
271
1-4
101
108

Moucha
R.
Forte
A.M.
Rowley
D.B.
Mitrovica
J.X.
Simmons
N.A.
Grand
S.P.
2009
Deep mantle forces and the uplift of the Colorado Plateau
Geophys. Res. Lett.
36
19
L19310

Mound
J.E.
Mitrovica
J.X.
1998
True polar wander as a mechanism for second-order sea-level variations
Science
279
5350
534
537

Müller
R.D.
Sdrolias
M.
Gaina
C.
Steinberger
B.
Heine
C.
2008
Long-term sea-level fluctuations driven by ocean basin dynamics
Science
319
5868
1357
1362

Peltier
W.
1974
The impulse response of a Maxwell Earth
Rev. Geophys. Space Phys.
12
4
649
669

Peltier
W.
1976
Glacial-isostatic adjustment-II. The inverse problem
Geophys. J. R. astr. Soc.
46
669
705

Piromallo
C.
Spada
G.
Sabadini
R.
Ricard
Y.
1997
Sea-level fluctuations due to subduction: the role of mantle rheology
Geophys. Res. Lett.
24
13
1587
1590

Pysklywec
R.N.
Mitrovica
J.X.
1998
Mantle flow mechanisms for the large-scale subsidence of continental interiors
Geology
26
8
687
690

Pysklywec
R.N.
Mitrovica
J.X.
1999
The role of subduction-induced subsidence in the evolution of the Karoo Basin
J. Geol.
107
2
155
164

Ricard
Y.
Wuming
B.
1991
Inferring the viscosity and the 3-D density structure of the mantle from geoid, topography and plate velocities
Geophys. J. Int.
105
3
561
571

Ricard
Y.
Spada
G.
Sabadini
R.
1993
Polar wandering of a dynamic earth
Geophys. J. Int.
113
2
284
298

Rowley
D.B.
Forte
A.M.
Moucha
R.
Mitrovica
J.X.
Simmons
N.A.
Grand
S.P.
2013
Dynamic topography change of the eastern United States since 3 million years ago
Science
340
1560
1563

Schmeling
H.
et al.
2008
A benchmark comparison of spontaneous subduction models—towards a free surface
Phys. Earth planet. Inter.
171
1–4
198
223

Spasojevic
S.
Gurnis
M.
2012
Sea level and vertical motion of continents from dynamic earth models since the late cretaceous
AAPG Bull.
96
11
2037
2064

Steinberger
B.
2007
Effects of latent heat release at phase boundaries on flow in the Earth's mantle, phase boundary topography and dynamic topography at the Earth's surface
Phys. Earth planet. Inter.
164
1–2
2
20

Steinberger
B.
O'Connell
R.W.
Mitrovica
J.X.
Vermeersen
L.
The convective mantle flow signal in rates of true polar wander
Ice Sheets, Sea Level and the Dynamic Earth, Geodynamics Series, Vol. 29
2002
233
256
AGU

Tromp
J.
Mitrovica
J.X.
1999
Surface loading of a viscoelastic earth - I. General theory
Geophys. J. Int.
137
847
855

Wen
L.
Anderson
D.L.
1997
Layered mantle convection: A model for geoid and topography
Earth planet. Sci. Lett.
146
3–4
367
377

Wu
P.
Peltier
W.R.
1982
Viscous gravitational relaxation
Geophys. J. Int.
70
435
485

Wu
P.
van der Wal
W.
2003
Postglacial sealevels on a spherical, self-gravitating viscoelastic earth: effects of lateral viscosity variations in the upper mantle on the inference of viscosity contrasts in the lower mantle
Earth planet. Sci. Lett.
211
1–2
57
68

Zhong
S.
Gurnis
M.
Hulbert
G.
1993
Accurate determination of surface normal stress in viscous flow from a consistent boundary flux method
Phys. Earth planet. Inter.
78
1–2
1
8

Zhong
S.
Gurnis
M.
Moresi
L.
1996
Free-surface formulation of mantle convection-I. Basic theory and application to plumes
Geophys. J. Int.
127
708
718

Zhong
S.
Paulson
A.
Wahr
J.
2003
Three-dimensional finite-element modelling of Earth's viscoelastic deformation: effects of lateral variations in lithospheric thickness
Geophys. J. Int.
155
679
695

Zhong
S.
McNamara
A.
Tan
E.
Moresi
L.
Gurnis
M.
2008
A benchmark study on mantle convection in a 3-D spherical shell using CitcomS
Geochem. Geophys. Geosyst.
9
10
1
32

APPENDIX A: EXPRESSIONS FOR |$\Delta \mathcal {SL}$| IN THE CASE OF SPHERICALLY SYMMETRIC EARTH MODELS

This appendix derives the relevant equations, based on Love number theory, for a prediction of the change in global sea level due to the surface mass loading of a spherically symmetric Earth model. The derivation expands the theory described in Kendall et al. (2005) and Dalca et al. (2013) to include dynamic topography and prescribed perturbations in the geoid. However, we only treat the isostatic equilibrium response, that is, we assume that the viscous response has fully decayed. We begin with the case of a non-rotating Earth.

A1 Non-rotating Earth model

The deformational and gravitational response of a spherically symmetric Earth model to a surface mass load can be written as a space–time convolution of the appropriate Green's function with the surface load. In particular, the response Δχ can be expressed as:
(A1)
where GF denotes the Green's function, θ and φ represent the colatitude and east longitude, and γ is the angle between the observation point (θ, φ) and the load point (θ, φ). This angle can be calculated as:
(A2)
If we denote the Green's function for the radial displacement of the solid surface and the gravitational potential perturbation by Γ(γ, t) and ϕ(γ, t), respectively, then these responses can be expressed as:
(A3)
(A4)
where ΔL is given by eq. (12). Using these expressions in eq. (16) gives
(A5)
The Green's functions in the above equations can be written in terms of the viscoelastic k and h Love numbers at spherical harmonic degree l (Peltier 1974; Tromp & Mitrovica 1999):
(A6)
(A7)
where the response is decomposed into an instantaneous elastic contribution and a discrete set of viscous decays.
As we have noted, our calculations will assume that the viscous response has fully decayed. In this case, our expression for the perturbation in global sea level can be simplified by expressing the surface mass load as having a single Heaviside step function time dependence, and performing the time convolution between the step function and the viscoelastic Green's functions analytically (where the latter is expressed in terms of the viscoelastic Love numbers). At infinite time, this convolution yields the following time-independent Green's functions:
(A8)
(A9)
where a and Me are the radius and mass of the Earth, respectively, and Pl is the Legendre polynomial at spherical harmonic degree l. In these expressions, |$h_l^{f}$| and |$k_l^{f}$| are the so-called fluid Love numbers that have the form (Peltier 1974, 1976):
(A10)
(A11)
Using these equilibrium forms, and the expression (12) for the change in surface mass load, yields the following version of eq. (A5) for the change in global sea level:
(A12)
where
(A13)
Note that under the assumption of isostatic equilibrium, computing the perturbation in sea level between any two time steps requires only the total difference in the load across the interval and thus, with no loss of generality, we consider only a single time step and drop the subscript j.
We can simplify the spatial convolution in eq. (A12) by introducing spherical harmonic basis functions. A general spherical harmonic decomposition of a scalar field χ can be written as:
(A14)
where the basis functions are normalized in the following way
(A15)
In this equation, the asterisk denotes the complex conjugate. In this case one can show that the following relationship holds (Mitrovica & Peltier 1991)
(A16)
We may now use the spherical harmonic formulation from eqs (A14) and (A16) to rewrite eq. (A12):
(A17)
where the subscript l, m refers to a spherical harmonic coefficient of degree l and order m,
(A18)
and we adopt the short form:
(A19)
The fully spectral form of eq. (A17) is
(A20)
This equation describes the spatially variable component to sea level change. The spatially invariant component that follows from volume conservation (19) can also be written in terms of a spherical harmonic expansion:
(A21)
Here we have made use of the following property
(A22)
We also used the following projections
(A23)
and
(A24)
where C represents the final ocean function.

A2 Rotating Earth model

Changes in the Earth's surface mass load (ice, ocean, sediments) and convective motions in the interior perturb the inertia tensor of the Earth, and thus the magnitude and orientation of the Earth's rotation axis. A change in the rotation vector, and the perturbation in the centrifugal potential that this change implies, impacts global sea level (Han & Wahr 1989; Milne & Mitrovica 1996, 1998; Mound & Mitrovica 1998). In this section, we review how this rotational feedback into global sea level can be incorporated into the Love number theory outlined in Appendix A1.

Convective perturbations to Earth rotation can be computed from the same viscous flow models of mantle convection used to compute dynamic topography and geoid anomalies, and a variety of approximate methods have been derived for this purpose (Ricard et al.1993; Steinberger & O'Connell 2002; Chan et al.2011). In ice age research, considerations of rotational stability deepen the integral nature of the sea level calculation since perturbations in the rotation vector are impacted by sea level changes and act, in turn, to perturb sea level. This rotational feedback can be included by augmenting eq. (A20) to include a term associated with the sea level response to the perturbed centrifugal potential (Milne & Mitrovica 1998). In the present case, this additional term will include a signal associated with convection, which is prescribed at the outset, and a surface mass loading term that will be successively improved through the iteration procedure described above. The latter requires a method for predicting the perturbation to the rotation vector associated with a surface mass loading (Mitrovica et al.2005). To be consistent with the assumption applied in Appendix A1, the ice age component of this calculation would treat the case of isostatic equilibrium.

The time-varying centrifugal potential Λ(θ, φ, t) can be decomposed into an initial value prior to the onset of loading and a perturbation from this state due to loading:
(A25)
The perturbation to the centrifugal potential, also known as the rotational driving potential, may be described as the difference between two ellipsoidal forms with distinct amplitude and orientation, and it will involve only spherical harmonic coefficients of degree zero and two. This term can therefore be written as (e.g. Milne & Mitrovica 1996):
(A26)
In analogy to the case of surface mass loading, fluid Love numbers exist that describe the deformation of the solid surface and the gravitational equipotential surface due to a changing rotational driving potential. These have the form (Peltier 1974, 1976):
(A27)
(A28)
Using these expressions, eq. (A20) is revised in the following manner to include rotational feedback:
(A29)
where
(A30)
All other equations described in Appendix A1 remain unchanged.