ISSM-SESAW v 1 . 0 : mesh-based computation of gravitationally consistent sea-level and geodetic signatures caused by cryosphere and climate driven mass change

A classical Green’s function approach for computing gravitationally consistent sea-level variations associated with mass redistribution on the earth’s surface employed in contemporary sea-level models naturally suits the spectral methods for numerical evaluation. The capability of these methods to resolve high wave number features such as small glaciers is limited by the need for large numbers of pixels and high-degree (associated Legendre) series truncation. Incorporating a spectral model into (components of) earth system models that generally operate on a mesh system also requires repetitive forward and inverse transforms. In order to overcome these limitations, we present a method that functions efficiently on an unstructured mesh, thus capturing the physics operating at kilometer scale yet capable of simulating geophysical observables that are inherently of global scale with minimal computational cost. The goal of the current version of this model is to provide high-resolution solidearth, gravitational, sea-level and rotational responses for earth system models operating in the domain of the earth’s outer fluid envelope on timescales less than about 1 century when viscous effects can largely be ignored over most of the globe. The model has numerous important geophysical applications. For example, we compute time-varying computations of global geodetic and sea-level signatures associated with recent ice-sheet changes that are derived from space gravimetry observations. We also demonstrate the capability of our model to simultaneously resolve kilometer-scale sources of the earth’s time-varying surface mass transport, derived from high-resolution modeling of polar ice sheets, and predict the corresponding local and global geodetic signatures.


Introduction
Earth system modeling of climate warming scenarios and their impact on society requires ever greater capacity to incorporate appropriate coupling of models that traditionally have operated in isolation from one another.One example is the necessity to couple the redistribution of earth surface mass and energy during secular and non-secular changes.The coupling of the major ice sheets to the earth's timevarying geoid was a main subject of Erich von Drygalski's PhD thesis (Drygalski, 1887), wherein, along with contemporary work of Woodward (1888), the ability of continental ice mass to attract ocean mass and alter sea level was first discussed and given theoretical treatment.However, the formal modern theory was only expounded almost a century later by Farrell and Clark (1976), who incorporated full accounting of elastic and viscous solid-earth deformation on a global scale.This theory is now fully incorporated into the literature for computing past and present-day sea-level variations (e.g., Wu and Peltier, 1983;Mitrovica and Peltier, 1991;Tushingham and Peltier, 1991;Mitrovica and Milne, 2003;Hagedoorn et al., 2007;Riva et al., 2010;Tamisiea, 2011) and for which contemporary software is available (Spada and Stocchi, 2007).
The importance of gravitational loading and self-attraction on earth system modeling is now demonstrated, for example, via coupling to ocean circulation (Kuhlmann et al., 2011) and ice-sheet models (e.g., Gomez et al., 2013;de Boer et al., 2014;Konrad et al., 2015) and analysis of the earth's rotational variability (Quinn et al., 2015).The recent focus of ice-sheet modeling is to incorporate more complicated physical processes that are demanded by new data from both space and terrestrial observations (e.g., Joughin et al., 2014;Schlegel et al., 2015;Pattyn, 2016).Simulation of changes in grounding line positions (e.g., Rignot et al., 2014) and velocity fields (e.g., Mouginot et al., 2014) of outlet glaciers requires kilometer-scale model resolution.Space gravimetric and ground-based crustal deformation data are also relevant to mass balance estimation over drainage basin to continental scales (e.g., Sasgen et al., 2013;Bouman et al., 2014) and should be carefully accounted for.Incorporating geopotential changes and solid-earth deformation into icesheet simulations, for example, allows grounding line positions to be influenced by changes in ocean-land contact position and the geometry of ice-shelf pinning points (e.g., Gomez et al., 2013;Adhikari et al., 2014;Konrad et al., 2015).Such incorporation implies global-scale simulation, and especially appropriate are solutions of the sea-level equation, as is often realized in global viscoelastic glacial isostatic adjustment (GIA) models (e.g., Clark, 1977;Hagedoorn et al., 2007).Our main goal here is to wed computations of global-scale solid-earth deformation and sea-level variation that are driven by kilometer-scale (high-order) mechanics of ice sheets on timescales less than 1 century when viscous effects can largely be ignored over most of the globe.
A major obstacle to efficiently coupling existing models has been their fundamentally different computational frameworks: 3-D ice-sheet models often operate on an unstructured mesh (e.g., Gagliardini and Zwinger, 2008;Larour et al., 2012), whereas self-gravitating sea-level models are mostly based on a pseudo-spectral method (e.g., Mitrovica and Peltier, 1991;Kendall et al., 2005) or a hybrid spectral/finite-element method (e.g., Martinec, 2000;Hagedoorn et al., 2007).While the spectral methods are now in widespread use (e.g., Spada and Stocchi, 2007;Whitehouse et al., 2007;Riva et al., 2010), the earliest developments for numerical solution of the sea-level equation involved discretization (e.g., Clark, 1977;Wu and Peltier, 1983;Tushingham and Peltier, 1991).However, the discrete methods were never fully examined for their potential computational advantages, nor were the generalizations to flexible adaptive meshing on the sphere examined and exploited.Models for ice-sheet coupling to ocean and solid earth require interfacing to variable adaptive meshing schemes with a hierarchy of scales.A systematic and efficient framework for capturing the kilometer-scale continuum physics of the grounding line interfaced to a globally resolved coastline and an open global ocean has not previously been provided.Here we present the new code, developed within the Jet Propulsion Laboratory's (JPL) massively parallelized and highly scalable Ice Sheet System Model (ISSM; Larour et al., 2012), to treat the full solid-earth deformation and sea-level coupling on a global scale, and yet retain the high-resolution and high-order capabilities of ice-flow models.This computational framework, which is termed ISSM's Solid Earth and Sea-level Adjustment Workbench (ISSM-SESAW), allows a straightforward and computationally less burdensome numerical approach to be realized.The main limitation of the current version is that the timescales must be on the order of a few decades for the elastic approximation to be valid.
In Sect.2, we briefly review the standard Green's function approach for solving the perturbation theory of relative sea level applied to an elastically compressible and density layered self-gravitating, rotating earth.In Sect.3, we provide our approach to evaluating key components of this theory on an anisotropic mesh and demonstrate its superiority (in terms of high-resolution capability, numerical accuracy, and computational efficiency) over contemporary pseudospectral methods.As example applications, in Sect.4, we produce computations of global geodetic and sea-level signatures associated with the recent evolution of polar ice sheets.The polar ice-sheet mass budget data are derived from space gravimetry observations, and hence are of relatively low resolution (on the order of 300 km).In order to demonstrate the high-resolution capability of our model, we provide in Sect. 5 sea-level fingerprints induced by high-resolution mass change of both polar ice sheets, as modeled by ISSM's core ice-flow capability.Finally, in Sect.6, we summarize key conclusions of this research and briefly outline its scope and limitations.

Theory of relative sea level
Redistribution of mass on the earth's surface caused by cryosphere and other climate driven phenomena, such as wind stress, ocean currents, and land water storage, perturbs the gravitational and rotational (centrifugal) potential of the planet.Due to the fundamental properties of self-gravitation, perturbation in these potentials induces sea-level change, solid-earth deformation, and polar motion.If magnitudes (or trends) of mass redistribution are known (e.g., from satellite observations), such important geodetic signatures can be computed using a simple model of relative sea-level variation.Following the seminal work of Farrell and Clark (1976), the so-called self-gravitating postglacial sea-level model for a viscoelastic, rotating earth has been discussed in several papers (e.g., Wu and Peltier, 1983;Mitrovica and Milne, 2003;Kendall et al., 2005;Hagedoorn et al., 2007;Spada and Stocchi, 2007;Martinec and Hagedoorn, 2014).Here, we briefly summarize the important (and relevant) components of this model.
For a viscoelastic earth, relative sea level at a given space on the earth's surface and time may be defined as the difference between the absolute sea level (i.e., sea surface without any dynamic effect of tides and ocean currents) and the solid-earth surface, assuming that these are measured relative to a common datum (Tamisiea, 2011).Small deviations in these variables from the respective initial states, following the mass conserving redistribution of the earth's surface or interior materials, may be written as follows: where S and N are respective changes in relative and absolute sea level, U is the associated radial displacement of the solid-earth surface, (θ, λ) are spatial coordinates (on the surface of a spherical earth) that represent colatitude and longitude, and t is time.In a physical sense, Eq. ( 1) implies that S is the exact variation of sea surface that would be observed on a measuring stick attached to the solid-earth surface (Farrell and Clark, 1976;Spada and Stocchi, 2007).
In what follows, we assume that the redistribution of surface mass is induced by transport of material into and out of the cryosphere and that there is an associated viscoelastic gravitational response of the solid earth.For the situation where it is mass transport between continental ice and oceans, it is most convenient to define a loading function, L, so that where H is the change in ice thickness, O is the so-called ocean function, ρ I is ice density, and ρ O is ocean water density.By definition, O = 1 for oceans and O = 0 otherwise (Munk and MacDonald, 1960).This function needs to be introduced in Eq. ( 2) because S is defined over the whole planet, including its continents.(Negative S in the continents, i.e., U − N, may be interpreted as the change in surface elevation.)Note that O is not an explicit function of time in our development, but would have to be included in cases of significant submergence, or emergence of coastal lands during mass transport (e.g., Johnston, 1993;Milne, 1998).
The mathematical description of the gravity and loading associated with mass transport requires perturbations in gravitational potential, , and rotational potential, , to enter the absolute sea level as follows: where g is the acceleration due to gravity.Spatial constants appearing above, E and C, are given by where R is the mean radius of the earth, A O is the surface area of oceans, S is the surface domain of a unit sphere, and C is a function of potentials and the associated deformation of the solid-earth surface and is given by C = + −gU .Note that eustatic terms E and C are essential to satisfy the mass conservation constraint (Farrell and Clark, 1976).In a hypothetical, non-gravitating (i.e., = 0), non-rotating (i.e., = 0), rigid (i.e., U = 0) earth, E solely describes S, and it is this metric that is often termed "sea-level equivalent" in order to (alternatively) quantify mass change in glaciers and ice sheets.Sometimes, E by itself is simply termed "eustatic sea level".
Similarly, the viscoelastic gravitational response of solid earth following redistribution of surface mass (Eq.2) may be partitioned for convenience as follows: where U and U are radial displacements of the solid-earth surface associated with perturbations in gravitational and rotational potentials, respectively.
In the following, we briefly present the fundamental concepts and mathematical descriptions of gravitational and rotational potentials, as well as the associated deformation of the solid-earth surface, required to fully define S (Eq.1).Contemporary models are mostly based on the same theory.

Gravitational potential and solid-earth deformation
The general model description presented above may be applied to any earth model, ranging from a simple rigid earth (e.g., Woodward, 1888) to a comprehensive 3-D viscoelastic earth with lateral heterogeneity and nonlinear rheology (e.g., Wu and van der Wal, 2003).Here, we consider the earth as a radially stratified elastic sphere, whose short-term responses are characterized by the so-called load Love numbers (Love, 1911;Longman, 1962) that are referred to the Legendre transform spectral representation of the spherical coordinates on the surface of a sphere.
In order to define and U , we employ Green's function approach to solving for interior earth responses at the surface, essentially following the load Love number formalism for a seismologically constrained elastic earth (e.g., Longman, 1962;Takeuchi et al., 1962).Let G and G U be the non-dimensional Green's functions for a radially stratified, spherically symmetric elastic earth that are, respectively, associated with and U .These functions may be represented in the domain of the Legendre transform as follows: where P l are Legendre polynomials of degree l (see Appendix A), k l and h l are the load Love numbers (Longman, 1962), and α is the arc length between the loading point and the evaluation point on the earth's surface.The load Love numbers appearing above have a simple physical interpretation: P l is the perturbation in degree l representation of non-dimensional gravitational potential in a Legendre transform space induced by the applied mass itself, whereas k l P l and h l P l are similar perturbations for non-dimensional gravitational potential and non-dimensional radial displacement of the solid-earth surface, respectively, caused by the elastic deformation of matter within the earth's interior.Intuitively, G = ∞ l=0 P l and G U = 0 for a rigid earth model.The terms 3gG / 4π R 2 ρ E and 3G U /(4π R 2 ρ E ) express the influence of point load of unit mass on the gravitational potential and radial displacement of the solid-earth surface, respectively (Farrell and Clark, 1976), where ρ E is the average density of the earth.Spatial convolution of these terms with the loading function (Eq.2) gives and U , and we may write where (θ , λ ) are the variable coordinates.These variable coordinates at which the loading function is defined are related to the fixed ones, (θ, λ), at which and U are evaluated via α according to the following cosine formula: cos α = cos θ cos θ + sin θ sin θ cos(λ − λ).

Rotational potential and solid-earth deformation
The surface mass redistribution and associated deformation of solid earth also induce changes in the earth's rotational vector (e.g., Munk and MacDonald, 1960;Lambeck, 1980;Sabadini et al., 1982).The corresponding change in rotational potential deforms both the solid-earth surface and the geoid, thus contributing to a relative sea-level signal.Although geological-timescale perturbations to the rotational vector such as true polar wander are governed by glacial isostatic adjustment (GIA) and mantle dynamics (e.g., Spada et al., 1992;Tsai and Stevenson, 2007;Martinec and Hagedoorn, 2014), short-timescale perturbations of decadal to centennial scale to polar motion are largely determined by cryosphere and other climate driven mass change (e.g., Gross, 2000;Chen et al., 2013;Mitrovica et al., 2015).In the context of mass exchange between continental ice and oceans (Eq.2), it is therefore important to account for rotational feedbacks in the sea-level computation.
In analogy with the description of and U based on the load Love number theory presented in Sect.2.1, we may express and U as follows (e.g., Lambeck, 1980;Milne, 1998): where Y 2mn are degree 2 spherical harmonics (SHs; see Appendix A), 2mn are the corresponding SH coefficients, and k 2 and h 2 are degree 2 tidal Love numbers (e.g., Peltier, 1974;Lambeck, 1980).These Love numbers parameterize the elastic response of the solid earth to a potential forcing that does not involve a direct loading on the earth's surface and have the following physical interpretation: k 2 Y 2mn is the perturbation in degree 2 order m representation of nondimensional rotational potential in a SH transform domain caused by the elastic deformation of matter within the earth's interior, and h 2 Y 2mn is the same for non-dimensional radial displacement of the solid-earth surface.For a rigid earth model, = 2 m=0 2 n=1 Y 2mn and U = 0.In order to define the perturbation 2mn , we consider body fixed right-handed Cartesian coordinates, x i , with the origin located at the center of mass (CM) of the initially equilibrium earth.(x 1 is aligned along the central meridian and x 3 is positive toward the North Pole.)In such a coordinate frame, the products of unperturbed inertia tensor vanish, i.e., I ij = 0 (for i = j = 1, 2, 3), and the moments of unperturbed inertia tensor for a (assumed) rotationally symmetric earth are given by I ii = A (for i = 1, 2) and I 33 = C, where A is the mean equatorial and C is the polar moment of inertia.Similarly, the components of (initially equilibrium, and unperturbed) an angular velocity vectors are given by ω i = δ i3 (for i = 1, 2, 3), where δ i3 are the Kronecker deltas and is the mean rotational velocity of the earth.Following the redistribution of mass (Eq.2), both I and ω are perturbed from their initial equilibrium states.Let I ij and m i be respective perturbation terms, where m i are non-dimensional and typically of order ≤ 10 −6 .Noting the normalization scheme (see Appendix A), 2mn may be defined as follows (Munk and MacDonald, 1960;Lambeck, 1980): When the rotational perturbations are small, m i (t) can be determined from the linearized Liouville equations.These then form the general equation of motion for an elastic rotating earth: where ] is the secular (fluid) Love number, and G is the universal gravitational constant.The Love number k s is a measure of the rotational deformation of a density stratified inviscid earth (Munk and MacDonald, 1960).The variables ψ i (for i = 1, 2, 3) appearing above are the so-called excitation functions and are given by Note that the first terms on the right-hand side are directly induced by mass redistribution, and hence are often called the "mass excitation functions" (Lambeck, 1980).
From the rotational theory presented above, it is clear that 2mn and hence and U can be evaluated completely if three perturbation parameters, namely I i3 , are known.In the present context of mass exchange between continental ice and oceans (Eq.2), we may write (Lambeck, 1980) where L 2mn are degree 2 SH coefficients of the loading function.

Methods
There are certain elements in the relative sea-level theory presented in Sect. 2 that would naturally favor the spectral methods for their numerical evaluation; expansion of nondimensional Green's functions in the form of an infinite sum of Legendre polynomials (Eq.6) is one such example.Indeed, contemporary sea-level models are mostly based on the so-called pseudo-spectral method (Mitrovica and Peltier, 1991) in which all variables appearing in the sea-level equation (SLE) are expanded in the form of SHs and individual SH coefficients are evaluated by satisfying the SLE itself (e.g., Milne, 1998;Mitrovica and Milne, 2003;Spada and Stocchi, 2007).An alternative to the employed Green's function approach to evaluating the (visco)elastic gravitational response of the solid earth (Sect.2.1) is to consider a comprehensive, 3-D finite-element (FE) modeling of the earth (e.g., Gasperini and Sabadini, 1990;Martinec, 2000;Wu, 2004).While the solid-earth response may then account for the lateral heterogeneity as imaged in seismology (e.g., Wu and van der Wal, 2003), it is mathematically and numerically far more cumbersome, a feature that does not diminish in the limit of an elastically responding earth model that we are interested in.We therefore find it convenient to validate or compare our results against the more common pseudo-spectral methods, for which contemporary software is also available (Spada and Stocchi, 2007).Despite the widespread application, one obvious disadvantage of pseudo-spectral methods is that these require large numbers of terms in the series expansion in order to accurately parameterize a slowly converging function such as G U (see Sect. 3.1).The associated basis functions (i.e., highdegree SHs) have short wavelength signals, which demand uniformly distributed high-resolution pixels over the whole planet.Need for high-degree series truncation in conjunction with high-resolution pixels naturally requires a high computational cost.The same statement applies for capturing high-resolution features such as rapid ice melting from an outlet glacier and adjacent sea-level changes: solutions must be evaluated at a large numbers of pixels (as pseudospectral methods require equal pixel size over the whole planet) and high-resolution signals can only be resolved with high-degree SHs (i.e., high-degree series truncation is essential).
Here we present a simple mesh-based computation of SLE on a self-gravitating, elastically compressible, rotating earth that exploits Green's representation of perturbation in gravitational potential and solid-earth deformation evaluated at the surface of the earth (Eq.6).We map the Green's functions on the mesh system with great accuracy, retaining kilometerscale resolution properties.The mesh architecture affords robust control on the number of elements and is thus computationally efficient.Figure 1 shows an example of the computational FE mesh on the solid-earth surface.This mesh is generated using Gmsh (Geuzaine and Remacle, 2009, http://geuz.org/gmsh/),along with anisotropic mesh refinement based on the Bidimensional Anisotropic Mesh Generator (BAMG) package developed by Hecht (2006).The mesh consists of 16 553 vertices and 33 102 elements.Element sizes are restricted to be in the range of about [60,1000] km.The mesh refinement metric used in this particular example is a function of the distance from the nearest coastline.Note that the method presented here is not a formal FE computation of elasto-gravitational deformation of the solid earth, but is an architectural design that allows for solving of the SLE on an unstructured mesh with arbitrary elemental shape or size.The relevant model and material parameters are listed in Table 1.

Evaluation of and U
Crucial to evaluating and U is to accurately sample Green's functions (Eq.6) that are given in the form of an infinite sum of Legendre polynomials.Since −k l decays exponentially but h l approaches very slowly (see Fig. 2a) toward a constant value as l → ∞, the following discussion focuses on accurate parameterization of G U .This discussion, however, equally applies to G as well.
In contemporary models, G U is evaluated by simply truncating the series at degree L, such that Typically, 60 < L < 600.For L = 128, for example, the approximation of G U is characterized by a systematic noise (blue in Fig. 2b) about the exact solution (to be defined later) with higher amplitudes near the loading point, and we may anticipate numerical difficulty in computing changes in bedrock slope or relative sea level near the position of a rapidly changing outlet glacier.It is important to note here that we consider the CM of the earth system reference frame (Blewitt, 2003) in our computations, so that degree 1 Love numbers are on the order h 1 = −1.29 and k 1 = −1.00.The same frame is used, for example, for computing gravity fields from the space geodetic satellites.A much better approximation than Eq. ( 14) would be the following: we may write (e.g., Schrama, 2005) where h ∞ is a constant that is reached asymptotically as l → ∞.We now assume that h l = h L for all l ∈ [L, ∞).
Since G U → −∞ as α → 0, Eq. ( 16) cannot be evaluated at the point of loading, i.e., at α = 0.In order to avoid this inherent singularity, we define the loading function (Eq.2) at the element centroids and evaluate Green's functions at the vertices so that α > 0 for nonzero element size.Let E and V be the total number of elements and vertices in the mesh (Fig. 1).For each vertex v ∈ [1, V], we compute G and G U due to unit loads that are centered at the individual elements e ∈ [1, E] as follows: This mesh is generated using Gmsh (Geuzaine and Remacle, 2009), with the BAMG mesh refinement algorithm (Hecht, 2006).The mesh refinement metric employed here is a function of the distance from the nearest coastline.
where variables with superscripts ve are matrices of size V × E, and those with e are vectors of size E × 1. Figure 2c illustrates how accurately our model samples the exact solution of G U (i.e., for L = 10 000) for an example vertex due to the nearby elemental unit loads.
Once G ve and G ve U are computed and L e (t) are given (to be discussed in Sect.3.3), we may perform the convolution integral (Eq.7) simply as follows: where A e are elemental areas.Note that variables with superscripts v are vectors of size V × 1.

Evaluation of and U
It is essential to compute L 2mn (Eq.13) and m i (Eqs.10-11) prior to evaluating and U .Degree 2 SH coefficients of a loading function can be approximated, according to Eq. (A6), as follows: where Y e 2mn are degree 2 SHs evaluated at elemental centroids (θ e , λ e ) of the mesh.Here, it is important to note that SH-based computations on an unstructured mesh are valid for low SH degrees, such as l = 2, which have very large signal wavelengths.
The system of non-homogeneous ordinary differential equations appearing in Eq. ( 10) can be solved for two unknowns, m 1 and m 2 .The solutions are given by , where c 1 and c 2 are constants to be determined from initial conditions, and We assume that time-dependent variables may be expressed as the sum of their incremental step changes.For instance, , where [δI i3 ] k is the incremental step change in moment of inertia over time t k ≤ t < t k+1 induced by the corresponding incremental step change in applied ice loads and associated sea-level variations (Eq.13), and H(t −t k ) is a Heaviside step function with magnitude of unity for t ≥ t k and zero otherwise.
If incremental step changes in parameters are known a priori (or computed) up to and including the Kth time, it is convenient to reset the time so that τ = t − t K .We may then write, for example, (For each time interval, we are essentially treating variables as if these can be expressed using a single Heaviside step function.)Substituting t by τ in Eqs. ( 20) and ( 21), a simplification of the latter equation follows Similarly, the following can be derived from Eq. ( 11) for m 3 : where c 3 is yet another constant to be determined from initial conditions.
If m i at τ = 0 − (i.e., at time t = t K , but just before imposing the Kth incremental change) are known, from Eqs. ( 20) and ( 23) we may set c i = m i (0 − ).Then m i can be evaluated for any time τ ≥ 0. Setting τ = t K+1 − t K , we can compute m i (0 − ) and hence c i for the subsequent, i.e., (K+1)th, incremental change.For the first incremental change, we impose m i (0 − ) = 0 as initial conditions, assuming the initial equilibrium state of unperturbed ω.
Once m i are computed at a given time t, 2mn can be easily obtained from Eq. ( 9), and the evaluation of and U becomes fairly straightforward as follows: where Y v 2mn are degree 2 SHs evaluated at vertices (θ v , λ v ) of the mesh.Recall that these two quantities are key to computing the rotational feedback.In this formulation we retain the time dependency as the mass changes directly enter the moment of inertia tensor, which, in turn, drive the solutions of m i (Eq.20).

Evaluation of other variables
As noted earlier, we define L at the elemental centroids, allowing evaluation of the entire set of variables, including S, at the vertices.Since L depends on S itself (Eq.2), it is necessary to map S v onto the elemental centroids of the mesh, and we do this by simply averaging the corresponding S v for individual elements.We may now write where v = 1, 2, 3 are the vertices of the eth (triangular) element.In order to evaluate the above equation, H e (t) and O e must be provided.We define O e using the Generic Mapping Tools (Wessel et al., 2013) for a set of coordinates (θ e , λ e ) that define the elemental centroids of the mesh.
Similarly, the eustatic terms appearing in Eq. ( 4) can be evaluated as follows: where The problem is avoided by using the approximation given by Eq. ( 16).(c) Demonstration of model ability to accurately parameterize G U .Solutions obtained by truncating the series (according to Eq. 16) at L = 128 are virtually the same as exact solutions (at least beyond ≈ 1 • from the load), and these are accurately sampled by our model at an example vertex due to unit loads applied at the elemental centroids of the mesh (circles in the figure).(See Sect.3.1 for the explanation of the exact solution.)For comparison, solutions associated with Eq. ( 14) illustrate the ability of contemporary pseudo-spectral models to parameterize G U .
Numerical discretization of all components of SLE is now complete, and these can be easily assembled to evaluate radial displacement of the solid-earth surface (Eq.5) as , and finally change in relative sea level (Eq. 1) as S v (t) = N v (t) − U v (t).

Solution algorithm and model performance
The computational algorithm used in our mesh-based model is similar to that of pseudo-spectral models (Mitrovica and Peltier, 1991) and is as follows: since S v (t) must be known while computing L e (t) and S v (t) itself, we employ a recursive scheme with an initial solution of S v (t) obtained by setting L e (t) = ρ I H e (t) (see Eq. 25).It is then feasible to iteratively compute L e (t) according to Eq. ( 25) with S v (t) obtained from the solution of a previous iteration.We iterate the simulation until a metric that quantifies the change in two successive solutions is minimal.It takes only a few iterations (typically, five to seven) to converge the solution so that the difference in successive solutions is to within the acceptable accuracy (typically, 5 orders of magnitude smaller than the solution itself).This is the standard algorithm for solving the SLE (e.g., Farrell and Clark, 1976;Mitrovica and Peltier, 1991;Spada and Stocchi, 2007), and does not require further explanation.
Once gravitationally consistent solutions for change in relative sea level, S, are obtained, several useful geodetic observables may be retrieved easily.Of particular interest, we may compute radial displacement of the solid-earth surface, U , from Eq. ( 5) and change in absolute sea level, N, from Eq. ( 1).Similarly, we can evaluate the following parameters related to the polar motion of the earth: mass excitation functions, χ i (for i = 1, 2), may be computed using the relationship χ i = I i3 / [C − A]; positions of the North Pole, (p 1 , p 2 ), in the right-handed Cartesian coordinates may be approximated as (p 1 , p 2 ) ≈ (m 1 , m 2 ); and change in the length of a day, D, is given by D = −m 3 D, where D ≈ 86 400 s is the length of a solar day.From some of these solutions, we may also infer other useful geodetic observables, such as changes in absolute gravity and geocentric motion of the earth.These will be further discussed in Sects.4.2 and 4.3.
Most of our computations are done at the vertices of the mesh.Therefore, we have to mainly deal with vectors of size V ×1.Evaluation of Green's functions, however, requires that matrices of size V ×E be considered (Eq.17), and it naturally demands more computer resources.Fortunately, we can compute Green's functions only once at the beginning as a part of the model initialization because these do not evolve as we simulate the model for an assumed elastic earth.For the mesh considered in this study (Fig. 1), which has V = 16 533 and E = 33 102, our Matlab ® code takes about 5 min for a serial run in a MacBook Pro (OS X 10.9.5) to compute Green's functions, and less than a minute to evaluate changes in sea level and associated geodetic parameters caused by instantaneous melting of a synthetic ice sheet.The employed unstructured mesh has elements of varying size in the range of about [60, 1000] km.We may vary this range, for instance, to capture high-resolution features in particular locations, say around the Amundsen Sea Sector (ASS), yet the computational cost is minimal as noted above as long as the new unstructured mesh consists of similar V and E. The lower limit of element size that our model can handle essentially depends on the degree at which series expansion of Green's functions (Eq.17) is truncated.We use L = 10 000 for all of our computations.Assuming P ≈ π R/L (Orszag, 1974), where P denotes the characteristic element size, it implies that our model can capture features of size as small as ≈ 2 km.
To our knowledge, there are no standard benchmark or model intercomparison experiments available in order to test and validate new relative sea-level models that operate on an elastic rotating earth.However, for a suitable set of experiments, we validate key components of our model by reproducing relevant published results as summarized in Appendix B. We now provide a brief comparison of our mesh model with contemporary pseudo-spectral models in terms of computational efficiency.In the latter models, as noted earlier, the SLE is discretized in the SH domain and individual SH coefficients are evaluated.For a chosen spatial resolution, P, the SH series expansion may be truncated at L ≈ π R/P, yielding about (π R/P + 1) 2 SH coefficients to be computed.In order to sample the SH signals uniformly over the whole planet, the pseudo-spectral approach also requires that pixels of equal size be considered.This implies that we must deal with matrices of size ≈ 4π R 2 /P 2 × (π R/P + 1) 2 while evaluating the SLE using pseudo-spectral models.If we are to compute the solutions in spectral formulation, for example, at 60 km resolution along the coastlines as considered in our mesh-based computation, we must deal with matrices of size about 141 685 × 111 947 (compared to vectors of size 16 533×1 required for our model) that certainly demand huge computer resources.In order to compare the model performance systematically, we consider a pseudo-spectral model that is essentially a Matlab ® version of SELEN (Spada and Stocchi, 2007) except the following: (1) the earth's surface is pixelized using MEALPix Toolbox that is a Matlab ® version of HEALPix (http://healpix.jpl.nasa.gov/);(2) the solid earth is treated as an elastic sphere; and (3) the SLE includes rotational feedback.This Matlab ® version of the model, coded by the authors (unpublished), is tested and validated against the original SELEN model for suitable experiments.In Table 2, we compare mesh-based and pseudospectral models in terms of numerical architecture and computational cost.The latter model already demands a large computer resource to capture even a moderate 51 km resolution.It becomes more cumbersome if we seek to deal with higher wave number features like smaller (kilometerscale) ice caps and ice fields using pseudo-spectral models, yet there is little degradation in computational efficiency using our mesh-based approach because we can refine the mesh (down to ≈ 2 km) wherever needed while maintaining the similar numbers of vertices in the mesh.

Some geodetic signatures of ice sheets
Of several climate driven phenomena of mass redistribution on the earth's surface, those related to the cryosphere may be of particular interest.Space-based observations have shown that ice sheets and glaciers expel a large volume of meltwater in an ongoing climate warming (e.g., Shepherd et al., 2012;Gardner et al., 2013), thus directly contributing to the sealevel rise (e.g., Hay et al., 2015).Such trends are likely to persist, if not amplify, throughout this century and beyond as increased atmospheric and oceanic temperatures are generally predicted for future climate change scenarios (e.g., Bamber and Aspinall, 2013;Jevrejeva et al., 2014).Here, as an example model application, we produce computations of some important geodetic signatures associated with the recent evolution of contemporary ice sheets.Observed from space or in terrestrial arrays, these signatures provide diagnostic information about strong shifts in climate (e.g., Chen et al., 2013).

The GRACE data
The twin Gravity Recovery and Climate Experiment (GRACE) satellites are now a way of monitoring and assessing earth's time-varying gravity field caused by the climate driven surface mass redistribution and transportation of materials within the earth's interior (e.g., Wouters et al., 2014).The GRACE data that are now available at an unprecedented resolution of a few hundreds of kilometers have revolutionized our approach to evaluating, for example, glaciers and ice-sheet mass balance (e.g., King et al., 2012;Ivins et al., 2013;Velicogna and Wahr, 2013;Sasgen et al., 2013) and terrestrial hydrological budget (e.g., Wahr et al., 1998;Swenson et al., 2003).
The GRACE data are distributed in the form of Stokes coefficients (Bettadpur, 2012) and, upon standard processing of these SH coefficients, with removal of the mantle GIA signal, we can express the relevant geophysical signals in terms of water height equivalent (WHE).In this analysis, we use GRACE Release-05 Level-2 GSM data products provided by the Center of Space Research, University of Texas at Austin (http://www.csr.utexas.edu/grace/RL05.html).The monthly time series of these data are available up to SH degree and order 60, and cover a period from April 2002 to March 2015 (hereafter referred to as the "GRACE period").There are only partial or no data available for a few months.We fill these data gaps through a simple linear scaling or interpolation between adjacent monthly data as appropriate.We replace degree 1 and degree 2 Stokes coefficients by the values obtained from the analysis of satellite laser ranging (SLR) observations of five passive geodetic satellites (Cheng et al., 2011(Cheng et al., , 2013a, b), b).This is particularly important as our computations predict a degree 1 field related to earth's geocentric motion and a degree 2 field related to polar motion.
We compute Stokes coefficient anomalies for further processing by subtracting the corresponding mean values (over the GRACE period) from individual Stokes coefficients.There may be several techniques of varying complexities to process these data (e.g., Swenson and Wahr, 2002;Chen et al., 2006;Sasgen et al., 2006;Chambers and Bonin, 2012), but all of these involve filtering the unphysical north-south striping patterns that are inherently due to the orbital geometry of the satellites (e.g., Ray and Luthcke, 2006), and reducing the so-called leakage effects that mainly operate between the adjacent sources of signal (e.g., Chen et al., 2006).Here, we generally follow the recipe of Ivins et al. (2013) for recovering WHE from these Stokes coefficient anomalies, except d Our model has variable element size.For L = 10 000, in principle, it can capture features of size as small as ≈ 2 km at minimal computational cost as listed above (provided that the new unstructured mesh consists of similar V and E).
that we employ a GIA computed by A et al. ( 2013) and select a rescaling such that the linear trends in ice mass loss from the Antarctic Ice Sheet (AIS) and Greenland Ice Sheet (GrIS) during the GRACE period are −90 and −240 Gt yr −1 , respectively.Throughout, we apply a Gaussian smoothing with a 300 km radius.This smoothing radius may be large enough to filter the short-wavelength noise, yet small enough to retain the actual geophysical signals, and is in the range of typical values recommended for variety of applications (Wahr et al., 1998;Swenson et al., 2003).
The temporal and spatial trends in the final products of the AIS and GrIS mass balance are shown in Fig. 3.The amplitudes of temporal variability are higher for the AIS, implying the large seasonal mass turnover there, but it could also be due to large signal amplitudes of the GIA model used (A et al., 2013).A kink is apparent in the AIS data (Fig. 3a); the ice sheet has lost more mass (−111 Gt yr −1 ) since 2007 than during the early period (−40 Gt yr −1 ).On the other hand, the GrIS was losing the mass more steadily with a little seasonal turnover, but at a much greater pace (−240 Gt yr −1 ).Spatial distributions of the (linear) mass balance trend are shown in Fig. 3b.The figure suggests that the AIS was mostly losing the mass from the ASS at a large rate (dH /dt ≤ −10 cm yr −1 ), and also from the Peninsula and Wilkes Land at a more modest rate.Note that the AIS has also gained mass in a few locations; there are clear signals of mass gain (at a rate of about 2 cm yr −1 ) in the northern East Antarctic Ice Sheet (EAIS).The GrIS lost mass at a much greater pace of about −20 cm yr −1 , mainly from the southeastern sector during the first half of the GRACE period, but losses also expanded toward the west later on.Figure 3b also suggests small mass changes in the north and interior of the GrIS during the GRACE period.All these features are generally consistent with other published solutions, for both AIS and GrIS (e.g., Velicogna and Wahr, 2013).However, it is important to note that our primary goal here is to demonstrate the predictive capabilities of our mesh-based model rather than computing precise mass budget solutions for the polar ice sheets.In Sect. 5 we describe some application of the methods to high-resolution ice-sheet models.

Sea level and other variables
The monthly time series of H (θ, λ, t) for both ice sheets are obtained from the GRACE data as discussed above.We force loading of the model by these mass balance solutions for the two major ice sheets.Our model computes monthly solutions for relative sea level, S, radial displacement of solidearth surface, U , and perturbation in absolute sea level, N .Figure 4 summarizes these solutions for a combined forcing of ice sheets, where we show the linear trends in variables obtained by fitting the corresponding monthly solutions in a least-square sense.Figure 4a depicts the trend in S with the following key features: a large rate of sealevel drop (dS/dt ≤ −3.0 mm yr −1 ) with large wavelengths around the GrIS; the same, but with relatively smaller wavelengths around the ASS; and a moderate rate of sea-level rise (1.5 mm yr −1 ) in the northern EAIS.The blue contours represent the global mean rate (GMR) of relative sea-level rise with a magnitude of about 0.91 mm yr −1 .The corresponding values for the AIS and GrIS are about 0.25 and 0.66 mm yr −1 , respectively.In the regions enclosed by these contours, sea level either falls or remains unchanged or rises more slowly than the GMR.In the exterior regions, on the other hand, sea level rises at a higher pace than the GMR.
From an ice-sheet modeling point of view, U and N are perhaps more important variables than S itself because these provide direct constraints to two of the important boundary conditions, namely the bedrock elevation and the sea surface height.Figure 4b shows the spatial distribution of the linear trend in U , with the same general features as observed for S (Fig. 4a) but of opposite sign.The solid-earth  , where X = N or S (for the Northern and Southern Hemisphere) and 0 = 1, . .., 7 (see Fig. 4a for their position on the global map).
surface uplift is predicted to occur at a relatively large rate (dU/dt ≥ 2.5 mm yr −1 ) around the ASS and GrIS, and subsides at a rather moderate rate of −0.5 mm yr −1 around the northern EAIS.The trend in N shows much greater variability in space (Fig. 4c).Here the computation is performed in the CM reference frame and the corresponding Green's function does not change monotonically (unlike those for S and U ) as a function of great-circle distance from the loading point.Note in Fig. 2b how G U , for example, increases monotonically as the evaluation point moves away from the load.The predicted change in absolute sea level drops at a rate of dN/dt ≤ −1.0 mm yr −1 around the ASS and GrIS, while a rise is predicted at a similar rate in the northern EAIS and the northern Pacific.Generally speaking, the relative sea level, sea surface height and solid-earth deformation are linearly related to each other (Eq.1).Therefore, sea-level drop is generally accompanied by the earth surface uplift and sea surface fall, and sea-level rise is by the earth surface subsidence and sea surface rise.
It may be useful to evaluate the corresponding changes in absolute gravity (gravity anomaly or disturbance), because this geodetic variable may be measured directly using absolute gravimeters (e.g., James and Ivins, 1998;Crossley et al., 2012) and space geodetic satellites.It may be possible to compute this variable on the same computational mesh as we use for solving the SLE, but it is more readily estimated from the solutions of N in the SH domain.Here we evaluate gravity anomaly, g, that could be measured by a gravimeter on the earth's surface as follows (Lambeck, 1980): where N lmn are the SH coefficients of N. The linear trend in g is shown in Fig. 4d.As expected, large negative trends are visible around the GrIS and ASS (d g/dt ≤ −1.0 µgal yr −1 ), where mass was being lost rapidly during the GRACE period (Fig. 3b).Similarly, in the regions with mass accumulation such as the northern EAIS (Fig. 3b), rising trends in gravity anomaly are predicted, but at a rather moderate rate of 0.5 µgal yr −1 .(A Gaussian filter of 100 km smoothing radius is employed for these plots.)  3 for their description).
An interesting exercise is to compare the relative contribution of individual ice sheets to the total sea-level change.
(Total sea-level change, in the context here, is due to the combined mass evolution of both AIS and GrIS.)For such analysis, we select 14 representative tide gauge stations, half of which are located in the Northern Hemisphere.The de-scriptions of these sites are given in Table 3 and their coordinates on the global map are shown in Fig. 4a.For two representative sites (one for each hemisphere), Fig. 5a shows the explicit evolution of sea-level change and the relative contribution of AIS and GrIS.In Honolulu, the total sea level rises faster (black line in the figure) than the GMR (red line) throughout the GRACE period.The GrIS contribution at this site is higher (light blue fill) than its contribution to the global mean value (blue line).The AIS influence at this site is similar.This is summarized in the figure inset, in which we compare the average trends in sea-level variation: the local contributions of GrIS (dS/dt = 0.82 mm yr −1 ) and AIS (0.33 mm yr −1 ) are both greater than the corresponding GMRs (i.e., 0.66 and 0.25 mm yr −1 , respectively), thus resulting in a much greater pace of local sea-level rise in Honolulu (1.15 mm yr −1 ) than the GMR (0.91 mm yr −1 ).By contrast, the total sea level falls at a rapid pace in the Pine Island Glacier (−2.88 mm yr −1 ) despite the positive and larger than global mean contribution of GrIS (0.82 mm yr −1 ), and it is mainly due to the local gravitational loss (Fig. 4d) associated with a strong loss in ice mass from the ASS (Fig. 3b).
Figure 5b summarizes a similar comparison for dS/dt at other tide gauge stations (see Table 3 and Fig. 4a).Although the figure is self-explanatory, we make brief remarks on some interesting features.On both the western coast (San Francisco; N4) and the eastern coast (Virginia Keys; N5) of the continental United States, contributions of AIS are greater and those of GrIS are smaller than the corresponding GMRs.Their combined effects, however, are somewhat contrasting: the pace of local sea-level rise (0.93 mm yr −1 ) is slightly higher than the GMR at San Francisco, and the opposite is true for Virginia Keys, where sea level rises at a rather modest rate of 0.82 mm yr −1 .Greatly contrasting rates are predicted at two closely located places, namely Reykjavik (N6) and Newlyn (N7).Since the AIS contributions are similar at both sites, differing signatures in total sea-level change are due to the GrIS: the ice sheet has a strong negative contribution at Reykjavik, causing the total sea level to drop at the great rate of −1.07 mm yr −1 .Its contribution at Newlyn, on the other hand, is minimal, and therefore the total sea level there rises much more slowly than the GMR (0.34 mm yr −1 ).Note that there are other interesting comparisons such as between the eastern (Casey; S5) and western (Rothera; S7) limits of the AIS.At Casey in the EAIS, the AIS contribution is minimal and the total sea-level rise (0.67 mm yr −1 ) is mainly due to the GrIS.However, at Rothera, both ice sheets have similar contributions but opposite signs, thus resulting in a virtually stagnant sea level (0.05 mm yr −1 ).In the end, it is also worthwhile reporting the overwhelmingly great pace of sea-level rise (1.54 mm yr −1 ) at Syowa in the northern EAIS (S4).The AIS contribution is large there, mainly due to the enhanced gravitational pulling (Fig. 4d) associated with the local mass gain (Fig. 3b).
What we have highlighted thus far are the global predictive features that can be efficiently extracted from our flex- ible FE mesh system.In the section that follows we apply those model predictive features to two important geodynamical observables associated with space geodesy: polar motion and geocentric motion of the earth.

Polar motion and geocentric motion
The redistribution of mass on the earth's surface (Eq.2) excites changes in the position of the rotational spin axis with respect to fixed positions in the crust.This is also known as changes in "earth orientation", or is commonly known as polar motion (Lambeck, 1980), and has important observational ties to the mass imbalance of the earth's two great ice sheets (Douglas et al., 1990).The GMR of sea-level rise caused by ice sheets is roughly 1 mm yr −1 since about 2005 (Shepherd et al., 2012).If rates increase during the next 80-100 years, sufficient to reach the upper end of the estimate of Jevrejeva et al. (2014) for 2100 (≈ 180 cm), then it will be important to have the state-of-the-art data assimilating ice-sheet mod-els, such as ISSM, incorporate such a profoundly fundamental observable.Here we compute some important attributes of polar motion induced by contemporary ice sheets and associated sea-level changes during the GRACE period.
Figure 6a and b shows 3-D plots for the monthly position of the North Pole.Complex interactions between the near-annual forcing and Chandler (433-day period) wobble results in a net polar wobble with varying amplitude.This can be seen in the figures for both ice sheets.While wobbling around the mean rotational axis, the pole also drifts away from its initial position, as indicated by trend lines in the figures.A kink in the drift direction is apparent for the AIS in about 2007, which may be linked to a similar feature observed in mass evolution of the ice sheet (see Fig. 3a).
In order to predict polar drift from our mesh-based computational framework, we evaluate classic mass excitation functions, χ 1 and χ 2 , associated with individual ice sheets (Fig. 6c) and the corresponding annual pole positions (Fig. 6d).The North Pole in Fig. 6d represents that of year 2002.The mass loss in the GrIS yields positive χ 1 and negative χ 2 in the employed right-handed Cartesian system, implying that the general drift direction is toward the fourth quadrant defined by x 1 > 0 and x 2 < 0 (see Fig. 6d).Since χ 1 and χ 2 are on a similar order of magnitudes, the GrIS induced drift is directed toward the ice sheet itself along ≈ 40 • W longitude.On the other hand, the combination of rapid mass loss in the ASS and mass gain in the northern EAIS yields positive excitation functions, thus causing the North Pole wander vector to be directed toward a Eurasian ≈ 60 • E longitude.Since the AIS and GrIS operate in-phase for χ 1 and outof-phase for χ 2 , their combined effects would amplify the former excitation function and shrink the latter one.Consequently, the pole would be drifting along ≈ 15 • W longitude during the GRACE period as shown in Fig. 6d.
These predictions are generally consistent with the report by Chen et al. (2013): glaciers and ice sheets explain a large fraction of observed χ 1 during 2005-2011, but their contributions to χ 2 are minimal.We find dχ 1 /dt = 2.87 ± 0.15 mas yr −1 for the same period, which is approximately half the trend attributed by Chen et al. (2013) to glaciers and ice sheets.Large rates of mass loss from other glaciated regions such as the high-altitude Himalayas (Gardner et al., 2013) may explain the discrepancy, although more rigorous effort is needed to justify this.Despite the minimal collective contributions of glaciers and ice sheets to the observed trend in χ 2 , it is important to highlight the significant role that the AIS is playing to counter the GrIS induced negative χ 2 (see Fig. 6c).Rapid mass loss from the ASS, aided by the mass gain in the northern EAIS, is responsible for drifting the pole along ≈ 15 • W longitude, which would otherwise be heading the GrIS.Observations of further eastward motion of the pole (along ≈ 15 • E longitude; Fig. 6d) may be explained by mass transport and other excitations unrelated to ice sheets.
From gravitationally consistent surface mass redistribution, the mesh model may also estimate the geocentric mo- tion of the earth.While observationally more elusive than polar motion, this is a fundamental parameter important to global reference frames.The geocentric motion is caused by the shift in relative position between the CM of the earth system and the center of figure (CF) of the solid-earth surface, and this information is essential to reconcile the geodetic data that are tracked from the ground stations using absolute gravimeters and also from the passive geodetic satellites using SLR.Let the CM-CF shift be denoted by the position vector X i (for i = 1, 2, 3) in the right-handed Cartesian system.The components of this vector can be computed from the degree 1 SH coefficients of a gravitationally consistent loading function (Eq.2) as follows (e.g., Wu et al., 2006): where δ 1i , δ 2i , and δ 3i are the Kronecker deltas, and l 1 is the degree 1 load Love number, just like h 1 and k 1 , and in the employed CM reference frame is on the order l 1 = −0.89(Blewitt, 2003).
The ice-sheet induced components of X during the GRACE period are plotted in Fig. 7. Since both ice sheets are located near the poles (i.e., with small |x i | for i = 1, 2 and large |x 3 |), their individual contributions to the CM-CF shift are naturally larger for X 3 , which is associated with a degree 1 zonal harmonic.However, the ice sheets have secular trends in this component of geocentric motion that oppose one another.This opposition mutes their combined signal (Fig. 7c), predicting a gradual shift toward the South Pole at a rate of −0.44 ± 0.03 mm yr −1 .Both seasonal amplitudes and secular trends in horizontal components of the geocentric motion, i.e., X 1 and X 2 , are quite minimal, compared with the corresponding solutions inferred from the SLR observations (Cheng et al., 2013b), despite the in-phase functioning of ice sheets for the latter component.Incorporation of degree 1 predictions into ice-sheet models such as ISSM will be important for considering geodetic reference frame stability, but it is unlikely to be relegated to the status of a data assimilation parameter due to problems with inferring decadal timescale CM-CF drift (Ries, 2013), especially when compared to other global space geodetic observables.

Sea-level fingerprints of high-resolution ice-sheet forcing
In Sect. 4 we have demonstrated the wide range of applications of ISSM-SESAW.Our loading functions (Eq.2) in these experiments are based on the GRACE inferred mass budget of polar ice sheets, which are inherently of low reso- lution (on the order of 250-350 km), allowing us to consider a relatively coarse mesh (with smallest elements on the order of 60 km).However, one of the major goals of ISSM-SESAW is to be able to resolve the kilometer-scale source of earth's time-varying surface mass transport (e.g., fast-flowing outlet glaciers, and rapid changes in grounding line position) and project it globally.In order to demonstrate this, we provide an example computation of sea-level fingerprint associated with high-resolution modeling of polar ice sheets.We model the dynamics of polar ice sheets over highresolution mass conserving beds (Morlighem et al., 2014, and personal communication).We mesh the ice sheets anisotropically from 50 km (AIS) and 15 km (GrIS) at the ice divide to 3 and 1 km in fast-flowing regions, respectively, with respective ice shelves meshed at 9 and 1 km.These meshes comprise a total of 102 308 (AIS) and 93 039 elements (GrIS).For both ice sheets, we perform a 200year control run using a 2-D shelfy-stream approximation (MacAyeal, 1993) for ice-flow mechanics.The models are initialized from steady state using basal friction inversion (Larour et al., 2012;Schlegel et al., 2015) and present-day InSAR surface velocities (Rignot et al., 2011).The models are initially relaxed to thermal steady state as well (Seroussi et al., 2013).The model forcings include a surface climatology provided by Regional Climate Model RACMO2.3 (Ettema et al., 2009) averaged during 1979-2010, and melt rates under Antarctic ice shelves obtained from mean annual simulation of ECCO (Estimating the Circulation and Climate of the Ocean) over 1992-2012(Heimbach, 2008)).The global sea-level model comprises the actual high-resolution meshes of both ice sheets, and is anisotropically adapted from 1 km in the highest-velocity areas of GrIS to 2000 km in the middle of the Pacific Ocean, for a total of 292 262 elements.Coastlines are resolved with the Generic Mapping Tool (Wessel et al., 2013), and the mesh itself is generated using Gmsh (Geuzaine and Remacle, 2009).We compute relative sea level on this global mesh by forcing ISSM-SESAW with high-resolution ice thickness change over a 200-year control run on both polar ice sheets.Solutions are shown in Fig. 8. Large negative changes in sea level are evident around both ice sheets, symptomatic of the overall loss of mass from both AIS and GrIS during the control runs.
These solutions are obtained by running ISSM-SESAW on NASA's Pleaides supercomputer.Due to the intensive nature of the convolution operation (i.e., a load defined at every elemental centroid contributes to sea level evaluated at every vertex of the mesh), a low-latency bandwidth network is required to distribute solutions across the network in a dense pattern.Pleaides has an InfiniBand ® interconnect, with all nodes connected in a partial hypercube topology, which makes this operation possible.Scaling is highly dependent upon the convolution loop, which needs to be carried out in a fragmented way, to space out the distribution of values across the cluster.This operation is CPU dependent, but is also highly limited by the interconnect speed.Without the discretization of SLE that operates on a flexible unstructured mesh system (presented in this paper) and use of a massively parallelized and fully scalable capability of ISSM-SESAW, such a high-resolution computation of relative sea level and associated global geodetic observables would have been virtually impossible.
While this experiment demonstrates the overall capability of our model development, we acknowledge that the earth has non-negligible creep deformation over the timescales of 200 years (e.g., Dietrich et al., 2010;Klemann et al., 2008), which the current version of ISSM-SESAW cannot account for.Note, however, that the dominant vertical displacement and potential field changes at wavelengths comparable to the ocean basins respond with much longer timescales than 200 years due to the stiffness of the mantle below ≈ 420 km (e.g., Lambeck and Johnston, 1998).The current version of the model also does not provide systematic solid-earth/sealevel feedbacks to the ice-flow mechanics, although we have fully embedded high-resolution ice-sheet models within the global mesh.We are currently working toward this two-way model coupling, and this may be presented in a future publication.

Conclusions
The motivation for this study is in concert with the rapid developments in ice-sheet modeling that have occurred over the past 5-10 years, wherein processes that occur at a kilometer scale must be captured along with local sea-level variability.Toward developing a coherent set of ice-sheet and solidearth/sea-level models that operates on a common computational architecture provided by JPL's ISSM (https://issm.jpl.nasa.gov/),we present a mesh-based approach to evaluating gravitationally consistent relative sea-level and associated geodetic variables.Unlike contemporary sea-level models that are mostly based on SH formulation, the model can function efficiently on an anisotropic unstructured mesh, thus capturing the physics operating at kilometer scale yet capable of simulating geophysical quantities that are inherently of global scale with minimal computational cost.
In order to explain the global model, we compute the evolution of sea-level fingerprints and other observables, such as sea surface height, gravity anomaly and solid-earth deformation, associated with GRACE inferred monthly mass balance of ice sheets for a period from April 2002 to March 2015 in a manner that is broadly familiar to the space geodesy and altimetry communities (e.g., Farrell and Clark, 1976;Riva et al., 2010;Mitrovica et al., 2011).We also evaluate the corresponding polar and geocentric motion of the earth and find that both ice sheets play a significant role in explaining the observed eastward drift of the North Pole since about 2005 (Chen et al., 2013), whereas the predicted influences on earth's geocentric motion are minimal compared with the SLR inferred estimates (Cheng et al., 2013b).One of the greatest strengths of our model is to resolve kilometer-scale sources of earth's surface mass transport and predict the corresponding local and global geodetic signatures.We demonstrate this by providing an example computation of sea-level fingerprint associated with high-resolution modeling of polar ice sheets.
Global geodetic and sea-level signatures that can be computed using the mesh model may have important implications for earth system modeling.Coupling a global sea-level model to a local mesh of a 3-D ice-sheet model, for example, enhances the realistic simulation of outlet glaciers, such as Pine Island Glacier, as it provides a direct constraint to two of the important boundary conditions, namely the bedrock elevation and the sea surface height, that would be consistent with global-scale climate driven mass redistribution.There may yet be several other applications that involve continentalscale gravitational and loading interaction.However, the current model development is strictly applied to an elastically compressible and density layered self-gravitating earth and, hence, suitable for short-timescale (monthly to decadal) evaluation of variables.For relatively long-timescale (centennial or longer) computations, the model should also account for viscoelastic response of the solid earth.It may be achieved through appropriate parameterization of long-term GIA response via time-dependent viscoelastic Love numbers.In the second experiment, we compute change in relative sea level caused by instantaneous melting of a synthetic ice sheet and compare it with the corresponding solution obtained from a Matlab ® version of SELEN (Spada and Stocchi, 2007).We assume that a disc-shaped ice sheet, which has an 18 • radius and is centered at 60 • north latitude and 85 • west longitude, melts uniformly.The pseudo-spectral model consists of 12 288 pixels and the SH expansion is truncated at L = 128.In order to match these settings as closely as possible, our mesh model also consists of structured (triangular) elements with 12 514 vertices (and 25 024 elements) and considers the traditional approximation of Green's functions (as in Eq. 14 rather than Eq.16) with L = 128.Solutions, shown in Fig. B2, show remarkable resemblance considering the differences in model setup.Notice, for example, the location of Z = 100 % contours.There are, however, minor (somewhat random) differences between the solutions (Fig. B2c), particularly around the deglaciated region, which may be considered as artifacts due to differing numerical approaches.In order to highlight the importance of high wave number signatures, we also obtain a solution by having L = 10 000 in our mesh model.Since all other computational settings are exactly the same, the difference in solutions shown in Fig. B2d illustrates the systematic errors associated with severity of the low-degree SH truncation at 128.The third experiment assesses the ability of our mesh model to provide rotational feedback to relative sea level, surface displacement and gravitational potential.We consider the exact same model setup as for the second experiment, but include the perturbation in rotational potential caused by polar motion and associated solid-earth deformation in SLE.Solutions are plotted in Fig. B3.Since the rotation induced sea level varies spatially over time in our formulation (see Sect. 3.2), this solution is retrieved by filtering out a dominant signal of 433-day Chandler wobbles.As the ice sheet loses mass, the earth's spin axis drifts toward a position aligned with the axis of the new principal moments of inertia, thus causing the sea level to drop locally with a characteristic SH degree 2 order 1 signature.While this is not a benchmark, per se, it is conceptually and quantitatively consistent with the rotational feedback reported in the numerical experiments (e.g., Mitrovica et al., 2009, Fig. 1c).

Figure 1 .
Figure 1.Example of unstructured mesh at earth surface.Both the (a) Northern and (b) Southern Hemisphere are shown, with continents depicted in cyan.This mesh is generated using Gmsh(Geuzaine and Remacle, 2009), with the BAMG mesh refinement algorithm(Hecht, 2006).The mesh refinement metric employed here is a function of the distance from the nearest coastline.

Figure 2 .
Figure 2. Parameterization of (elastic) solid-earth deformation caused by surface loading.(a) Load Love numbers, h l , up to 1800 Legendre degree.As l → ∞, h l converges slowly toward a constant value.(b) Non-dimensional Green's function, G U , computed by truncating the series at L = 128.The conventional approximation (Eq.14) produces noise, with greater amplitudes near the loading point.The problem is avoided by using the approximation given by Eq. (16).(c) Demonstration of model ability to accurately parameterize G U .Solutions obtained by truncating the series (according to Eq. 16) at L = 128 are virtually the same as exact solutions (at least beyond ≈ 1 • from the load), and these are accurately sampled by our model at an example vertex due to unit loads applied at the elemental centroids of the mesh (circles in the figure).(See Sect.3.1 for the explanation of the exact solution.)For comparison, solutions associated with Eq. (14) illustrate the ability of contemporary pseudo-spectral models to parameterize G U .

Figure 3 .
Figure 3. Summary of the GRACE data used in this study.(a) Ice mass change in the AIS and GrIS from April 2002 to March 2015.(b) Spatial distribution of rate of change in ice thickness, averaged over the GRACE period.See Sect.4.1 for a detailed discussion of the GRACE data processing.

Figure 4 .
Figure 4. Some important geodetic signatures of ice sheets during the GRACE period.Rates of change in (a) relative sea level, (b) solid-earth deformation, (c) absolute sea level, and (d) absolute gravity.(Notice the different scale and color order in the color bars.)These results are obtained by linearly fitting the corresponding monthly solutions in a least-square sense.The blue contours in (a) represent the trend in the global mean value, with magnitude dS/dt = 0.91 mm yr −1 .Annotations are supplied for 14 locations in (a) (see Table3for their description).

Figure 5 .
Figure 5. Magnitudes and trends of sea-level change at 14 selected locations.(See Table 3 for their description and Fig. 4a to locate their position on the global map.)(a) Change in relative sea level over time at two representative sites.Average rates of relative sealevel change are also shown in the insets.(b) Rate of change in relative sea level, averaged over the GRACE period, at 12 other locations.Cyan and blue colors (both fill and line) represent the contribution from the AIS and GrIS, respectively.The combined contributions of ice sheets are shown in red.Black lines in Fig. 5a represent the local values and all others denote the global mean values.

Figure 6 .
Figure 6.Polar motion during the GRACE period.Monthly pole positions (with respect to April 2002), (p 1 , p 2 ), excited by (a) AIS and (b) GrIS mass loss in the right-handed coordinate system (x 1 , x 2 ) (see d). Red lines show the average drift directions.A kink is apparent for the AIS.(c) Mass excitation functions, χ 1 and χ 2 , associated with individual ice sheets.(d) Annual pole positions, after removing Chandler wobbles, with respect to the 2002 position.For comparison, the observed long-term (green arrow; Mitrovica et al., 2006) and recent (2005-2011; blue arrow; Chen et al., 2013) drift directions are also shown.Note that 1 mas ≈ 3.09 cm.

Figure 7 .
Figure 7. Geocentric motion during the GRACE period.The CM-CF shift, with respect to the April 2002 position, along the (a) x 1 , (b) x 2 , and (c) x 3 directions in a right-handed Cartesian system.(SeeFig.6dfor a positive sense of horizontal axes; the vertical axis, x 3 , is positive out of the North Pole.)Note that different scales are used (in the right y axis) for the SLR-based estimates(Cheng et al., 2013b).

Figure 8 .
Figure 8. Sea-level fingerprints of high-resolution ice-sheet forcing.The upper frame shows relative sea level computed by ISSM-SEASAW at the end of a 200-year control run of both polar ice sheets.The corresponding ice forcings (e.g., Schlegel et al., 2015) are shown in the lower frame.Close-ups on Greenland (frames a and c) and Antarctica (b, d) are provided for convenience.We do not account for rotational feedback in this particular set of computations.

Figure B2 .
Figure B2.Change in relative sea level induced by deglaciation on an elastic non-rotating earth.Normalized sea level with respect to global mean value (see Eq. B1), computed by (a) our mesh model and (b) a pseudo-spectral model for similar computational settings.Z = 100 % contours are plotted (in blue) to facilitate comparison.(c) Difference between the two model solutions.Some (random) discrepancies appearing particularly around the deglaciated region may be considered as numerical artifacts due to differing model algorithms.(d) Difference between two mesh models with L = 128 and L = 10 000 (Eq. 14).Large discrepancies appear around the deglaciated region, implying the systematic errors associated with low-degree SH representation of solid-earth response.

Figure B3 .
Figure B3.Rotational feedback to relative sea level.(a) Normalized sea level computed by ISSM-SESAW as in Fig. B2a, but with rotational feedback included.Notice the shape of Z = 100 % (blue contours) with respect to the one drawn in Fig. B2a.(b) Difference between (a) and Fig. B2a, illustrating the (SH degree 2 order 1) signature of rotation induced sea-level variation.

Table 2 .
Comparison of pseudo-spectral and mesh-based computations of the SLE.We denote element (or pixel) size by P, number of elements (or pixels) by E, number of vertices by V, and the degree at which (associated) Legendre series is truncated by L. (V is not relevant for the pseudo-spectral approach.)It is the model initialization time.For pseudo-spectral models, most of this time is required to compute SHs that would be of size E × (L + 1) 2 .Our model mainly utilizes this time to compute Green's functions, which are of size V × E. b It is the CPU time required to solve the SLE following instantaneous melting of a hypothetical ice sheet.Pseudo-spectral models once again deal with matrices of size E × (L + 1) 2 .Our model mostly deals with vectors of size V × 1. c Both pseudo-spectral and mesh-based models are coded in Matlab ® and simulated in a MacBook Pro (OS X 10.9.5).We employ the Parallel Computing Toolbox ™ of Matlab ® with four local workers in parallel runs.Serial runs use a single worker. a

Table 3 .
Location of selected tide gauge stations, with Global Sea Level Observing System (GLOSS) ID.
• W * ID has the following format: X0