A mass conserving formalism for ice sheet, solid Earth and sea level interaction

Polar ice sheets are important components of any Earth System model. As the domains of land, ocean, and ice sheet change, they must be consistently defined 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 consistent and mass conserving descriptions of evolving land and ocean domains, grounded and floating ice masks, coastlines and grounding lines, and bedrock and geoid height as viewed from space. Here we present a geometric description of an evolving ice sheet margin and its 5 relations to sea level change, the position and loading of the solid Earth and include the ice shelves and adjacent ocean mass. We generalize the formulation so that it is applied to arbitrarily distributed ice, bedrock and adjacent ocean, and their interactive evolution. The formalism simplifies computational strategies that seek to conserve mass in Earth System models.

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. Defining geometry is a first order step for proper construction of models for tracking grounding line migration (e.g., Gudmundsson et al., 2012) and in this paper we propose a simple set of definitions that the field may find useful for future debate and reconciliation.
Much of the basic configuration setup for mechanical analysis of ice sheet evolution at the ice-bedrock-ocean interfaces has 5 been treated with great rigor, for example, in Chapter 3 of Hutter (1983). Similar setups are also familiar in the development of glacial isostatic adjustment (GIA) theory for ice, sea level and bedrock evolution following the Last Glacial Maximum (LGM) with migrating coast lines (e.g., Milne, 1998;Lambeck et al., 2003;Mitrovica and Milne, 2003). We refer to these as traditional configurations. Modern satellite techniques have allowed us to gain knowledge of both the present-day locations and migration rates of the grounding line (e.g., Rignot et al., 2011;Milillo et al., 2019). However, both observation and numerical simula-10 tion of subtle change in grounding line positions are complicated by the presence of complex geometric features including floating ice shelves, ice rises and rumples, and retrograde bedrock slopes. These features complicate the required geometrical simplifications used in traditional theory for ice-bedrock-ocean interface changes, especially if the interactions are two-way (e.g., Lingle and Clark, 1985;Gomez et al., 2013;de Boer et al., 2014). For example, estimating sea level contribution of an evolving ice sheet in such cases becomes non-trivial (Goelzer et al., 2019). 15 In the following, we begin by presenting a generalized definition of land, ocean and ice sheet domains and that of coastline and grounding line positions (Section 2). We consider a distributed system of ice domains comprised of glaciers and ice sheets of arbitrary geometric configurations 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. We then explore the implication of this new geometrical setup for sea level and solid Earth loading studies (Section 3). Finally, in Section 4, we 20 summarize the key conclusions of the study.

Land, ocean and ice sheet domains and their boundaries
To begin our discussion of a proposed geometrical setup, we consider a spherical planet whose surface is divided into land and ocean. Distributed ice domains including glaciers and ice sheets exist on both land and ocean ( Figure 1). We define B(ω, t) and S(ω, t) to be the land/bedrock/seafloor and sea surface elevations, respectively, measured relative to the same datum, 25 preferably consistent with the International Terrestrial Reference Frame (e.g., Altamimi et al., 2016). Here ω denotes 2-D spatial coordinates on the planetary surface at time t. Depending upon the spatial scale, we interchangeably use ω to represent for geographic coordinates (θ, φ) or Cartesian coordinates (x, y), assuming that an appropriate coordinate transformation is applied. We define H(ω, t) as the thickness of ice bounded between its surface and the base. When the base of ice is in contact with the underlying bedrock, ice is grounded. Ice may float on subglacial or proglacial lakes or on the ocean, with its base 30 above the bedrock. Our focus here is on marine portions of the ice sheet, with B(ω, t) < S(ω, t).
It is important to note that S(ω, t) may be highly variable both in space and time due to relatively short-term dynamic processes such as tides, wind stress, atmospheric pressure variability, and associated ocean circulation. At timescales on the Lakes are considered as part of the land. Ice can have multiple domains, shown here with blue sheds. The land-ocean 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 the ice sheet, we assume that all of ice on land (gridded portions of blue sheds) is grounded. Consequently, flotation of ice on subglacial and proglacial environments are not considered in this study. order of 20 years or longer, changing interactions between sea level and ice sheet may become important (e.g., Hillenbrand et al., 2017;Larour et al., 2019). Hence, any long-term change in S(ω, t) caused, for example, by a sustained water and heat exchange between the land ice and the ocean is of central interest to our formulation. Here, S(ω, t) strictly refers to this quasistatic component of sea surface that is free from ocean dynamic signal and high-frequency signals of waves or meteorologically driven fluctuations. We hereafter refer S(ω, t) to as "sea level" for brevity. Sea level as defined above represents an equipotential 5 surface whose spatial pattern mimics the geoid (Tamisiea, 2011;Gregory et al., 2019). Under this definition, we may include areas that are non-oceanic, including the interior of marine ice sheets. Figure 2a shows a cross-sectional view of land, ocean, and the grounded and floating portions of ice sheet domains. We introduce a function F (ω, t) that is used to define these domains in terms of the principles of hydrostatic equilibrium: 10 where ρ o and ρ i are the densities of ocean water and ice, respectively. Note that, by definition, R(ω, t) = S(ω, t) − B(ω, t) where R(ω, t) is sea surface relative to the seafloor, termed the "relative sea level" (Gregory et al., 2019). By design, the position with F (ω, t) = 0 represents a boundary between the land and the ocean, provided that ω is in direct contact with the open ocean. This land-ocean boundary is generally known as a coastline. It follows from equation (1) that ice thickness at the coastline, H C (ω, t), is given by: and that the coastline is free of ice where B(ω, t) = S(ω, t). Since ice thickness cannot be negative, 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 then is replaced by the term grounding line.
Given the definition of coastline and grounding line, we may define the ocean domain, O(ω, t), as follows: = 0 otherwise, except when ω is precisely a coastline or a grounding line. (3)

5
The complementary land domain is simply given by L(ω, t) = 1 − O(ω, t). Note that neither O(ω, t) nor L(ω, t) is defined at the coastline or grounding line. The concept of "connectivity" is introduced because the criterion F (ω, t) < 0 in and of itself is not always sufficient to define the ocean domain. One obvious example where ω may not belong to ocean in spite of having F (ω, t) < 0 is a deep continental trough with bathymetry well below sea level. Unless this trough or any such locations is physically connected to the open ocean via oceanic water, we consider this to be part of land rather than the ocean.

10
This assumption ensures an interconnected system of Earth's oceanic waters, termed global ocean, as has been considered traditionally in physical oceanography and sea level studies.
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 sub-domains: where I i (ω, t) represents the i-th ice domain. Individual ice sheets and glaciers can be thought of individual ice domains. As defined above, the grounding line within a given ice domain is where F (ω, t) = 0, provided again that ω is in direct contact with the open ocean. Using equations (3) and (4), 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). Note that F (ω, t) ≤ 0 can exist even within the grounded ice mask. A condition where this may occur is when the interior of a marine ice sheet resides on a deep trough (see Figure 2a). If a column of ice is shielded 20 by surrounding bathymetric highs from the open ocean, we expect it to remain grounded. The ice sheet domain as defined by equation (4) can accommodate complexities such as a pinning point that can modulate marine ice sheet instability on retrograde slopes (e.g., Matsuoka et al., 2015;Whitehouse et al., 2017). This is where the new setup diverges from traditional theory. The employed assumption in our description, however, limits us from capturing the floating ice on subglacial and proglacial lakes that are not part of the global ocean (see Figure 1). We believe that the localized processes of ice-lake interactions are of 25 secondary importance, at least, for the purpose of capturing large-scale interplay between the continental ice sheets, solid Earth and sea level in the first generation of Earth System models.
Our definition of coastline and grounding line, and hence that of land, ocean and ice domains, 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 sheet driven global and regional sea level change. Although a distributed system of ice domains is an integral part of the 30 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.
In order to estimate sea level contribution of an ice sheet, we define a flotation height of ice so that the ice thickness in excess of H 0 (ω, t) represents the so-called "height above flotation", H F (ω, t). Mathematically, Here we invoke a grounded ice mask G(ω, t) to ensure no net sea level contribution from the floating ice shelf. Clearly, for the grounded ice sheet that rests on the bedrock whose elevation is at or above sea level. For grounded portions of the marine ice sheet, H F (ω, t) < H(ω, t) and it can be negative (see regime 4 in Figure 2a), implying that this portion can take up ocean water when it is connected to the open ocean and hence contribute to sea level inversely.
The evolving ice sheet geometry is usually described in terms of evolving ice thickness and ice sheet domain. Indeed, 10 prognostic simulations of ice sheet models track the transport of mass in terms of equivalent ice thickness distribution. This redistribution of ice mass and associated relative sea level change induce a solid Earth response in terms of its gravitation, rotation and viscoelastic deformation. This, in turn, modulates the bedrock topography as well as the geoid. We may therefore describe the evolving ice sheet geometry in terms of ice thickness, bedrock elevation, and sea level (Figure 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 ice sheet 15 geometry at time t + ∆t, equation (6) gives , 0 is given by equation (5). Similarly, we may rewrite equation (1) for F (ω, t + ∆t) and define the new ocean domain O(ω, t + ∆t), land domain L(ω, t + ∆t), ice sheet domain I(ω, t + ∆t), and grounded ice mask G(ω, t + ∆t) as described in Section 2. Note that ∆B(ω, ∆t) and ∆S(ω, ∆t) may have components that are forced by external processes such as farfield ice melting or tectonics. 20 In what follows, we assume that the net change in ice sheet mass directly affects the ocean mass. Quantifying the fraction of ice mass change that contributes to sea level is a non-trivial problem. It cannot be approached analytically and is beyond the scope of this study. Despite the assumption, not all of ∆H(ω, ∆t) contributes to sea level change all the time. On the other hand, in response to the externally forced bedrock and sea level, the ice sheet may still contribute to sea level change even when ∆H(ω, ∆t) = 0. Traditionally, change in ice thickness above flotation, i.e. ∆H F (ω, ∆t) = H F (ω, t + ∆t) − H F (ω, t), is used 25 to calculate sea level contribution of an ice sheet (e.g., Bindschadler et al., 2013). As we show below this simplistic approach yields some error, particularly when evolving bedrock and sea level are considered (Larour et al., 2019;Goelzer et al., 2019).
We define ∆H S (ω, ∆t) be a portion of ice thickness that contributes to sea level over the period ∆t. The following relationship holds for generalized ice geometries and bedrock and/or sea level forcings: ∆H S (ω, ∆t) = ∆H(ω, ∆t) L(ω, t) L(ω, t + ∆t) + ∆H F (ω, ∆t) |∆L(ω, ∆t)| + 1 − ρ w ρ o H 0 (ω, t + ∆t) ∆L(ω, ∆t), (7) 30 where ρ w is the fresh water density. The change in land domain is simply given by ∆L(ω, ∆t) = L(ω, t + ∆t) − L(ω, t). Since ∆H F (ω, ∆t), unlike H 0 (ω, t + ∆t), can track the direction of margin migration (i.e., advance or retreat), absolute value of and H0(ω, t + ∆t) in regime 2 (see equation 7), and is zero in regime 3. We zoom in around the grounding line to assess ∆HS(ω, ∆t) in different scenarios: (c) when ice thickness changes but relative sea level does not, typically assumed in stand-alone ice sheet models; (d) when (externally forced) sea level changes but ice thickness does not; and (e) when both ice thickness and (ice driven or externally forced) relative sea level change. In b and d, only sea level, not the bedrock, is changed for simplicity. Sketches are not to scale.
∆L(ω, ∆t) is used in the second term on the right-hand side. This method requires bookkeeping the global land domain. This is crucial for considering distributed ice-ocean 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, often it may be possible to replace all L's appearing in equation (7) by corresponding G's, the grounded ice masks. However, 5 in an exceptional case with H 0 (ω, t + ∆t) > 0 holding in the areas of on-land ice margin migration, using G's rather than L's incorrectly predict a non-zero contribution from the last term in equation (7). Figure 2 shows a few schematics of evolving ice sheet geometries and their sea level contributions. We consider all plausible cases by combining scenarios of ice thickness and relative sea level changes. In reference to these sketches and equation (7), we outline three distinct regimes of an evolving ice sheet and their sea level contribution over the period ∆t: 1. Where ice remains grounded on the bedrock at the elevations above, at, or below corresponding sea level at both times t and t + ∆t, all of ∆H(ω, ∆t) contributes to sea level (first term on the right-hand side). It turns out that ∆H(ω, ∆t) = 5 ∆H F (ω, ∆t) only in marine portions of the grounded ice and only when evolving bedrock and sea level are considered.
Goelzer et al. (2019) present a method to backtrack ∆H(ω, ∆t) from ∆H F (ω, ∆t) in such situations, assuming that ∆B(ω, ∆t) and ∆S(ω, ∆t) are known. Over the timescale of a few decades or shorter, it may be that relative sea level does not evolve significantly, and in these cases ∆H F (ω, ∆t) ≈ ∆H(ω, ∆t). Stand-alone ice sheet models typically inherit this assumption even though simulation timescales can be on the order of centuries (e.g., Bindschadler et al., 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 recently deglaciated areas. In both cases, all of ∆H(ω, ∆t) contributes to 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 15 yield nonzero ∆H F (ω, ∆t) and possibly modulate the ice flow dynamics.
2. Where ice has transitioned from the grounded to floating state or vice versa over the period ∆t, sea level contribution from this regime depends on both ∆H F (ω, ∆t) and H 0 (ω, t + ∆t) (the last two terms on the right-hand side). In the absence of externally forced bedrock and sea level, it follows that |∆H F (ω, ∆t)| ≤ |∆H(ω, ∆t)| in this regime.
The last term in the equation accounts for the fact that fresh water density evolves during the accretion and ablation of 20 ice, whereas the average ocean water density in the vicinity of the grounding line acts to determine the flotation height.
For typical values of ρ w and ρ o , we find that about 2.5% of H 0 (ω, t+∆t) contributes to sea level change. This particular contribution has been considered by Goelzer et al. (2019).
This is the only regime in which responses to the externally forced ∆B(ω, ∆t) and/or ∆S(ω, ∆t) contribute to sea level change, even when ∆H(ω, ∆t) = 0 (see Figure 2d).

25
3. Where ice is floating at both times t and t + ∆t, there is no net direct contribution from this regime of the ice sheet to sea level change. However, change in ice shelf thickness can greatly affect the interior ice sheet dynamics via modulation of buttressing force (e.g., Gudmundsson et al., 2019).
In Figure 3, we present a quantitative assessment of ∆H S (ω, ∆t) and ∆H F (ω, ∆t) over the next 350 years for Pine Island and Thwaites glaciers. Notice the systematic error associated with a customary approach of using ∆H F (ω, ∆t) to quantify sea 30 level contribution from an ice sheet (Figure 3c). The discrepancy in regime 1 is due to the evolving bedrock and sea level, and that in regime 2 is due to the missing fraction of newly grounded or floating ice. These results are based on the recent work Estimated ice thickness change that directly contributes to sea level, given by equation (7). (c) Error associated with a conventional method that assumes that sea level contribution of an ice sheet is given by ∆HF (ω, ∆t). In the latter two panels, we point out three regimes of the ice sheet as defined in Section 3 (also see Figures 2b-e) in which the following holds true: ∆HF (ω, ∆t) = ∆HS(ω, ∆t) = ∆H(ω, ∆t) in marine portions of regime 1; ∆HF (ω, ∆t) = ∆HS(ω, ∆t) = ∆H(ω, ∆t) in regime 2; and ∆HF (ω, ∆t) = ∆HS(ω, ∆t) = 0 in regime 3. In this particular example, the conventional method underpredicts sea level contribution from the ice sheet by the amount that corresponds to panel c.
of Larour et al. (2019), who provide consistent solutions of evolving H(ω, t), B(ω, t) and S(ω, t) for Antarctica 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. The effect of farfield ice mass change, extrapolated into the next 500 years based on the contemporary trend, on the evolution of bedrock and sea level in Antarctica is also accounted for.

5
The barystatic sea level change, defined as the ocean mass related global mean sea level change (Gregory et al., 2019), due to ice mass change over the period ∆t is given by The numerator in the equation owes to the net change in ice mass over ∆t, and the integral in denominator represents the ocean surface area at time t + ∆t. Assume that an ice sheet collapses instantaneously and that all of melt water contributes to sea 10 level. Resulting ∆R I C (∆t) represents the "potential sea level" of the ice sheet at time t, and it can be readily derived by setting ∆H(ω, ∆t) = −H(ω, t) and G(ω, t + ∆t) = 0 in the limit of ∆t → 0 in equation (7).
Ocean surface area changes as grounding line and coastline positions migrate. Such a migration over ∆t is controlled by a number of processes that contribute to ∆R(ω, ∆t). A comprehensive synopsis of these processes can be stated as follows: ∆R(ω, ∆t) = ∆R I C (ω, ∆t) + ∆R L C (ω, ∆t) + ∆R I P (ω, ∆t) + ∆R L P (ω, ∆t) + ∆R D (ω, ∆t) + ∆R V (ω, ∆t). The first four terms on the right-side of the equation represent the processes that contribute to sea level by changing the ocean mass (termed the barystatic components), while ∆R D (ω, ∆t) represents the density related steric sea level change and ∆R V (ω, ∆t) represents a component due to vertical land motion that is not captured by the barystatic processes. We use the superscript I to refer to the ice sheet under consideration and superscript L to the other parts of the land, including farfield ice sheets, glaciers and hydrological basins, that contribute to barystatic sea level change. When these sources of freshwater Assuming that the contemporaneous period ∆t is on the order of 10 years, we may interpret ∆R(ω, ∆t) as ongoing change in 10 relative sea level monitored by satellite gravimetry and altimetry (e.g., WCRP Global Sea Level Budget Group, 2018). Along the same lines, we may interpret ∆R I C (ω, ∆t) + ∆R L C (ω, ∆t) and ∆R I P (ω, ∆t) + ∆R L P (ω, ∆t) as ongoing sea level change driven by contemporary global 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 ice sheet evolves, it not only contributes to sea level but also loads the underlying solid Earth. For the period [t, t + ∆t], 15 we may define a mass conserving field ∆M (ω, ∆t) that describes the net change in mass per unit area on the solid Earth surface as follows: ∆M (ω, ∆t) = ρ i ∆H S (ω, ∆t) + ρ w ∆R I C (ω, ∆t) O(ω, t + ∆t).
Here ∆H S (ω, ∆t) is given by equation (7) and ∆R I C (ω, ∆t) represents the associated relative sea level change, whose spatial pattern is dictated by the perturbation in Earth's gravitational and rotational potentials and associated viscoelastic deformation 20 of the solid Earth (Farrell and Clark, 1976;Milne and Mitrovica, 1998). This component of sea level is precisely the same as the first term on the right-side of equation (9), whose ocean averaged value is given by equation (8). Because sea level is defined globally, including on land, we ensure mass conservation by use of an ocean mask in equation (10).

Conclusions
We have divided the Earth's surface into complementary domains of land and ocean, which are separated by coast lines. While 25 there may be multiple land domains, we maintain a single global ocean as in traditional studies of physical oceanography and sea level. Distributed bodies of ice intersect land and 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, descriptions given in this paper handle the complex geometries of both the coastlines and evolving grounding lines, complementary to those of farfield land, ocean and ice domains and their respective overall evolutionary history. 30 The main importance of this paper is that a unified method is outlined to determine the exact fraction of ice thickness change that contributes to sea level change, ∆H S (ω, ∆t), over the period ∆t. Along with its obvious application to estimate global mean sea level, this field 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 fine-scale complexities in geometry. In the most simplified case when bedrock and sea level do not evolve, our method reduces to the traditional approach that assumes that change in ice height above flotation, ∆H F (ω, ∆t), directly contributes to sea level. We recommend that the ice sheet modeling community considers ∆H S (ω, ∆t) rather than ∆H F (ω, ∆t) as a metric to quantify the sea level contribution of evolving ice sheets. This is 5 especially appropriate for model analyses 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).