A kinematic formalism for tracking ice-ocean mass exchange on the Earth’s surface and estimating sea-level change

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 ice-sheet dynamics, solid-Earth deformation, and sea-level adjustment requires both geodetically consistent and mass conserving 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 5 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 sea-level 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 time scales. 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. 10 Copyright statement. © 2020 California Institute of Technology. Government sponsorship acknowledged.


Introduction
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 sea-level 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 ice-sheet 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 kilometer-scale 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-bedrockocean 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 level-set method, which has been previously applied for tracking grounding lines (e.g., Seroussi et al., 2014) and calving-front 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 ice-sheet dynamics via the modulation of bedrock slope, grounding-line 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 observation-based (e.g., Barletta et al., 2018;Kingslake et al., 2018) and the model-based studies (e.g., Lingle and Clark, 1985;Gomez et al., 2010Gomez et al., , 2015Adhikari et al., 2014;Konrad et al., 2015;Larour et al., 2019). The inclusion of evolving bedrock and sea level in a dynamical ice-sheet model, however, requires a modification to the common method of estimating sea-level 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 sea-level 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 stand-alone ice-sheet models and those that account for isostatic bedrock adjustment (e.g., Le Meur and Huybrechts, 1996;Pattyn, 2017) or a self-consistent 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 sea-level 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 sea-level change and mass conservation in the Earth system. Finally, in Sect. 5, we summarize the key conclusions.

Land, ocean and ice domains and their interfaces
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 2-D 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 high-frequency 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 dif- Figure 1. Conceptual depiction of land, ocean and ice domains in the Earth system. Gridded areas represent land and the rest represents the ocean. Lakes are considered part of the land. Ice can have multiple domains, shown here with blue shading. The landocean boundary is generally defined as the coastline, which is called grounding line when it is part of the ice domain. Because our focus is on grounding line migration in marine portions of an icesheet/shelf system, we assume that all of ice on land (gridded portions of blue shading) is grounded. Consequently, floatation of ice on subglacial and proglacial environments is not considered in this study.
ference, 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 so-called barystatic sea-level change, which is the globalmean sea-level (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 self-consistent 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).

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 (ω, t) = 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 (ω, t) = 0 satisfies the hydrostatic equilibrium between ice and the oceanic water in the marine sectors where B(ω, t) < S(ω, 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 (ω, t) ≤ 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(ω, t) = S(ω, 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 zero-level set may consist of several simple curves, T i (F ), that divide the planetary surface into several nonoverlapping regions, i (F ). Let − i (F ) 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 sea-level studies, only the largest amongst − i (F ) forms the ocean domain. Smaller − i (F ), 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 (ω, t) < 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 − S (F ) 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 zero-level set represents the land-ocean boundaries. Here is a positive number to ensure F(ω, t) > 0 at ω ∈ T S (F ).
The land-ocean boundaries are generally known as coastlines. Given the definition of the level-set function (Eq. 3), coastlines are free of ice where B(ω, t) = S(ω, t). No coastline exists with B(ω, t) > S(ω, t). Only in the marine sectors where B(ω, t) < S(ω, t) does a coastline have finite ice thickness and is then replaced by the term grounding line.

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 L( 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 O(ω, t) nor L(ω, t) is defined at the coastlines or grounding lines. These interfaces rather form their own level set, T (F), 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 Ice-sheet and Sea-level 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 I(ω, t) to be a globally distributed system of ice domains, such that For many applications, it may be useful to decompose I(ω, t) into a number of subdomains: I(ω, t) = {I 1 , I 2 , . . ., I i , . . .}, where I 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 level-set T (F). Using Eqs. (5) and (6), we may define the grounded ice mask simply as G(ω, t) = I(ω, t)L(ω, t) and the floating ice mask as I(ω, t)O(ω, t). Equation (6) is a simple, and perhaps the most generic, definition for ice domains, which can accommodate any geometric features such as kilometer-scale 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 large-scale 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 ice-driven global and regional sea-level 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 far-field ice.

Sea-level contribution from an ice sheet
The estimation of the sea-level contribution from an evolving ice sheet, featuring marine-based 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.

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 so-called 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 land-ice thickness that can potentially contribute to sea-level change (Fig. 2a). It is clear from the equations that H F (ω, t) = H (ω, 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 F (ω, t) < H (ω, t) and, in fact, H F (ω, t) can take negative values. If the ice sheet were to disintegrate, for instance, regions with H F (ω, t) < 0 such as Sector D in Fig. 2a can take up ocean water and contribute to sea-level fall. The evolving ice-sheet geometry is usually described in terms of ice thickness and ice-sheet margins. Indeed, prognostic simulations of ice-sheet models track the transport of mass in terms of equivalent ice-thickness 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 far-field ice, or GIA, or tectonics. We may describe the evolving ice-sheet 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)  where H 0 (ω, t + t) = ρ o /ρ i max S(ω, t) + S(ω, t) − B(ω, t) − B(ω, t) , 0 is given by Eq. (7). Similarly, we may rewrite Eq. (1) for F (ω, t + t) and define the new ocean domain O(ω, t + t), land domain L(ω, t + t) and grounded ice domain G(ω, t + 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 H (ω, t) = 0 as we count floating ice in the ocean mass (see Appendix A). The stand-alone ice-sheet models evaluate the change in HAF in order to calculate the sea-level 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 oceanic-water volume (rather than the freshwater volume) by spatially integrating −[ρ i /ρ o ] H F (ω, t) and divide it by the ocean surface area to estimate the GMSL change. Apart from this water-density-related error, the HAF method in absence of evolving bedrock and MSL yields the correct estimates of the sea-level 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 stand-alone ice-sheet 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).

A new field for estimating sea-level 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 sea-level 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 sea-level 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 H (ω, t) = H F (ω, 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 ice-sheet margins. When ice margin advances over the period t, newly glaciated areas must satisfy H (ω, t) = 0 and H (ω, t) > 0. When it retreats, H (ω, t) = −H (ω, t) must hold in the recently deglaciated areas. In both cases, all of H (ω, t) contributes to the sea-level 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 ice-flow dynamics. While the effects of far-field 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 sea-level 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 | H F (ω, t)| < | H (ω, t)| (see Appendix A2). The change in ice thickness in excess of the change in HAF (right-side 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 H (ω, t) = 0. Specific examples are given in Appendix A2.
-Regime 3: where ice remains floating at both times t and t + t.
Since H F (ω, t) = 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 (right-side term in Eq. 12) contributes to the sea-level change (see Appendix A3).
Given the new field for estimating sea-level contribution (Eq. 10), we may readily calculate the GMSL change by spatially integrating −[ρ i /ρ w ] H S (ω, 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 H (ω, t) = −H (ω, t) and G(ω, t + t) = 0 in the limit of t → 0. Note that H (ω, t) and G(ω, t + t) are implicit via Eqs. (9), (11) and (12).

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 high-resolution dynamical ice-flow model (Larour et al., 2012) that is fully coupled with a global solid-Earth deformation and sea-level adjustment model (Adhikari et al., 2016) under the present-day surface climatology and a realistic sub-ice-shelf melting scenario (Seroussi et al., 2017). They also account for the effects of far-field ice-mass 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 space-gravimetry-based 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 off-line 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 ice-thickness change in excess of the change in HAF (see Appendix A2). We also show the time series of the total Antarctic ice-volume change that is attributable to the GMSL change ( Fig. 3d and e). We find that the new method predicts systematically larger sea-level 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 sea-level contribution than the HAF method by about 5 %. Note that the HAF method usually converts the ice-volume 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.

Sea-level change and mass conservation in the Earth system
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 land-surface 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 ice-sheet dynamics, solid-Earth deformation and sea-level adjustment may we find self-consistent 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, R(ω, t) = S(ω, t) − B(ω, 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 T (F), 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 space-based 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 sea-level 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 far-field ice and hydrological basins. When these sources of freshwater contribute to sea-level change over a contemporaneous period [t, t + 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 (−∞, t], and the induced viscous response of the solid Earth may still contribute to the RSL change over the period [t, t + 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 gravime- try and altimetry (WCRP Global Sea Level Budget Group, 2018). We may interpret R I C (ω, t) + R L C (ω, t) and R I P (ω, t) + R L P (ω, 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 sea-level adjustment. Because the GRD effect is applied to the entire column of the ocean water, only  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 longer-term (decadal or longer timescale) sea-level 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 "sea-level equation" on a self-gravitating, viscoelastically compressible, rotating Earth (Farrell and Clark, 1976;Milne and Mitrovica, 1998). To this end, we must consider a mass-conserving 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 R I C (ω, 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 sea-level equation, in essence, means that we load the solid Earth by the mass-conserving surface load (Eq. 14) and let its GRD response dictate the self-consistent patterns of RSL, MSL and bedrock changes as well as the new positions of coastlines and grounding lines. The ice-sheet 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 R I C (ω, 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).

Conclusions
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 far-field 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 global-mean sea-level 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 kilometer-scale 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 sea-level contribution from the Antarctic Ice Sheet. We recommend that the ice-sheet modeling community consider the proposed method as a metric to quantify the sea-level 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). Figure A1. Scenarios of ice-thickness and RSL change and sea-level contributions. (a) Since the column of ice remains grounded at both times and its thickness does not change over the period, the ice column does not contribute to GMSL change. The HAF method incorrectly predicts GMSL drop because of a positive value of H F in response to the imposed drop in RSL. (b) For the grounded ice to float in the case of fixed RSL, it must thin sufficiently, such that H < H F < 0 holds true, contributing to the GMSL rise. (c) Significant rise in RSL may cause the grounded ice to float even if its thickness does not change. As a result, the column of ice contributes to the GMSL rise. (d) Melting of the floating ice produces freshwater that occupies slightly larger volume than that of the ocean water that was replaced by the ice. This excess volume, the hatched portion of the freshwater column, causes the GMSL to rise. Sketches are not to scale.
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: R(ω, t) ≥ 0 and H (ω, t) > 0. In this case, the only condition for the floating ice to be grounded is through its sufficient thickening such that 0 < H F (ω, t) < H (ω, t). Both terms appearing in Eq. (A3) take positive values, causing the GMSL to fall.
-Case A2.5: R(ω, t) < 0 and H (ω, t) = 0. Here the condition that H (ω, t) = 0 < H F (ω, 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: R(ω, t) < 0 and H (ω, t) = 0. Only when H (ω, t) < 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 sea-level contribution does not depend on the evolving RSL. In this scenario, the ASL is given in terms of freshwater equivalent height by [ρ i /ρ w −ρ i /ρ o ] H (ω, t), whose ice equivalent height can be deduced from Eq. (12) by setting H F (ω, t) = 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. A The zero-level set of F that satisfies the hydrostatic equilibrium between ice and the oceanic water T The zero-level set of F that represents the coastlines and grounding lines t Time t Time period ω 2-D spatial coordinates, geographic (θ, φ) or Cartesian (x, y), on the planetary surface Code and data availability. The data and code used to produce Fig. 3 are available online at https://doi.org/10.7910/DVN/9LUJTD (Adhikari et al., 2020).