the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A kinematic formalism for tracking ice–ocean mass exchange on the Earth's surface and estimating sealevel change
Surendra Adhikari
Erik R. Ivins
Eric Larour
Lambert Caron
Helene Seroussi
Polar ice sheets are important components of the Earth system. As the geometries of land, ocean and ice sheets evolve, they must be consistently captured within the lexicon of geodesy. Understanding the interplay between the processes such as icesheet dynamics, solidEarth deformation, and sealevel adjustment requires both geodetically consistent and massconserving descriptions of evolving land and ocean domains, grounded ice sheets and floating ice shelves, and their respective interfaces. Here we present mathematical descriptions of a generic level set that can be used to track both the grounding lines and coastlines, in light of ice–ocean mass exchange and complex feedbacks from the solid Earth and sea level. We next present a unified method to accurately compute the sealevel contribution of evolving ice sheets based on the change in ice thickness, bedrock elevation and mean sea level caused by any geophysical processes. Our formalism can be applied to arbitrary geometries and at all timescales. While it can be used for applications with modeling, observations and the combination of two, it is best suited for Earth system models, comprising ice sheets, solid Earth and sea level, that seek to conserve mass.
© 2020 California Institute of Technology. Government sponsorship acknowledged.
Recently there has been intense interest in defining the physics involved in determining multidecadal change in the location and the migration rate of the grounding line, a boundary separating a grounded ice sheet from its floating extension, usually a floating ice shelf (e.g., Nowicki and Wingham, 2008; Schoof, 2012; Sergienko and Wingham, 2019). Indeed, how well a numerical model of marine ice sheets predicts the sealevel contribution largely depends on its ability to capture the subtle migration of grounding lines. The nonequilibrium thermodynamics, fluid dynamics, plastic failure criteria and conditions governing nonlinear stability of ice sheets are, quite generally, up for lively debate. In order to better tackle the difficult nonlinear physics and to better address the associated numerical challenges (e.g., Schoof, 2007; Durand et al., 2009; Sayag and Grae Worster, 2013; Favier et al., 2016; Seroussi and Morlighem, 2018) as well as to define proper observational criteria for locating the grounding lines and their migrations (e.g., Hogg et al., 2018; Milillo et al., 2019), it is important to agree on some of the baseline variables and boundary conditions. Direct interactions with the ocean (e.g., Seroussi et al., 2017; Nakayama et al., 2018) and the solid Earth (e.g., Gomez et al., 2010; Larour et al., 2019) are now seen as critical elements that must be incorporated into projections, or retrospective paleoclimate simulations, of the rate of grounding line retreat in a warming climate (e.g., Jones et al., 2015; Whitehouse et al., 2017). Given the computational complexity of this problem, however, it is essential to properly define the simple geometrical parameters, primarily moving boundaries at the ice–bedrock–ocean interfaces, for there to be rationally organized intercomparison among various research teams and their results.
A general description of mechanical analysis of icesheet evolution at the ice–bedrock–ocean interfaces has been given for a set of simplified geometries (for example in Chapter 3 of Hutter, 1983), owing to the lack of constraining data or computational resources. A similar geometric approach is also familiar in the development of glacial isostatic adjustment (GIA) theory for ice sheets, sea level and bedrock evolution following the Last Glacial Maximum with migrating grounding lines and coastlines (e.g., Milne, 1998; Lambeck et al., 2003; Mitrovica and Milne, 2003). Modern satellite techniques have allowed us to gain knowledge of both the present locations and migration rates of the grounding line (e.g., Rignot et al., 2011; Milillo et al., 2019). However, both observations and numerical simulations of subtle change in grounding line positions are complicated by the presence of kilometerscale geometric features, such as ice rises and rumples, rugged fjord geometries, and uneven bedrock topography. These features complicate the required geometrical simplifications used in the previous studies of ice–bedrock–ocean interface changes, especially when the system of ice sheets, solid Earth and sea level is fully interactive (e.g., Lingle and Clark, 1985; Gomez et al., 2013; de Boer et al., 2014; Konrad et al., 2015; Larour et al., 2019). Here we consider a simple levelset method, which has been previously applied for tracking grounding lines (e.g., Seroussi et al., 2014) and calvingfront positions (e.g., Bondzio et al., 2017), and generalize it to facilitate a precise tracking of both the grounding lines and coastlines of arbitrary geometries in a seamless manner. The method is very generic and can be used for applications based on modeling, observations, or combination of models and observations.
Evolving bedrock and sea level impact the icesheet dynamics via the modulation of bedrock slope, groundingline positions and gravitational driving stress. For the marine portions of the ice sheet having retrograde bedrock slopes, this effect promotes the stability, as has been demonstrated by both the observationbased (e.g., Barletta et al., 2018; Kingslake et al., 2018) and the modelbased studies (e.g., Lingle and Clark, 1985; Gomez et al., 2010, 2015; Adhikari et al., 2014; Konrad et al., 2015; Larour et al., 2019). The inclusion of evolving bedrock and sea level in a dynamical icesheet model, however, requires a modification to the common method of estimating sealevel contribution. The method, based on the concept of ice height above floatation (e.g., Bindschadler et al., 2013), yields inaccurate results for the marine portions of the ice sheet. Goelzer et al. (2020) recently provide appropriate corrections for the effects of bedrock elevation change and externally forced sea level. Our goal here is to formulate a unified method to calculate the exact fraction of ice thickness that contributes to the sealevel change over a given period by considering evolving bedrock and sea level driven by any geophysical processes. In conjunction with observational data, the method can be applied to a variety of models such as standalone icesheet models and those that account for isostatic bedrock adjustment (e.g., Le Meur and Huybrechts, 1996; Pattyn, 2017) or a selfconsistent GRD (gravitational, rotational, deformational) response of the solid Earth (e.g., Gomez et al., 2013; Larour et al., 2019). In the latter set of models, the presented formalism ensures mass conservation in the Earth system by exchanging mass between the land and the ocean, accounting for the induced GRD response of solid Earth, and adjusting the ocean area through migration of grounding lines and coastlines, simultaneously.
In the following, we begin by presenting a generalized description of land, ocean and ice domains and their respective interfaces (Sect. 2). We consider a global ocean, composed of an interconnected system of oceanic basins, and distributed system of ice domains, comprising glaciers, ice sheets and ice shelves, that can be straightforwardly employed in any Earth system model in order to track the global mass transport and assess the evolution of a dynamic system of ice sheets, solid Earth and sea level. In Sect. 3, we briefly review the common method of estimating sealevel contribution of ice sheets and present a new method, wherein we isolate mass and volume contributions to the ocean, which is critical to accurately drive the GRD response of the solid Earth. In Sect. 4, we assess our formalism in a broader context of sealevel change and mass conservation in the Earth system. Finally, in Sect. 5, we summarize the key conclusions.
To begin our discussion, we consider a spherical planet whose surface is divided into complementary domains of land and ocean. The ocean may be thought of as an interconnected system of oceanic basins – just like Earth's ocean that also includes fjords and marginal seas such as the Mediterranean – that are able to freely exchange and redistribute mass between them. This assumption simplifies what would otherwise be an arduous task for mass attribution and conservation in the Earth system. Distributed ice domains including glaciers, ice sheets and ice shelves exist on the land or float on the ocean (Fig. 1). In order to present mathematical descriptions of these domains and their interfaces at time t, we denote 2D spatial coordinates on the planetary surface by ω. Depending upon the spatial scale (e.g., the ocean vs. glaciers), we interchangeably use ω to represent geographic coordinates (θ,ϕ) or Cartesian (x,y), assuming that an appropriate coordinate transformation is applied. The entire formalism presented in this study can be derived from three field variables: the solid Earth surface (i.e., land surface or sea floor) or simply bedrock B(ω,t), mean sea level (MSL) S(ω,t) and ice thickness H(ω,t). The first two fields must be defined relative to the same reference ellipsoid (e.g., Altamimi et al., 2016).
Our definition of MSL complies with that given by Gregory et al. (2019): time mean of sea surface over a sufficiently long period so that the effects of waves, tides or meteorologically driven highfrequency fluctuations are eliminated. The period of time mean may be on the order of 20 years or longer, a timescale over which interactions between sea level and ice sheet may become important (e.g., Hillenbrand et al., 2017; Larour et al., 2019). One key difference, however, is that the change in MSL in the present context does not account for the steric component that is due to the change in the ocean density. Here, in the strict sense of the word, the change in global mean of MSL is given by the socalled barystatic sealevel change, which is the globalmean sealevel (GMSL) change due to the exchange of water between the land and the ocean, and the evolving spatial pattern of MSL is dictated by the GRD response of the solid Earth to land–ocean (water or sediment) mass exchange and the tectonic activities. This definition of evolving MSL is familiar in GIA modeling wherein there is a requirement to solve for a gravitationally selfconsistent solution of evolving bedrock and (nonsteric) MSL driven by the ice–ocean mass exchange following the Last Glacial Maximum (Farrell and Clark, 1976; Milne and Mitrovica, 1998). The MSL as defined above represents an equipotential surface whose spatial pattern matches the geoid (Tamisiea, 2011).
2.1 Coastlines and grounding lines as a seamless interface
We develop our formalism based on the principle of hydrostatic equilibrium for a system of ice and ocean. Since H(ω,t) may be considered a globally defined field, with $H(\mathit{\omega},t)=\mathrm{0}$ outside the ice domains, this concept can be generalized to deduce a criterion for delineating boundaries between the land and the ocean, as well as the floating and the grounded ice. We define
such that $F(\mathit{\omega},t)=\mathrm{0}$ satisfies the hydrostatic equilibrium between ice and the oceanic water in the marine sectors where $B(\mathit{\omega},t)<S(\mathit{\omega},t)$. Here ρ_{i} and ρ_{o} are the average densities of ice and ocean water, respectively. Our goal here is to use Eq. (1) as a basis for defining the land–ocean boundaries consistently. The equation, at first glance, suggests that the ocean (land) takes negative (positive) values of F(ω,t) and their interfaces have zero values. However, a few aspects should be further clarified. To simplify a mathematical description of the land–ocean boundaries, we ensure the absence of marine ice cliffs that have larger thickness than the floatation height (i.e., negative of the second term on the right side of the equation) by assuming that $F(\mathit{\omega},t)\le \mathrm{0}$ at the “ice front”. The ice front that satisfies the equality (inequality) here represents the calving face of a tidewater glacier (an ice shelf). Along the same lines, we assume that the terrestrial ice cliffs are not present where $B(\mathit{\omega},t)=S(\mathit{\omega},t)$. These assumptions are generally valid, as the ice flow is diffusive on the timescale (decades or longer) we are interested in.
We now define a level set of the function F(ω,t) such that
This zerolevel set may consist of several simple curves, T_{i}(F), that divide the planetary surface into several nonoverlapping regions, Ω_{i}(F). Let ${\mathrm{\Omega}}_{i}^{}\left(F\right)$ denote the regions in which the function F(ω,t) takes negative values and which are therefore the candidates of the ocean domain. Since we consider the ocean to be an interconnected water volume, termed the global ocean, as in traditional physical oceanography and sealevel studies, only the largest amongst ${\mathrm{\Omega}}_{i}^{}\left(F\right)$ forms the ocean domain. Smaller ${\mathrm{\Omega}}_{i}^{}\left(F\right)$, if there are any, and their boundaries T_{i}(F) are considered to be part of the land, meaning they are unable to freely exchange mass with the global ocean by GRD processes. One obvious example of the region that does not belong to the ocean in spite of having $F(\mathit{\omega},t)<\mathrm{0}$ is a continental trough with bathymetry below MSL. Unless this trough is physically connected to the global ocean via oceanic water, we consider this to be part of the land rather than the ocean. Let ${\mathrm{\Omega}}_{S}^{}\left(F\right)$ be the union of all these small nonoceanic regions and T_{S}(F) be the union of corresponding boundaries. We modify F(ω,t) to define a new function,
so that its zerolevel set
represents the land–ocean boundaries. Here ϵ is a positive number to ensure $\mathcal{F}(\mathit{\omega},t)>\mathrm{0}$ at ω∈T_{S}(F).
The land–ocean boundaries are generally known as coastlines. Given the definition of the levelset function (Eq. 3), coastlines are free of ice where $B(\mathit{\omega},t)=S(\mathit{\omega},t)$. No coastline exists with $B(\mathit{\omega},t)>S(\mathit{\omega},t)$. Only in the marine sectors where $B(\mathit{\omega},t)<S(\mathit{\omega},t)$ does a coastline have finite ice thickness and is then replaced by the term grounding line.
2.2 Definitions of land, ocean and ice domains
Given the definition of coastlines and grounding lines, we may define the ocean domain as follows:
The land domain is simply given by $\mathcal{L}(\mathit{\omega},t)=\mathrm{1}\mathcal{O}(\mathit{\omega},t)$. Surface areas of these complementary domains together make up
the total area of the planetary surface, a necessary condition for mass conservation in the Earth system. Note that neither 𝒪(ω,t)
nor ℒ(ω,t) is defined at the coastlines or grounding lines. These interfaces rather form their own level set,
𝒯(ℱ), as defined in Sect. 2.1. For practical purposes, however, one may carry all these masks and level sets as a single
field. For example, the Icesheet and Sealevel System Model (ISSM; https://issm.jpl.nasa.gov/, last access: 31 August 2020) uses the field md.mask.ocean_levelset
, which takes −1 in the ocean, 1 in land, and 0 at the coastlines and grounding lines.
We define ℐ(ω,t) to be a globally distributed system of ice domains, such that
For many applications, it may be useful to decompose ℐ(ω,t) into a number of subdomains: $\mathcal{I}(\mathit{\omega},t)=\mathit{\{}{\mathcal{I}}_{\mathrm{1}},{\mathcal{I}}_{\mathrm{2}},\mathrm{\dots},{\mathcal{I}}_{i},\mathrm{\dots}\mathit{\}}$, where ℐ_{i}(ω,t) represents the ith ice domain. Individual ice sheets and glaciers can be thought of individual ice domains. As defined in Sect. 2.1, the grounding line within a given ice domain is given by the levelset 𝒯(ℱ). Using Eqs. (5) and (6), we may define the grounded ice mask simply as $\mathcal{G}(\mathit{\omega},t)=\mathcal{I}(\mathit{\omega},t)\mathcal{L}(\mathit{\omega},t)$ and the floating ice mask as $\mathcal{I}(\mathit{\omega},t)\mathcal{O}(\mathit{\omega},t)$. Equation (6) is a simple, and perhaps the most generic, definition for ice domains, which can accommodate any geometric features such as kilometerscale pinning points or rugged fjords that can modulate marine icesheet instability on retrograde slopes (e.g., Matsuoka et al., 2015; Whitehouse et al., 2017). The employed definition of floating ice mask, however, limits us from capturing the floating ice on subglacial and proglacial lakes that are not part of the global ocean (Fig. 1). We believe that the localized processes of ice–lake interactions are of secondary importance, at least, for the purpose of capturing largescale interplay between the continental ice sheets, solid Earth and sea level in the current Earth system models.
Our definition of coastlines and grounding lines, and hence that of the land and ocean and the grounded and floating ice, facilitates direct evaluation of the interaction between a dynamic system of ice sheets, solid Earth and sea level, as well as the estimation and interpretation of icedriven global and regional sealevel change by conserving mass in the Earth system. Although a distributed system of ice domains is an integral part of the Earth system, in the following we consider, for brevity, a single domain as an ice sheet, while other ice domains are collectively referred to as farfield ice.
The estimation of the sealevel contribution from an evolving ice sheet, featuring marinebased grounded and floating ice, is not trivial, particularly in light of evolving bedrock and MSL. Here we review the common method and its limitations and present a new method that is applied to arbitrary ice geometries, all kinds of bedrock and MSL forcings, and at all timescales.
3.1 Change in ice height above floatation
We use the bedrock and MSL to define a floatation height for ice:
such that the ice thickness in excess of H_{0}(ω,t) represents the socalled height above floatation (HAF). For convenience, we define HAF only in the grounded ice domain:
so that we may interpret it as the fraction of landice thickness that can potentially contribute to sealevel change (Fig. 2a). It is clear from the equations that ${H}_{\mathrm{F}}(\mathit{\omega},t)=H(\mathit{\omega},t)$ for the grounded ice sheet that rests on the bedrock whose elevation is at or above MSL. For grounded portions of the marine ice sheet, ${H}_{\mathrm{F}}(\mathit{\omega},t)<H(\mathit{\omega},t)$ and, in fact, H_{F}(ω,t) can take negative values. If the ice sheet were to disintegrate, for instance, regions with ${H}_{\mathrm{F}}(\mathit{\omega},t)<\mathrm{0}$ such as Sector D in Fig. 2a can take up ocean water and contribute to sealevel fall.
The evolving icesheet geometry is usually described in terms of ice thickness and icesheet margins. Indeed, prognostic simulations of icesheet models track the transport of mass in terms of equivalent icethickness distribution. The transport of mass within the ice domain and ice–ocean mass exchange induce a GRD response of the solid Earth, which further redistributes the mass in the Earth system. This modulates the bedrock topography as well as MSL. These evolving fields may also have components that are forced by external processes such as contemporaneous melting of farfield ice, or GIA, or tectonics. We may describe the evolving icesheet geometry in terms of ice thickness, bedrock elevation and MSL (see Fig. 2b–e). We denote ΔH(ω,Δt), ΔB(ω,Δt) and ΔS(ω,Δt) to be the change in respective fields over the time interval Δt. For the new icesheet geometry at time t+Δt, Eq. (8) gives ${H}_{\mathrm{F}}(w,t+\mathrm{\Delta}t)$ = $\mathcal{G}(\mathit{\omega},t+\mathrm{\Delta}t)\left[H\right(\mathit{\omega},t)$ + ΔH(ω,Δt) − ${H}_{\mathrm{0}}(\mathit{\omega},t+\mathrm{\Delta}t)]$, where ${H}_{\mathrm{0}}(\mathit{\omega},t+\mathrm{\Delta}t)$ = ρ_{o}∕ρ_{i} $max\left[\mathit{\right\{}S(\mathit{\omega},t)+\mathrm{\Delta}S(\mathit{\omega},\mathrm{\Delta}t)B(\mathit{\omega},t)\mathrm{\Delta}B(\mathit{\omega},\mathrm{\Delta}t)\mathit{\}},\mathrm{0}]$ is given by Eq. (7). Similarly, we may rewrite Eq. (1) for $F(\mathit{\omega},t+\mathrm{\Delta}t)$ and define the new ocean domain $\mathcal{O}(\mathit{\omega},t+\mathrm{\Delta}t)$, land domain $\mathcal{L}(\mathit{\omega},t+\mathrm{\Delta}t)$ and grounded ice domain $\mathcal{G}(\mathit{\omega},t+\mathrm{\Delta}t)$ as described in Sect. 2.
In what follows, we assume that the net change in grounded ice mass results in the equivalent change in ocean mass, ensuring mass conservation in the Earth system. On the one hand not all of ΔH(ω,Δt) contributes to change in mass of oceanic water, but on the other hand, in response to the externally forced bedrock and MSL change, the ice sheet may still contribute to change in ocean mass even when $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{0}$ as we count floating ice in the ocean mass (see Appendix A). The standalone icesheet models evaluate the change in HAF in order to calculate the sealevel contribution of an ice sheet, termed the HAF method for brevity (e.g., Bindschadler et al., 2013; Nowicki et al., 2016):
These models generally (but incorrectly) calculate the equivalent oceanicwater volume (rather than the freshwater volume) by spatially integrating $[{\mathit{\rho}}_{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{o}}]\mathrm{\Delta}{H}_{\mathrm{F}}(\mathit{\omega},\mathrm{\Delta}t)$ and divide it by the ocean surface area to estimate the GMSL change. Apart from this waterdensityrelated error, the HAF method in absence of evolving bedrock and MSL yields the correct estimates of the sealevel contribution (see Appendix A). In fact, the effects of ΔB(ω,t) and ΔS(ω,t) may be negligible over the timescale of a few decades or shorter. The standalone icesheet models typically inherit this assumption, even though simulation timescales can be on the order of centuries. Over such relatively longer timescales, this simplistic approach yields some error, especially in the marine portions of an ice sheet (Larour et al., 2019; Goelzer et al., 2020).
3.2 A new field for estimating sealevel contribution
In order to overcome the limitations of the HAF method, we define a unified field, ΔH_{S}(ω,Δt), that contributes to sealevel change by modulating both the mass and volume of oceanic water over the period Δt. This field captures the effects of evolving bedrock and MSL induced by any geophysical processes and is applied to arbitrary ice geometries and at all timescales. We find it convenient to partition ΔH_{S}(ω,Δt) upfront into two components:
such that the first component ΔH_{M}(ω,Δt) modulates both the mass and volume of the oceanic water, while the second component ΔH_{V}(ω,Δt) can only modulate the ocean volume. The following relationships hold for generalized ice geometries, as well as bedrock and MSL forcings:
where ρ_{w} is the freshwater density. For grounded ice sheets, the mass component makes up about 97 % of ΔH_{S}(ω,Δt), which loads the solid Earth and induces its GRD response and sealevel adjustment, which will be further discussed in Sect. 4.
While a detailed interpretation of individual terms appearing in Eqs. (11) and (12) is given in Appendix A by considering all possible scenarios of evolving ice thickness, bedrock elevation and MSL, Fig. 2 illustrates a few representative scenarios. In reference to this figure and Eqs. (11) and (12), we outline three distinct regimes:

Regime 1: where ice remains grounded at both times t and t+Δt.
All of ΔH(ω,Δt) in this regime contributes to sealevel change by modulating both the mass and volume of the ocean (first term on the right side of Eq. 11), irrespective of the elevation of bedrock upon which the ice is grounded. It turns out $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)\ne \mathrm{\Delta}{H}_{\mathrm{F}}(\mathit{\omega},\mathrm{\Delta}t)$ only in the marine portions of the regime and only when evolving bedrock and MSL are considered (see Appendix A1). Goelzer et al. (2020) present a method to backtrack ΔH(ω,Δt) from ΔH_{F}(ω,Δt) in such situations, assuming that ΔB(ω,Δt) and ΔS(ω,Δt) are known.
This regime also includes land areas covered by the evolving icesheet margins. When ice margin advances over the period Δt, newly glaciated areas must satisfy $H(\mathit{\omega},t)=\mathrm{0}$ and $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)>\mathrm{0}$. When it retreats, $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)=H(\mathit{\omega},t)$ must hold in the recently deglaciated areas. In both cases, all of ΔH(ω,Δt) contributes to the sealevel change.
Externally forced ΔB(ω,Δt) or ΔS(ω,Δt) does not affect the estimate of ΔH_{S}(ω,Δt) in this regime although it may alter bedrock slope or gravitational driving stress and possibly modulate the iceflow dynamics. While the effects of farfield ice melting and associated ocean loading may be negligible due to their relatively long wavelength imprints, ΔB(ω,Δt) due to large earthquakes beneath the ice sheet may have some impact on ice dynamics.

Regime 2: where ice transitions from grounded to floating, or the reverse, over the period Δt.
The sealevel contribution from this regime mainly depends on the change in HAF (second term on the right side of Eq. 11), which modulates both the mass and volume of oceanic water. In the absence of externally forced bedrock and MSL, it follows that $\left\mathrm{\Delta}{H}_{\mathrm{F}}\right(\mathit{\omega},\mathrm{\Delta}t\left)\right<\left\mathrm{\Delta}H\right(\mathit{\omega},\mathrm{\Delta}t\left)\right$ (see Appendix A2). The change in ice thickness in excess of the change in HAF (rightside term in Eq. 12) nominally modulates the volume of the oceanic water. This is due to the difference in volume between the freshwater that would be produced when ice melts and the oceanic water that would be displaced when it floats.
This is the only regime where, in response to the externally forced ΔB(ω,Δt) or ΔS(ω,Δt), an ice sheet may modulate both the mass and volume of the ocean even when $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{0}$. Specific examples are given in Appendix A2.

Regime 3: where ice remains floating at both times t and t+Δt.
Since $\mathrm{\Delta}{H}_{\mathrm{F}}(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{0}$ holds true in this regime, the change in ice thickness does not modulate the ocean mass itself, but it releases or takes up freshwater that has slightly larger volume than the oceanic water upon which it floats. This minor difference in water volume (rightside term in Eq. 12) contributes to the sealevel change (see Appendix A3).
Given the new field for estimating sealevel contribution (Eq. 10), we may readily calculate the GMSL change by spatially integrating $[{\mathit{\rho}}_{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{w}}]\mathrm{\Delta}{H}_{\mathrm{S}}(\mathit{\omega},\mathrm{\Delta}t)$, which yields the total freshwater volume being added to the ocean over the period Δt, and dividing it by the ocean surface area at time t+Δt. Assume that an ice sheet collapses instantaneously and that all of the meltwater makes it to the ocean. The resulting GMSL change represents the potential sea level of the ice sheet at time t, and it can be readily derived from Eq. (10) by setting $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)=H(\mathit{\omega},t)$ and $\mathcal{G}(\mathit{\omega},t+\mathrm{\Delta}t)=\mathrm{0}$ in the limit of Δt→0. Note that ΔH(ω,Δt) and $\mathcal{G}(\mathit{\omega},t+\mathrm{\Delta}t)$ are implicit via Eqs. (9), (11) and (12).
3.3 Quantitative comparison of the two methods
Here we present a case study to demonstrate the level of improvements possible by employing the new method (Eq. 10) over the HAF method (Eq. 9). For a quantitative comparison, we rely on the recent work of Larour et al. (2019), who provide consistent solutions of evolving H(ω,t), B(ω,t) and S(ω,t) for the Antarctic Ice Sheet over the next 500 years. They simulate a highresolution dynamical iceflow model (Larour et al., 2012) that is fully coupled with a global solidEarth deformation and sealevel adjustment model (Adhikari et al., 2016) under the presentday surface climatology and a realistic subiceshelf melting scenario (Seroussi et al., 2017). They also account for the effects of farfield icemass change on the evolution of bedrock and MSL in Antarctica. To this end, they consider mass balance of the Greenland Ice Sheet and global glaciers, extrapolated into the next 500 years based on the spacegravimetrybased measurements. The ongoing change in bedrock and MSL due to the viscous response of the solid Earth to the global deglaciation since the Last Glacial Maximum is also accounted for through an offline coupling of a GIA model (Caron et al., 2018).
In Fig. 3, we compare the two methods both in terms of their spatial and temporal patterns. We show ΔH_{S}(ω,Δt) and ΔH_{F}(ω,Δt) computed at 2350 CE relative to 2000 CE for Pine Island and Thwaites glaciers (Fig. 3a–c). To facilitate the interpretation, we separate the model domain into three regimes as in Fig. 2b: Regime 1 (Regime 3) consists of regions that are grounded (floating) at both times; and Regime 2 consists of regions that transition from grounded to floating over the course of simulation. Only ΔH(ω,Δt) contributes to ΔH_{S}(ω,Δt) in Regime 1 (see Eq. 11). In the marine portions of this regime, ΔH(ω,Δt) and hence ΔH_{S}(ω,Δt) differ from ΔH_{F}(ω,Δt) due to the effects of evolving bedrock and MSL in the latter field. The difference between ΔH_{S}(ω,Δt) and ΔH_{F}(ω,Δt) in Regime 2 and Regime 3 are due to ΔH_{V}(ω,Δt) (Eq. 12), which accounts for the volumetric contribution of icethickness change in excess of the change in HAF (see Appendix A2). We also show the time series of the total Antarctic icevolume change that is attributable to the GMSL change (Fig. 3d and e). We find that the new method predicts systematically larger sealevel contribution, compared to the HAF method, throughout the model simulation. The difference in the first 100 years is more than 15 %, and in the last 100 years it is about 8 %–10 %. We isolate the mass (Eq. 11) and the volume component (Eq. 12) of the new method to show that the former component alone, which drives the GRD response of the solid Earth, consistently predicts larger sealevel contribution than the HAF method by about 5 %. Note that the HAF method usually converts the icevolume change (Fig. 3d) into the equivalent oceanic water (rather than the freshwater) volume change (e.g., Bindschadler et al., 2013) and hence systematically underpredicts the amplitude of GMSL change by an additional 2 %–3 %, which is not accounted for in the above analysis.
Delineation of evolving coastlines (and grounding lines) and estimation of ΔH_{S}(ω,Δt) require the knowledge of ΔH(ω,Δt), ΔB(ω,Δt) and ΔS(ω,Δt) along with accurate information of the solid Earth surface (e.g., bedrock topography in ice domains, and landsurface topography and ocean bathymetry in the vicinity of coastlines). Parts of ΔB(ω,Δt) and ΔS(ω,Δt) are induced by ΔH_{S}(ω,Δt) and associated ocean mass change themselves. These fields are therefore intertwined with each other, and only by using a massconserving Earth system model that can capture icesheet dynamics, solidEarth deformation and sealevel adjustment may we find selfconsistent solutions. We find it convenient to treat the change in bedrock and MSL collectively in terms of the change in relative sea level (RSL), which by definition is the MSL relative to the sea floor or bedrock (Gregory et al., 2019). Mathematically, $\mathrm{\Delta}R(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{\Delta}S(\mathit{\omega},\mathrm{\Delta}t)\mathrm{\Delta}B(\mathit{\omega},\mathrm{\Delta}t)$. Since ΔR(ω,Δt) may be induced by processes other than ΔH_{S}(ω,Δt), we must consider them as they impact the estimate of evolving 𝒯(ℱ), and hence the ocean surface area, and ΔH_{S}(ω,Δt) itself. In fact, we should also consider the change in steric MSL, which is not accounted for in ΔS(ω,Δt) (see Sect. 2.1). The inclusion of this component, however, must be accompanied by the spatially and temporally varying ocean density (see, for example, Eq. 1), which modulates the coastlines, and hence the ocean surface area, but does not affect the buoyant force on the ice, groundingline migration, and the estimate of ΔH_{S}(ω,Δt).
To further diagnose ΔR(ω,Δt), especially in light of spacebased observations and existing GRD models, we present a synopsis of contributing processes as follows:
The first four terms on the right side of the equation represent the processes that exchange water between the land and the ocean and contribute to sealevel change by inducing the GRD response of the solid Earth (termed, for brevity, the barystatic components). We use the superscript I to refer to the ice sheet under consideration and superscript L to refer to other parts of the land, including farfield ice and hydrological basins. When these sources of freshwater contribute to sealevel change over a contemporaneous period $[t,t+\mathrm{\Delta}t]$, corresponding changes in RSL are denoted with the subscript C. The land–ocean water exchange may have occurred in the past, i.e., over the period $(\mathrm{\infty},t]$, and the induced viscous response of the solid Earth may still contribute to the RSL change over the period $[t,t+\mathrm{\Delta}t]$. These components are denoted with the subscript P. The last term appearing in the equation captures other nonbarystatic processes that may or may not induce GRD response of the solid Earth but at least modulate the ocean bathymetry or coastal geometry. These processes include earthquakes, landslides, sediment transport and coastal subsidence, amongst others. Assuming that the contemporaneous period Δt is on the order of 10 years, we may interpret ΔR(ω,Δt) as the nonsteric component of ongoing RSL change monitored by the satellite gravimetry and altimetry (WCRP Global Sea Level Budget Group, 2018). We may interpret $\mathrm{\Delta}{R}_{\mathrm{C}}^{\mathrm{I}}(\mathit{\omega},\mathrm{\Delta}t)+\mathrm{\Delta}{R}_{\mathrm{C}}^{\mathrm{L}}(\mathit{\omega},\mathrm{\Delta}t)$ and $\mathrm{\Delta}{R}_{\mathrm{P}}^{\mathrm{I}}(\mathit{\omega},\mathrm{\Delta}t)+\mathrm{\Delta}{R}_{\mathrm{P}}^{\mathrm{L}}(\mathit{\omega},\mathrm{\Delta}t)$ as ongoing RSL change driven by contemporary global surface mass redistribution (e.g., Adhikari et al., 2019) and by global GIA processes (e.g., Peltier et al., 2015; Caron et al., 2018), respectively.
As defined in Eq. (10), only part of ΔH_{S}(ω,Δt) potentially contributes to the ocean mass change, loads the underlying solid Earth, and induces its GRD response and contributes to sealevel adjustment. Because the GRD effect is applied to the entire column of the ocean water, only $[{\mathit{\rho}}_{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{o}}]\mathrm{\Delta}{H}_{\mathrm{M}}(\mathit{\omega},\mathrm{\Delta}t)\approx \mathrm{0.892}\times \mathrm{\Delta}{H}_{\mathrm{M}}(\mathit{\omega},\mathrm{\Delta}t)$ induces the barystatic component of the RSL change, which in reference to Eq. (13) is equivalent to $\mathrm{\Delta}{R}_{\mathrm{C}}^{\mathrm{I}}(\mathit{\omega},\mathrm{\Delta}t)$. The freshwater equivalent of other parts of ΔH_{S}(ω,Δt), i.e., the sum of $[{\mathit{\rho}}_{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{w}}{\mathit{\rho}}_{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{o}}]\mathrm{\Delta}{H}_{\mathrm{M}}(\mathit{\omega},\mathrm{\Delta}t)\approx \mathrm{0.025}\times \mathrm{\Delta}{H}_{\mathrm{M}}(\mathit{\omega},\mathrm{\Delta}t)$ and $[{\mathit{\rho}}_{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{w}}]\mathrm{\Delta}{H}_{\mathrm{V}}(\mathit{\omega},\mathrm{\Delta}t)$, contributes to the RSL change by modulating oceanicwater density, and hence it may be considered part of the steric MSL. For grounded ice sheets, this component of RSL change is about 97 % smaller than the barystatic component. The remaining terms of Eq. (13) are what we collectively refer to as the “externally forced” RSL change. In other words, these are the RSL components not directly induced or contributed by the ice sheet under consideration over the period Δt.
To solve for the spatial pattern of the steric MSL due to ΔH_{S}(ω,Δt), one must consider a dynamic ocean circulation model. Such computations are generally not warranted in the longerterm (decadal or longer timescale) sealevel studies, owing to their smaller amplitudes compared to those of the barystatic component. The spatial pattern of the barystatic RSL due to ΔH_{S}(ω,Δt) can be obtained by solving the socalled “sealevel equation” on a selfgravitating, viscoelastically compressible, rotating Earth (Farrell and Clark, 1976; Milne and Mitrovica, 1998). To this end, we must consider a massconserving field that describes the net change in mass per unit area on the solid Earth surface:
such that its global integral is zero. Here ΔH_{M}(ω,Δt) is given by Eq. (11) and $\mathrm{\Delta}{R}_{\mathrm{C}}^{\mathrm{I}}(\mathit{\omega},\mathrm{\Delta}t)$ is precisely the same as the first term on the right side of Eq. (13). Because the RSL is defined globally, including in land, we must invoke the ocean mask in the equation. Solving the sealevel equation, in essence, means that we load the solid Earth by the massconserving surface load (Eq. 14) and let its GRD response dictate the selfconsistent patterns of RSL, MSL and bedrock changes as well as the new positions of coastlines and grounding lines. The icesheet models that account for local or regional isostatic adjustment of bedrock (e.g., Le Meur and Huybrechts, 1996; Bueler et al., 2007; Pattyn, 2017) generally do not consider $\mathrm{\Delta}{R}_{\mathrm{C}}^{\mathrm{I}}(\mathit{\omega},\mathrm{\Delta}t)$ as part of the surface load. As a result, these models violate mass conservation in the Earth system and capture incomplete signals of ΔB(ω,Δt) and ΔS(ω,Δt) in the estimation of ΔH_{S}(ω,Δt).
We have divided the Earth's surface into complementary domains of land and ocean, which are separated by coastlines. While there may be multiple land domains, we maintain a single global ocean of interconnected oceanic water as in the majority of studies in physical oceanography and sea level. Distributed bodies of ice intersect the land and the ocean to form glaciers, ice sheets and ice shelves. Grounding lines are defined as the coastlines that belong to the ice domains. The set of generic, and quite simple, mathematical descriptions presented here can handle the complex geometries of both the coastlines and grounding lines, complementary to those of farfield land, ocean and ice domains and their respective overall evolutionary history.
Based on this formalism of evolving coastlines and grounding lines, we present a unified method to calculate the exact fraction of ice thickness that contributes to the sealevel change over a given period. The method is a function of evolving ice thickness, bedrock elevation and mean sea level driven by any geophysical processes. Along with its obvious application to estimate the globalmean sealevel change, it is absolutely critical to track the global mass transport, and assess the response of a dynamic system of ice sheets, solid Earth and sea level, while accounting for kilometerscale features in ice–bedrock–ocean geometries. Our method requires bookkeeping of the global land and ocean domains. This is crucial for considering distributed ice and land domains with complex geometries in Earth system models. For the treatment of an individual ice sheet, however, it is sufficient to track a continental or regional land domain, provided that there is an understanding that ocean water may recede from, or reinundate, continental land (e.g., Johnston, 1993; Milne, 1998). In fact, it may often be possible to use grounded ice masks in place of the land domains. In the most simplified case when the bedrock and mean sea level do not evolve, our method reduces to the common method that is based on the concept of ice height above floatation. For an example model simulation considered in this study (Larour et al., 2019), we find that the new method systematically yields 10 %–15 % more sealevel contribution from the Antarctic Ice Sheet. We recommend that the icesheet modeling community consider the proposed method as a metric to quantify the sealevel contribution of evolving ice sheets. This is especially appropriate for model analysis that is informed by ice and ocean mass monitoring from space assets, such as ocean and ice altimetry, radar interferometry, and space gravimetry (e.g., Bentley and Wahr, 1998).
Here we provide an indepth comparison between the HAF method (Eq. 9) and ours (Eq. 10) in light of evolving ice thickness, bedrock elevation and MSL. The latter two fields may be treated collectively in terms of the relative sea level (RSL) which, by definition, is the MSL relative to the bedrock or sea floor. In the following, we consider all plausible scenarios by combining the change in ice thickness, ΔH(ω,Δt), and relative sea level, ΔR(ω,Δt), over the period Δt.
A1 Where ice remains grounded at both times t and t+Δt
In our method, all of ΔH(ω,Δt) irrespective of ΔR(ω,Δt) contributes to sealevel change by modulating both the mass and volume of the oceanic water (see the first term on the right side of Eq. 11). The same is true for the HAF method as long as ice remains grounded on the bedrock whose elevation is at or above MSL at both times t and t+Δt, in which case $\mathrm{\Delta}{H}_{\mathrm{F}}(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)$. There is nonetheless a minor densityrelated difference between the two methods: our method evaluates the freshwater equivalent height $[{\mathit{\rho}}_{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{w}}]\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)$, whereas the HAF method generally evaluates the oceanicwater equivalent height $[{\mathit{\rho}}_{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{o}}]\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)$. As a result, the HAF method systematically underestimates the amplitude of the globalmean sealevel (GMSL) change by about 2 % to 3 %.
If the ice is grounded on the marine bedrock (whose elevation is below MSL) at least at time t or t+Δt, the HAF method generally yields incorrect solution in addition to the densityrelated error noted above. In this case, $\mathrm{\Delta}{H}_{\mathrm{F}}(\mathit{\omega},\mathrm{\Delta}t)\ne \mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)$ generally holds true, and depending upon the relative amplitudes and signs of ΔH(ω,Δt) and ΔR(ω,Δt), the HAF method may over or underpredict the sealevel contribution compared to our method. Two special cases are worth mentioning:

Case A1.1: $\mathrm{\Delta}R(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{0}$. When the effect of evolving RSL is not considered, we find that $\mathrm{\Delta}{H}_{\mathrm{F}}(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)$ and, consequently, the two methods are equivalent.

Case A1.2: $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{0}$. When the thickness of grounded ice does not change, the HAF method may incorrectly predict nonzero sealevel contribution in cases when $\mathrm{\Delta}R(\mathit{\omega},\mathrm{\Delta}t)\ne \mathrm{0}$ (see Fig. A1a). In this case, the HAF method systematically overestimates (underestimates) GMSL change when ΔR(ω,Δt) is greater (less) than zero.
A2 Where ice transitions from grounded to floating, or the reverse, over the period Δt
Here the working principle of both methods is same and as follows. We derive the potential sea level (PSL) contributions of ice thicknesses at time t and t+Δt. We then compute their difference to derive the actual sea level (ASL) contribution over the period Δt. We may define PSL and ASL in terms of freshwater equivalent height as follows (see Fig. A1):
Ice equivalent height of the terms appearing on the right side of Eq. (A3) are reported in the main text: the second term on the right side of Eq. (11) and the term on the right side of Eq. (12), respectively. Since the HAF method deals with the oceanicwater density, rather than the freshwater density, both PSL and ASL in this method can be deduced from the above equations by replacing ρ_{i}∕ρ_{w} by ρ_{i}∕ρ_{o}. The second terms in these equations vanish, and the ASL is given in terms of oceanicwater equivalent height by [ρ_{i}∕ρ_{o}]ΔH_{F}(Δt), whose ice equivalent height is reported in the main text (Eq. 9).
When ice transitions from grounded to floating, ${H}_{\mathrm{F}}(t+\mathrm{\Delta}t)=\mathrm{0}$ and the first term appearing in Eq. (A3) always takes a negative value. Three distinct scenarios of evolving ice thickness and RSL may be of interest:

Case A2.1: $\mathrm{\Delta}R(\mathit{\omega},\mathrm{\Delta}t)\le \mathrm{0}$ and $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)<\mathrm{0}$. In this case, the only condition for the grounded ice to float is through its sufficient thinning such that $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)<\mathrm{\Delta}{H}_{\mathrm{F}}(\mathit{\omega},\mathrm{\Delta}t)<\mathrm{0}$ (see Fig. A1b). Both terms appearing in Eq. (A3) take negative values, causing the GMSL to rise.

Case A2.2: $\mathrm{\Delta}R(\mathit{\omega},\mathrm{\Delta}t)>\mathrm{0}$ and $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{0}$. Here the condition $\mathrm{\Delta}{H}_{\mathrm{F}}(\mathit{\omega},\mathrm{\Delta}t)<\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{0}$ holds true (see Fig. A1c). Since the first term appearing in Eq. (A3) is about 97 % larger in magnitude than the second term (which takes a positive value), it causes the GMSL to rise. In other words, the externally forced RSL rise causes the ice to contribute to GMSL rise even though its thickness does not change.

Case A2.3: $\mathrm{\Delta}R(\mathit{\omega},\mathrm{\Delta}t)>\mathrm{0}$ and $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)\ne \mathrm{0}$. Only when $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)>\mathrm{0}$ and its amplitude is significantly larger (by about a factor of 35) than that of ΔH_{F}(ω,Δt), the second term appearing in Eq. (A3) that takes a positive value dominates and causes the GMSL to fall. Otherwise, the GMSL rises even when ice thickens over the period Δt.
When ice transitions from floating to grounded, H_{F}(t)=0 and the first term appearing in Eq. (A3) always takes a positive value. Three distinct scenarios of evolving ice thickness and RSL may be of interest:

Case A2.4: $\mathrm{\Delta}R(\mathit{\omega},\mathrm{\Delta}t)\ge \mathrm{0}$ and $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)>\mathrm{0}$. In this case, the only condition for the floating ice to be grounded is through its sufficient thickening such that $\mathrm{0}<\mathrm{\Delta}{H}_{\mathrm{F}}(\mathit{\omega},\mathrm{\Delta}t)<\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)$. Both terms appearing in Eq. (A3) take positive values, causing the GMSL to fall.

Case A2.5: $\mathrm{\Delta}R(\mathit{\omega},\mathrm{\Delta}t)<\mathrm{0}$ and $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{0}$. Here the condition that $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{0}<\mathrm{\Delta}{H}_{\mathrm{F}}(\mathit{\omega},\mathrm{\Delta}t)$ holds true. Since the first term appearing in Eq. (A3) is about 97 % larger in magnitude than the second term (which takes a negative value), it causes the GMSL to fall. In other words, the externally forced RSL drop causes the ice to further contribute to GMSL drop even though its thickness does not change.

Case A2.6: $\mathrm{\Delta}R(\mathit{\omega},\mathrm{\Delta}t)<\mathrm{0}$ and $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)\ne \mathrm{0}$. Only when $\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)<\mathrm{0}$ and its amplitude is significantly larger (by about a factor of 35) than that of ΔH_{F}(ω,Δt), the second term appearing in Eq. (A3) that takes a negative value dominates and causes the GMSL to rise. Otherwise, the GMSL falls even when ice thins over the period Δt.
A3 Where ice remains floating at both times t and t+Δt
We may evaluate PSL and ASL contributions from this region based on Eqs. (A1)–(A3). Since H_{F}=0 at both times t and t+Δt, the evaluation of the sealevel contribution does not depend on the evolving RSL. In this scenario, the ASL is given in terms of freshwater equivalent height by $[{\mathit{\rho}}_{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{w}}{\mathit{\rho}}_{\mathrm{i}}/{\mathit{\rho}}_{\mathrm{o}}]\mathrm{\Delta}H(\mathit{\omega},\mathrm{\Delta}t)$, whose ice equivalent height can be deduced from Eq. (12) by setting $\mathrm{\Delta}{H}_{\mathrm{F}}(\mathit{\omega},\mathrm{\Delta}t)=\mathrm{0}$ (see Fig. A1d). When the ice thins (thickens), it causes the GMSL to rise (fall) by modulating the volume (not mass) of the oceanic water. The ASL contribution in the HAF method can be deduced by replacing ρ_{i}∕ρ_{w} by ρ_{i}∕ρ_{o} and is effectively zero and, therefore, does not appear explicitly in Eq. (9). This suggests that the HAF method systematically underestimates the amplitude of GMSL change.
B  Solid Earth surface (i.e., land surface or sea floor), or simply bedrock elevation 
ΔB  Change in bedrock elevation over the period Δt 
F  A function such that F=0 satisfies the hydrostatic equilibrium between ice and the oceanic water 
ℱ  A function such that ℱ=0 represents the grounding lines or coastlines 
𝒢  Mask of the grounded portions of an ice sheet 
H  Ice thickness 
ΔH  Change in ice thickness over the period Δt 
H_{0}  Floatation height for ice 
H_{F}  Ice height above floatation (HAF) defined for grounded ice 
ΔH_{F}  Change in HAF over the period Δt 
ΔH_{M}  A component of ΔH_{S} that modulates both the mass and volume of the ocean over the period Δt 
ΔH_{S}  A new field for estimating the sealevel contribution of ice sheets over the period Δt 
ΔH_{V}  A component of ΔH_{S} that can only modulate the ocean volume over the period Δt 
ℐ  A globally distributed system of ice domains 
ℒ  A globally distributed system of land domains 
ΔM  Change in mass per unit area on the solid Earth surface over the period Δt 
𝒪  The global ocean domain 
ΔR  Nonsteric component of the relative sea level (RSL) change over the period Δt 
$\mathrm{\Delta}{R}_{\mathrm{C}}^{\mathrm{I}}$  A component of ΔR due to the contemporary icesheet mass change 
$\mathrm{\Delta}{R}_{\mathrm{C}}^{\mathrm{L}}$  A component of ΔR due to the contemporary change in farfield ice and land water storage 
$\mathrm{\Delta}{R}_{\mathrm{P}}^{\mathrm{I}}$  A component of ΔR due to the past icesheet mass change 
$\mathrm{\Delta}{R}_{\mathrm{P}}^{\mathrm{L}}$  A component of ΔR due to the past change in farfield ice and land water storage 
ΔR_{O}  A component of ΔR due to other nonsteric processes 
ρ_{i}  Density of ice 
ρ_{o}  Mean density of the oceanic water 
ρ_{w}  Density of freshwater 
S  Mean sea level (MSL), excluding its steric component 
ΔS  Change in S over the period Δt 
T  The zerolevel set of F that satisfies the hydrostatic equilibrium between ice and the oceanic water 
𝒯  The zerolevel set of ℱ that represents the coastlines and grounding lines 
t  Time 
Δt  Time period 
ω  2D spatial coordinates, geographic (θ,ϕ) or Cartesian (x,y), on the planetary surface 
The data and code used to produce Fig. 3 are available online at https://doi.org/10.7910/DVN/9LUJTD (Adhikari et al., 2020).
SA conceived and conducted the research. All authors contributed to the theory. ERI helped write the first draft of the manuscript. EL helped produce Fig. 3. All authors contributed to the writing and editing of the manuscript.
The authors declare that they have no conflict of interest.
This research was carried out at the Jet Propulsion Laboratory (JPL), California Institute of Technology, under a contract with National Aeronautics and Space Administration (NASA). We acknowledge three anonymous reviewers for their thorough and constructive reviews that improved the manuscript.
This research has been funded through the JPL Research, Technology and Development programs (grant nos. 01STCRR.17.235.118, 2017–2019; and 01STCRR.19.021.241, 2019), the NASA SeaLevel Change Team (grant no. 16SLCT160015, 2018–2020), and the NASA Earth Surface and Interior Program (grant no. 19ESI190020, 2020–2022).
This paper was edited by Alexander Robinson and reviewed by three anonymous referees.
Adhikari, S., Ivins, E. R., Larour, E., Seroussi, H., Morlighem, M., and Nowicki, S.: Future Antarctic bed topography and its implications for ice sheet dynamics, Solid Earth, 5, 569–584, https://doi.org/10.5194/se55692014, 2014. a
Adhikari, S., Ivins, E. R., and Larour, E.: ISSMSESAW v1.0: meshbased computation of gravitationally consistent sealevel and geodetic signatures caused by cryosphere and climate driven mass change, Geosci. Model Dev., 9, 1087–1109, https://doi.org/10.5194/gmd910872016, 2016. a
Adhikari, S., Ivins, E. R., Frederikse, T., Landerer, F. W., and Caron, L.: Sealevel fingerprints emergent from GRACE mission data, Earth Syst. Sci. Data, 11, 629–646, https://doi.org/10.5194/essd116292019, 2019. a
Adhikari, S. Ivins, E. Larour, E. Caron, L., and Seroussi, H.: Sample data for computing sealevel contribution from ice sheets, https://doi.org/10.7910/DVN/9LUJTD, Harvard Dataverse, V1, 2020. a
Altamimi, Z., Rebischung, P., Metivier, L., and Collilieux, X.: ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions, J. Geophys. Res.Sol. Ea., 121, 6109–6131, 2016. a
Barletta, V. R., Bevis, M., Smith, B. E., Wilson, T., Brown, A., Bordoni, A., Willis, M., Khan, S. A., RoviraNavarro, M., Dalziel, I., Smalley Jr., R., Kendrick, E., Konfal, S., Caccamise II, D. J., Aster, R. C., Nyblade, A., and Wiens, D. A.: Observed rapid bedrock uplift in Amundsen Sea Embayment promotes icesheet stability, Science, 360, 1335–1339, 2018. a
Bentley, C. R. and Wahr, J. M.: Satellite gravity and the mass balance of the Antarctic ice sheet, J. Glaciol., 44, 207–203, 1998. a
Bindschadler, R. A., Nowicki, S. M. J., AbeOuchi, A., Aschwanden, A., Choi, H., Fastook, J., Granzow, G., Greve, R., Gutowski, G., Herzfeld, U., Jackson, C., Johnson, J., Khroulev, C., Levermann, A., Lipscomb, W. H., Martin, M. A., Morlighem, M., Parizek, B. R., Pollard, D., Price, S. F., Ren, D., Saito, F., Sato, T., Seddik, H., Seroussi, H., Takahashi, K., Walker, T., R., and Wang, W. L.: Icesheet model sensitivities to environmental forcing and their use in projecting future sealevel (The SeaRISE Project), J. Glaciol., 59, 195–224, 2013. a, b, c
Bondzio, J. H., Morlighem, M., Seroussi, H., Kleiner, T., Ruckamp, M., Mouginot, J., Moon, T., Larour, E. Y., and Humbert, A.: The mechanisms behind Jakobshavn Isbrae's acceleration and mass loss: A 3D thermomechanical model study, Geophys. Res. Lett., 44, 6252–6260, 2018. a
Bueler, E., Lingle, C. S., and Brown, J.: Fast computation of a viscoelastic deformable Earth model for icesheet simulations, Ann. Glaciol., 46, 97–105, 2007. a
Caron, L., Ivins, E. R., Larour, E., Adhikari, S., Nilsson, J., and Blewitt, G.: GIA model statistics for GRACE hydrology, cryosphere and ocean science, Geophys. Res. Lett., 45, 2203–2212, 2018. a, b
de Boer, B., Stocchi, P., and van de Wal, R. S. W.: A fully coupled 3D icesheet–sealevel model: algorithm and applications, Geosci. Model Dev., 7, 2141–2156, https://doi.org/10.5194/gmd721412014, 2014. a
Durand, G., Gagliardini, O., de Fleurian, B., Zwinger, T., and Le Meur, E.: Marine ice sheet dynamics: Hysteresis and neutral equilibrium, Geophys. Res. Lett., 114, F03009, https://doi.org/10.1029/2008JF001170, 2009. a
Farrell, W. E. and Clark, J. A.: On postglacial sea level, Geophys. J. Roy. Astr. S., 46, 647–667, 1976. a, b
Favier, L., Pattyn, F., Berger, S., and Drews, R.: Dynamic influence of pinning points on marine icesheet stability: a numerical study in Dronning Maud Land, East Antarctica, The Cryosphere, 10, 2623–2635, https://doi.org/10.5194/tc1026232016, 2016. a
Goelzer, H., Coulon, V., Pattyn, F., de Boer, B., and van de Wal, R.: Brief communication: On calculating the sealevel contribution in marine icesheet models , The Cryosphere, 14, 833–840, https://doi.org/10.5194/tc148332020, 2020. a, b, c
Gomez, N., Mitrovica, J. X., Huybers, P., and Clark, P. U.: Sea level change as a stabilizing influence on marine ice sheets, Nat. Geosci., 3, 850–853, 2010. a, b
Gomez, N., Pollard, D., and Mitrovica, J. X.: A 3D coupled ice sheet – sea level model applied to Antarctica through the last 40 ky, Earth Planet. Sc. Lett., 384, 88–99, 2013. a, b
Gomez, N., Pollard, D., and Holland, D.: Sealevel feedback lowers projections of future Antarctic IceSheet mass loss, Nat. Comm., 6, 8798, https://doi.org/10.1038/ncomms9798, 2015. a
Gregory, J. M., Griffies, S. M., Hughes, C. W., Lowe, J. A., Church, J. A., Fukumori, I., Gomez, N., Kopp, R. E., Landerer, F., Le Cozannet, G., Ponte, R. M., Stammer, D., Tamisiea, M. E., and van de Wal, R. S. W.: Concepts and Terminology for Sea Level: Mean, Variability and Change, Both Local and Global, Surv. Geophys., 40, 1251–1289, 2019. a, b
Hillenbrand, C.D., Smith, J. A., Hodell, D. A., Greaves, M., Poole, C. R., Kender, S., Williams, M., Joest Andersen, T., Jernas, P. E., Elderfield, H., Klages, J. P., Roberts, S. J., Gohl, K., Larter, R. D., and Kuhn, G.: West Antarctic Ice Sheet retreat driven by Holocene warm water incursions, Nature, 547, 43–48, https://doi.org/10.1038/nature22995, 2017. a
Hogg, A. E., Shepherd, A., Gilbert, L., Muir, A., and Drinkwater, M. R.: Mapping ice sheet grounding lines with CryoSat2, Adv. Space Res., 62, 1191–1202, 2018. a
Hutter, K.: Theoretical Glaciology, Reidel Publishing Co., Dordrecht, Netherlands, pp. 510, 1983. a
Johnston, P.: The effect of spatially nonuniform water loads on predictions of sea level change, Geophys. J. Int., 114, 615–634, 1993 a
Jones, R. S., Mackintosh, A. N., Norton, K. P., Golledge, N. R., Fogwill, C. J., Kubik, P. W., Christl, M., and Greenwood, S. L.: Rapid Holocene thinning of an East Antarctic outlet glacier driven by marine ice sheet instability, Nat. Commun., 6, 8910, https://doi.org/10.1038/ncomms9910, 2015. a
Kingslake, J., Scherer, R. P., Albrecht, T., Coenen, J., Powell, R. D., Reese, R., Stansell, N. D., Tulaczyk, S., Wearing, M. G., and Whitehouse, P. L.: Extensive retreat and readvance of the West Antarctic Ice Sheet during the Holocene, Nature, 558, 430–434, 2018. a
Konrad, H., Sasgen, I., Pollard, D., and Klemann, V.: Potential of the solidEarth response for limiting longterm West Antarctic Ice Sheet retreat in a warming climate, Earth Planet. Sc. Lett., 432, 254–264, 2015 a, b
Lambeck, K., Purcell, A., Johnston, P. J., Nakada, M., and Yokoyama, Y.: Waterload definition in the glaciohydroisostatic sealevel equation, Quaternary Sci. Rev., 22, 309–318, 2003. a
Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res., 117, F01022, https://doi.org/10.1029/2011JF002140, 2012. a
Larour, E., Seroussi, H., Adhikari, S., Ivins, E. R., Caron, L., Morlighem, M., and Schlegel, N.: Slowdown in Antarctic mass loss from solid Earth and sealevel feedbacks, Science, 364, eaav7908, https://doi.org/10.1126/science.aav7908, 2019. a, b, c, d, e, f, g, h, i
Le Meur, E. and Huybrechts, P.: A comparison of different ways of dealing with isostasy: examples from modeling the Antarctic Ice Sheet during the last glacial cycle, Ann. Glaciol., 23, 309–317, 1996. a, b
Lingle, C. and Clark, J.: A numerical model of interactions between a marine ice sheet and the solid earth: application to a west Antarctic ice stream, J. Geophys. Res., 90, 1100–1114, 1985. a, b
Matsuoka, K., Hindmarsh, R. C. A., Moholdt, G., Bentley, M. J., Pritchard, H. D., Brown, J., Conway, H., Drews, R., Durand, G., Goldberg, D., Hattermann, T., Kingslake, J., Lenaerts, J. T. M., Martin, C., Mulvaney, R., Nicholls, K., Pattyn, F., Ross, N., Scambos, T., and Whitehouse, P. L.: Antarctic ice rises and rumples: Their properties and significance for icesheet dynamics and evolution, EarthSci. Rev., 150, 724–745, 2015. a
Milillo, P., Rignot, E., Rizzoli, P., Scheuchl, B., Mouginot, J., BuesoBello, J., and PratsIraola, P.: Heterogeneous retreat and ice melt of Thwaites Glacier, West Antarctica, Sci. Adv., 10, eaau3433, https://doi.org/10.1126/sciadv.aau3433, 2019. a, b
Milne, G. A.: Refining models of the glacial isostatic adjustment process, PhD Thesis, Physics Dept., University of Toronto, Canada, 1998. a, b
Milne, G. A. and Mitrovica, J.: Postglacial sealevel change on a rotating Earth, Geophys. J. Int., 133, 1–19, 1998. a, b
Mitrovica, J. X. and Milne, G. A.: On postglacial sea level: I. General theory, Geophys. J. Int., 154, 253–267, 2003. a
Nakayama, Y., Menemenlis, D., Zhang, H., Schodlok, M., and Rignot, E.: Origin of Circumpolar Deep Water intruding onto the Amundsen and Bellingshausen Sea continental shelves, Nat. Commun., 9, 3403, https://doi.org/10.1038/s41467018058131, 2018. a
Nowicki, S. M. and Wingham, D. J.: Conditions for a steady ice sheetice shelf junction, Earth Planet. Sc. Lett., 265, 246–255, 2008. a
Nowicki, S. M. J., Payne, A., Larour, E., Seroussi, H., Goelzer, H., Lipscomb, W., Gregory, J., AbeOuchi, A., and Shepherd, A.: Ice Sheet Model Intercomparison Project (ISMIP6) contribution to CMIP6, Geosci. Model Dev., 9, 4521–4545, https://doi.org/10.5194/gmd945212016, 2016. a
Pattyn, F.: Sealevel response to melting of Antarctic ice shelves on multicentennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISh v1.0), The Cryosphere, 11, 1851–1878, https://doi.org/10.5194/tc1118512017, 2017. a, b
Peltier, W. R., Argus, D. F., and Drummond, R.: Space geodesy constrains ice age terminal deglaciation: The global ICE6G_C (VM5a) model, Geophys. Res. Lett., 122, 450–487, 2015. a
Rignot, E., Mouginot, J., and Scheuchl, B.: Antarctic grounding line mapping from differential satellite radar interferometry, Geophys. Res. Lett., 38, L10504, https://doi.org/10.1029/2011GL047109, 2011. a
Sayag, R. and Grae Worster, M.: Elastic dynamics and tidal migration of grounding lines modify subglacial lubrication and melting, Geophys. Res. Lett., 40, 5877–5881, 2013. a
Schoof, C.: Ice sheet grounding line dynamics: Steady states, stability, and hysteresis, J. Geophys. Res., 112, F03S28, https://doi.org/10.1029/2006JF000664, 2007. a
Schoof, C.: Marine ice sheet stability, J. Fluid Mech., 698, 62–72, 2012. a
Sergienko, O. V. and Wingham, D. J.: Grounding line stability in a regime of low driving and basal stresses, J. Glaciol., 65, 833–849, 2019. a
Seroussi, H. and Morlighem, M.: Representation of basal melting at the grounding line in ice flow models, The Cryosphere, 12, 3085–3096, https://doi.org/10.5194/tc1230852018, 2018. a
Seroussi, H., Morlighem, M., Larour, E., Rignot, E., and Khazendar, A.: Hydrostatic grounding line parameterization in ice sheet models, The Cryosphere, 8, 2075–2087, https://doi.org/10.5194/tc820752014, 2014. a
Seroussi, H., Nakayama, Y., Larour, E., Menemenlis, D., Morlighem, M., Rignot, E., and Khazendar, A.: Continued retreat of Thwaites Glacier, West Antarctica, controlled by bed topography and ocean circulation, Geophys. Res. Lett., 44, 6191–6199, 2017. a, b
Tamisiea, M. E.: Ongoing glacial isostatic contributions to observations of sea level change, Geophys. J. Int., 186, 1036–1044, 2011. a
WCRP Global Sea Level Budget Group: Global sealevel budget 1993–present, Earth Syst. Sci. Data, 10, 1551–1590, https://doi.org/10.5194/essd1015512018, 2018. a
Whitehouse, P. L., Bentley, M. J., Vieli, A., Jamieson, S. S. R., Hein, A. S., and Sugden, D. E.: Controls on Last Glacial Maximum ice extent in the Weddell Sea embayment, Antarctica, J. Geophys. Res., 122, 371–397, 2017. a, b
 Abstract
 Copyright statement
 Introduction
 Land, ocean and ice domains and their interfaces
 Sealevel contribution from an ice sheet
 Sealevel change and mass conservation in the Earth system
 Conclusions
 Appendix A: Interpretation of ΔH_{F}(ω,Δt) and ΔH_{S}(ω,Δt)
 Appendix B: Notation
 Code and data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Copyright statement
 Introduction
 Land, ocean and ice domains and their interfaces
 Sealevel contribution from an ice sheet
 Sealevel change and mass conservation in the Earth system
 Conclusions
 Appendix A: Interpretation of ΔH_{F}(ω,Δt) and ΔH_{S}(ω,Δt)
 Appendix B: Notation
 Code and data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References