the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Sealevel fingerprints emergent from GRACE mission data
Surendra Adhikari
Erik R. Ivins
Thomas Frederikse
Felix W. Landerer
Lambert Caron
The Gravity Recovery and Climate Experiment (GRACE) mission data have an important, if not revolutionary, impact on how scientists quantify the water transport on the Earth's surface. The transport phenomena include land hydrology, physical oceanography, atmospheric moisture flux, and global cryospheric mass balance. The mass transport observed by the satellite system also includes solid Earth motions caused by, for example, great subduction zone earthquakes and glacial isostatic adjustment (GIA) processes. When coupled with altimetry, these space gravimetry data provide a powerful framework for studying climaterelated changes on decadal timescales, such as ice mass loss and sealevel rise. As the changes in the latter are significant over the past two decades, there is a concomitant selfattraction and loading phenomenon generating ancillary changes in gravity, sea surface, and solid Earth deformation. These generate a finite signal in GRACE and ocean altimetry, and it may often be desirable to isolate and remove them for the purpose of understanding, for example, ocean circulation changes and postseismic viscoelastic mantle flow, or GIA, occurring beneath the seafloor. Here we perform a systematic calculation of sealevel fingerprints of onland water mass changes using monthly Release06 GRACE Level2 Stokes coefficients for the span April 2002 to August 2016, which result in a set of solutions for the timevarying geoid, seasurface height, and vertical bedrock motion. We provide both spherical harmonic coefficients and spatial maps of these global field variables and uncertainties therein (https://doi.org/10.7910/DVN/8UC8IR; Adhikari et al., 2019). Solutions are provided for three official GRACE data processing centers, namely the University of Texas Austin's Center for Space Research (CSR), GeoForschungsZentrum Potsdam (GFZ), and Jet Propulsion Laboratory (JPL), with and without rotational feedback included and in both the centerofmass and centeroffigure reference frames. These data may be applied for either study of the fields themselves or as fundamental filter components for the analysis of oceancirculation and earthquakerelated fields or for improving ocean tide models.
© 2019 California Institute of Technology. Government sponsorship acknowledged.
Geodesists have long understood that the ocean mean sea surface follows the shape of the Earth's geoid (Rapp, 1983) and that changes in onland water storage are a source of timevarying gravity (Lambert and Beaumont, 1977). The fundamental relationship of changes in land ice and water, solid Earth, and seasurface height is essential to the study of past and present relative sea level (e.g., Peltier, 1982; Clark et al., 2002; Tamisiea, 2011). Our recent gain in confidence for monitoring the geographic locations and amplitudes of both seasonal and supraseasonal changes in global glacier and ice sheet mass, dating to the beginning of the radar interferometry and altimetry era of the early 1990s (e.g., Rignot et al., 2011; Shepherd et al., 2018), strengthens our ability to effectively harness this information to construct informative models about global sealevel variability associated with a selfattraction and loading phenomenon (Spada and Galassi, 2016; Larour et al., 2017; Mitrovica et al., 2018). The mathematical formalism relating changes in gravitational, rotational, and solid Earth deformation responses to land ice and hydrological mass change has now niched itself into contemporary studies of sealevel change: the prediction of “sealevel fingerprints”. Sealevel fingerprints are a consequence of the fact that the water elements being transported laterally between land and oceans carry mass, gravitational attraction, and the ability to change the radial stress at the solid Earth surface. These are characterized, for example, as changes in relative sea level encircling areas of intense ice mass loss such as Patagonia, coastal Alaska, the Amundsen Sea sector of West Antarctica, and the Greenland Ice Sheet (e.g., Mitrovica et al., 2001; Tamisiea et al., 2014; Riva et al., 2010; Adhikari and Ivins, 2016).
To date, space gravimetric measurements using Gravity Recovery and Climate Experiment (GRACE) monthly gravity fields and the subpolar ocean altimetry measurements from TOPEX/Poseidon and Jason each have multiple geophysical signals and respective noise floors that are generally high enough that clear detection of these contemporary landmassdriven fingerprints in the oceans has remained elusive. However, it is believed that these signals will eventually emerge in these data systems. Such a belief springs, in part, from the fact that amplitudes of internal ocean variability in intra and interannual mass that GRACE observes are relatively mute in comparison to onland hydrology, twoway landtoocean transport, and secular trends in land ice changes (Chambers, 2006; Chambers and Willis, 2008; Watkins et al., 2015; Wiese et al., 2016; Save et al., 2016). In fact, Hsu and Velicogna (2017) have used oceanbottom pressure data, in conjunction with space geodetic data, to claim that fingerprints associated with decadalscale onland mass changes are detected. Furthermore, Davis and Vinogradova (2017) have shown that the fingerprints of Greenland ice mass loss have had measurable influences on tide gauge records along the eastern coast of the US since the mid1990s. Galassi and Spada (2017) have noted that the influence of a landmassinduced fingerprint may be reflected in tide gauge records of relative sea level at the northern Antarctic Peninsula, as there is a distinct change in trend at about the year 2000 CE, possibly reflecting increased regional ice mass loss. Each of these observations might be considered both intriguing and preliminary in terms of providing the community with unambiguous detection of sealevel fingerprints.
The effects of the fingerprints are nonetheless important to disentangle from many geophysical and ocean circulation models and dataset. New insights into the regional and global sealevel budgets are sought through explicitly combining ocean altimetry with the space gravimetry information, and a key part of this combination is to account for the details of sealevel fingerprints (e.g., Rietbroek et al., 2012; Frederikse et al., 2018). Consideration of landice and waterdriven fingerprints is also necessary when using geodetic data to search for glacial isostatic adjustment (GIA) signals residing at or beneath the seafloor (e.g., Simon et al., 2018) or examining ice loss on land when ocean water surrounds the region, such as in Graham Land of the Antarctic Peninsula (Ivins et al., 2011; Sterenborg et al., 2013). Future applications of sealevel fingerprints in geophysical geodesy should include the study of great earthquakes (M_{w}≥8.0) at subduction zones (Han et al., 2016) and at ocean rifts beneath the open ocean (Han et al., 2015) or adjacent to the Antarctic Ice Sheet (King and SantamaríaGómez, 2016).
This paper describes a dataset of monthly changes in relative sea level, geoid height, and vertical bedrock motion induced by mass redistribution from land to ocean. These are derived from the Release06 GRACE Level2 monthly Stokes coefficients for the period April 2002 to August 2016. The GRACE mission data have been instrumental to the study of the Earth's climate system (Tapley et al., 2019) and have helped us resolve numerous longstanding questions in oceanography, hydrology, cryosphere, and geodesy. The GRACE gravity solutions are now employed for providing new insights into changes in ocean circulation (Johnson and Chambers, 2013; Landerer et al., 2015; Saynisch et al., 2015; Mazloff and Böning, 2016). The terrestrial water storage is now rigorously quantified for continents (Rodell et al., 2015; Hirschi and Seneviratne, 2017) as is the global cryospheric mass balance (Velicogna, 2009; Jacob et al., 2012; Ivins et al., 2013; Luthcke et al., 2013; Schrama et al., 2014). Land mass change and its exchange with the global oceans, in fact, makes it possible to successfully reconstruct subtle changes in the position of Earth's spin axis on interannual timescales (Adhikari and Ivins, 2016), thus providing a confidence in the robustness of GRACEbased estimates of global surface mass transport.
Relative sea level is defined as the height of the ocean water column bounded by two surfaces: solid Earth surface and sea surface. Change in relative sea level, ΔS, at a geographical location described by colatitude and longitude (θ,ϕ) over the time interval Δt may be expressed as follows:
where ΔN and ΔU are corresponding changes in seasurface height and bedrock elevation, respectively. Both of these variables are usually expressed relative to the reference ellipsoid, which in turn is defined relative to either the centerofmass (CM) of the total Earth system or the centeroffigure (CF) of the solid Earth surface. Tide gauges provide direct measurements of ΔS, whereas satellite altimetry measures ΔN in the CM reference frame.
Mass redistributed on Earth's surface provides a direct perturbation to the Earth's gravitational and rotational potentials, causing a corresponding perturbation in the geoid height. Since the geoid height on a realistic Earth does not necessarily have to coincide with the seasurface height, we write
where ΔΦ is the net perturbation in Earth's surface potential, ΔC is a spatial invariant that explains the discrepancy between the seasurface height and geoid height (Tamisiea, 2011), measured with respect to the same reference ellipsoid, and g is the mean gravitational acceleration at Earth's surface. As will be further discussed in Sect. 3, ΔC is essential to conserve mass. Spacebased gravity missions, such as GRACE and GRACE FollowOn (GRACEFO), provide direct measurements of the nonrotational part of ΔΦ (to be defined explicitly in Sect. 3) in the CM frame; the satellite system cannot measure the rotational part of ΔΦ because it retrieves data in an inertial reference frame.
In geodetic applications, global field variables are typically expanded in a spherical harmonic (SH) domain. Most of the GRACE data processing centers – including the University of Texas Austin's Center for Space Research (CSR), GeoForschungsZentrum Potsdam (GFZ), and Jet Propulsion Laboratory (JPL) – provide monthly solutions for normalized SH coefficients of the gravitational potential termed “Stokes coefficients”. Stokes coefficient anomalies – the values that deviate from the mean (static) field – can be used to readily retrieve changes in onland ice and water storage or ocean bottom pressure. The goal of this paper is to provide Stokes coefficient anomalies (i.e., SH coefficients of ΔΦ minus rotational centrifugal potential) associated with the sealevel fingerprint of monthly changes in onland ice and water storage, which are derived from CSR, GFZ, and JPL GRACE Stokes coefficients themselves. As we shall further clarify below, we provide these new fingerprint coefficients and their corresponding spatial maps computed in both CM and CF reference frames, with and without rotational feedback included. We also provide solutions for change in relative sea level ΔS and bedrock elevation ΔU. Corresponding solutions for seasurface height, ΔN, may be retrieved using Eq. (1). For brevity, we hereafter drop the Δ symbol and assume that variables imply “change” in respective fields – not the absolute fields – implicitly.
Here we briefly summarize the fundamental concept and a numerical technique of solving the socalled “sealevel equation”. Much of the background and supporting materials may be found, for example, in Farrell and Clark (1976), Mitrovica and Peltier (1991), Adhikari et al. (2016), and Spada (2017). Let $L(\mathit{\theta},\mathit{\varphi},t)$ be the global, massconserving load function so that
where $H(\mathit{\theta},\mathit{\varphi},t)$ is the onland change in water equivalent height over the time period t, and $S(\mathit{\theta},\mathit{\varphi},t)$ is the corresponding change in relative sea level on the oceanic domain 𝒪(θ,ϕ). By definition, 𝒪=1 for the oceans and 0 otherwise. For ease of discussion, we write $F(\mathit{\theta},\mathit{\varphi},t)=H(\mathit{\theta},\mathit{\varphi},t)\left[\mathrm{1}\mathcal{O}(\mathit{\theta},\mathit{\varphi})\right]$ so that $F(\mathit{\theta},\mathit{\varphi},t)$ defines the model “forcing” function.
The net change in onland (water) mass directly affects the relative sea level, hence conserving mass on a global scale. Such a redistribution of mass on Earth's surface perturbs its gravitational and rotational potentials and further redistributes the ocean mass. The net result of these perturbations is the sealevel fingerprint: a unique spatial pattern of relative sea level that is consistent with fundamental physical features of a realistic Earth. For a selfgravitating elastically compressible rotating Earth, we compute the sealevel fingerprint by satisfying the following sealevel equation:
The physical interpretation of the righthand side terms is as follows:

The barystatic term, E(t), directly follows from the mass conservation principle. This spatial invariant describes S that would be resulted in by distributing the net change in land water storage uniformly over the oceans.

Changes in the Earth's surface potential, $\mathrm{\Phi}(\mathit{\theta},\mathit{\varphi},t)$, and the solid Earth surface, $U(\mathit{\theta},\mathit{\varphi},t)$, may be partitioned as follows:
$$\begin{array}{}\text{(5)}& {\displaystyle}\left\{\begin{array}{c}\mathrm{\Phi}(\mathit{\theta},\mathit{\varphi},t)\\ U(\mathit{\theta},\mathit{\varphi},t)\end{array}\right\}=\left\{\begin{array}{c}{\mathrm{\Phi}}_{\mathrm{g}}(\mathit{\theta},\mathit{\varphi},t)+{\mathrm{\Phi}}_{\mathrm{r}}(\mathit{\theta},\mathit{\varphi},t)\\ {U}_{\mathrm{g}}(\mathit{\theta},\mathit{\varphi},t)+{U}_{\mathrm{r}}(\mathit{\theta},\mathit{\varphi},t)\end{array}\right\},\end{array}$$where Φ_{g} and U_{g} are the respective signals associated with the perturbation in gravitational potential. We may compute Φ_{g} and U_{g} by convolving L (Eq. 3) with the respective Green's functions. Similarly, Φ_{r} and U_{r} are associated with the perturbation in rotational potential. The change in Earth orientation driven by shift in the inertia tensor causes both solid Earth deformation and sealevel change (Lambeck, 1980). The net effects of the change in orientation of Earth's spin axis thus provide a rotational feedback (e.g., Milne and Mitrovica, 1998). We may compute Φ_{r} and U_{r} based on the perturbation in Earth's inertia tensor due to the global surface mass redistribution described by L (Eq. 3). We define all of the terms appearing in Eqs. (4) and (5) explicitly in Appendix A.

The last term in Eq. (4) represents the oceanaveraged value of $(\mathrm{\Phi}/gU)$. This spatial invariant is essential to ensure that the global mean relative sealevel change is the same as the barystatic term.
To solve for the sealevel fingerprint in a conventional SH domain (e.g., Mitrovica and Peltier, 1991) and isolate useful SH coefficients noted in Sect. 2, we express Eq. (4), using Eq. (5), in the following form:
where $X={\mathrm{\Phi}}_{\mathrm{g}}/g$, $Y={\mathrm{\Phi}}_{\mathrm{r}}/g$, $P={U}_{\mathrm{g}}$, $Q={U}_{\mathrm{r}}$, and $C=E<\mathrm{\Phi}/gU>$. By default, we account for the rotational feedback, which when excluded, Eq. (6) takes a reduced form with Y=0, Q=0, and $C=E<{\mathrm{\Phi}}_{\mathrm{g}}/g{U}_{\mathrm{g}}>$. We now multiply both sides of Eq. (6) by the ocean function, 𝒪, to get the following:
where $\widehat{S}=\mathcal{O}S$, $\widehat{X}=\mathcal{O}X$, and so on. In the employed spectral methods (Appendix B), we find it more straightforward to solve Eq. (7) rather than Eq. (6). Since all of the righthand side terms appearing in Eq. (7) depend on $\widehat{S}$ itself (see Eqs. 3 and 4 and Appendix A), we solve the equation recursively until the desired solution convergence is achieved (see Appendix B5). We consider the barystatic sea level (Eq. A1) as the starting solution; i.e. $\widehat{S}=\widehat{E}$, where $\widehat{E}=\mathcal{O}E$. Once Eq. (7) is solved for $\widehat{S}$, all of the terms appearing in Eq. (6) may be retrieved easily.
We expand all of the terms appearing in Eqs. (6) and (7) in the SH domain (as in Eq. B1). Inserting these SH expansions into Eq. (7) and equating the corresponding (degree l, order m) SH coefficients, we find the following for any rth recursion:
where $r=\mathrm{1},\mathrm{2},\mathrm{\dots}{r}_{\mathrm{max}}$ is the recursion counter, and r_{max} is the value of r for which the desired convergence is attained. Note that dependence of righthand side terms on ${\widehat{S}}_{lm}$ itself is explicitly stated. For r=1, we set ${\widehat{S}}_{lm}^{\mathrm{0}}={\widehat{E}}_{lm}$. Since ${\widehat{E}}_{lm}$ does not depend on ${\widehat{S}}_{lm}$, it does not evolve during the recursion. We define ${\widehat{E}}_{lm}$ (Eq. B17) and all of the other hatted coefficients appearing above (Eq. B18 and so on) in Appendix B. The hatted coefficients depend on corresponding nonhatted coefficients, which are explicitly defined in Eqs. (B20)–(B23) and (B25).
Once we obtain the final solution for ${\widehat{S}}_{lm}^{r}$ (after iteration r=r_{max}), denoted for simplicity by ${\widehat{S}}_{lm}$, final solutions for all of the nonhatted (degree p, order q) coefficients are obtained as well. These nonhatted coefficients automatically satisfy the sealevel equation itself (Eq. 6) in the SH domain; i.e.,
Note that all of the SH coefficients appearing above are only a function of time t. With final solutions achieved for all of the terms appearing in Eq. (9), SH coefficients of geoid height change for a selfgravitating Earth are given by X_{pq}(t) and those for a selfgravitating rotating Earth by [X_{pq}(t)+Y_{pq}(t)]. Similarly, SH coefficients of bedrock elevation change are given by −P_{pq}(t) for a selfgravitating Earth and by $\left[{P}_{pq}\left(t\right)+{Q}_{pq}\left(t\right)\right]$ for a selfgravitating rotating Earth.
Here we give a brief summary of the steps undertaken to develop sealevel fingerprints and complementary data products. First, we note that the GRACE processing centers, including CSR, GFZ, and JPL, have a variety of methods employed to reduce noise, but the system has an inherent resolution limit of about 300 to 400 km in radius at the Earth's surface. Hence, the Stokes coefficients for the potential field provided by the official centers are truncated at a varying degree and order, from 60 to 96. We employ a truncation at degree and order 60, as many months may be much noisier than others.
We use GRACE Level2 Release06 data products provided by all three premier (and official) data processing centers (available at ftp://podaac.jpl.nasa.gov/allData/grace/L2/, last access: 7 May 2019) that are available for the spans April 2002 through August 2016 (CSR and JPL) and January 2003 through November 2014 (GFZ). The Release06 GSM files represent the total gravity variability due to land surface hydrology, cryospheric changes, episodic seismogenic processes, and GIA. We assume that all mass transport information is contained within the postprocessed GSM files in which background models for the mass changes in atmosphere and oceans having periodicities shorter than 1 month are removed (Dobslaw et al., 2017). GSM datasets are also corrected for solid Earth and ocean tides by the processing centers (see Stammer et al., 2014; Bettadpur, 2018). We also assume continuous transfer of net mass to and from the oceans takes place on all timescales. This includes a trend that supplies the mass term of sealevel rise. To do this correctly, we derive degree 1 coefficients from JPL Release06 data products following the methods of Swenson et al. (2008). We replace degree 2 order 0 coefficients by those derived from satellite laser ranging analysis (Cheng et al., 2011) that are compatible with Release06 data products (available at ftp://podaac.jpl.nasa.gov/allData/grace/docs/TN11_C20_SLR.txt, last access: 7 May 2019). The physical origins motivating this replacement are well known: there is far greater sensitivity to changes in degree 2 order 0 coefficients that can be retrieved from higher orbiting satellites tracked by terrestrial laser stations than for GRACE (Cheng and Ries, 2017). We apply GIA correction coefficient by coefficient using the expected values from a Bayesian analysis (Caron et al., 2018), available at https://vesl.jpl.nasa.gov/solidearth/gia/ (last access: 7 May 2019). Finally, for all coefficients, we remove corresponding 11year (January 2003–December 2013) mean values to retrieve Strokes coefficient anomalies.
By combining GSM Stokes coefficient anomalies with GIA and lowdegree coefficients as noted above, we may derive corresponding coefficients for land water storage anomalies, H_{lm}(t), as follows (Wahr et al., 1998):
where ρ_{w} is the water density, ρ_{e} is the Earth's mean density, ${k}_{l}^{\prime}$ are the load Love numbers of degree l, a is the Earth's mean surface radius, r is the Gaussian smoothing window, and ${C}_{lm}^{\mathrm{gsm}}$ and ${C}_{lm}^{\mathrm{gia}}$ are the GSM and GIA Stokes coefficient anomalies, respectively. The term enclosed by braces is the Gaussian smoothing filter. We consider r=300 km to comply with the socalled gain factors that are used to restore the attenuated signals (detailed below). An asterisk associated with GSM coefficients is meant to imply that these solutions are corrected for more accurate lowdegree Stokes coefficients as noted above.
Monthly land water storage fields, $H(\mathit{\theta},\mathit{\varphi},t)$, may be generated by assembling the coefficients (Eq. 10) in an SH domain, as in Eq. (B1). Gaussian smoothing aimed at removing the data noise also attenuates the signals. An appropriate scaling of the fields is therefore essential. For the ice sheet and peripheral glaciers in Greenland, three nonoverlapping subdomains of Antarctica, and 15 regions of global glaciers and ice caps, we compare our estimates of average rate of regional mass change during February 2003 through June 2013 with those computed by Schrama et al. (2014) and derive the scaling factors – unique for CSR, GFZ, and JPL data products – for each of these 19 cryospheric domains. As for the noncryospheric continental domains, Landerer and Swenson (2012) analyzed monthly land water storage signals obtained from the GRACE observations and the Noah land surface model, simulated within the Global Land Data Assimilation System (GLDASNoah), and derived global gridded gain factors. We combine these factors to scale $H(\mathit{\theta},\mathit{\varphi},t)$ for the entire continents. Our estimates of barystatic time series are comparable to, for example, JPL mascon solutions for both trends and seasonal amplitudes (Fig. 1).
A detailed description of scaling may be found in Adhikari and Ivins (2016), who used the same recipe to postprocess the CSR Release05 GRACE Level2 Stokes coefficients for robust reconstruction of interannual variability in position of Earth's spin axis. This gives us extraordinary confidence that the procedure for generating land water storage fields and corresponding fingerprints are not only sound but highly robust at long wavelengths. The effects of scaling on SH coefficients of select fields are shown in Fig. 2. Our model solutions are also robust. For example, our estimates of relative sealevel change (Fig. 3), vertical bedrock motion (not shown), and geoid height change (not shown) are consistent with the respective solutions computed using a wellvalidated sealevel solver (Adhikari et al., 2016) that operates on an unstructured global mesh of the Ice and Sealevel System Model (ISSM; https://issm.jpl.nasa.gov/, last access: 7 May 2019). We find that GFZ solutions are slightly different from CSR and JPL solutions (Fig. 4), although the difference in sealevel fingerprints is generally within 1σ uncertainties (compare Figs. 3d and 4d). We show the origin of discrepancies by plotting the degree variance spectrum.
Based on CSR, GFZ, and JPL Stokes coefficients, we provide with this article monthly SH coefficients of

model forcing function, F_{lm}(t);

geoid height change, [X_{lm}(t)+Y_{lm}(t)];

vertical bedrock motion, $\left[{P}_{lm}\left(t\right)+{Q}_{lm}\left(t\right)\right]$; and

relative sealevel change, S_{lm}(t),
computed in both CM and CF reference frames with and without the rotational feedback included. Effects of Earth's rotation and the reference frame origin on select fields are shown in Fig. 5. The SH coefficients for seasurface height may be obtained by summing coefficients for bedrock motion and those for relative sealevel change (see Eq. 1). While one may readily assemble these coefficients in an SH domain to retrieve the corresponding monthly fields, we also supply $\mathrm{0.5}{}^{\circ}\times \mathrm{0.5}{}^{\circ}$ gridded solutions for user convenience.
We provide uncertainty associated with monthly fields as well, both in terms of spatial maps and SH coefficients. Quantification of the uncertainty is determined by the following recipe. Based on the JPL Release06 (GIA uncorrected) mascon solutions and associated standard errors (Watkins et al., 2015; Wiese et al., 2016), we use a Monte Carlo approach to generate 5000 ensemble members of monthly land water storage solutions. We apply a unique GIA correction, computed by Caron et al. (2018), to each of these ensemble members. Next we solve the sealevel equation to derive an equivalent number of solutions for $S(\mathit{\theta},\mathit{\varphi},t)$, $U(\mathit{\theta},\mathit{\varphi},t)$, and other fields. Finally, we quantify the standard errors associated with each field, weighted by the likelihood of each GIA model (Caron et al., 2018). Figure 3 shows our estimates of standard errors associated with the trends in land water storage and relative sea level.
The utility of the data we provide is that they may be used to rigorously remove those patterns that are attributable to geoid height change and bedrock motions caused by onland mass changes from ocean altimetry, bottompressure, and tide gauge studies. Such removal is essential for future studies of the patterns of sealevel change owing to internal variability of the climate system which drives changes in ocean density, fresh water fluxes, and circulation (e.g., Bilbao et al., 2015; Fasullo and Nerem, 2018).
As we supply sealevel fingerprints and complementary data products with and without rotational feedback, we owe the readers some additional words of caution and recommendations. First, from the Eulerian equations of rotational motion, we solve for the feedback consistently designed for periods longer than 434 d (the period of the Chandler wobble). The rationale is that both the solid Earth and ocean pole tide (Haubrich and Munk, 1959) are removed from the satellite solutions for GRACE gravimetry and TOPEX/Poseidon and Jason altimetry on a routine basis (e.g., Wahr, 1985; Desai, 2002; Desai and Yuan, 2006). The improvements in the ocean pole tide, in fact, have been accomplished by many years of assimilation of the altimetric mission data. Hence, at periods near, or less than, 434 d, the paths to unambiguously generating solutions to the sealevel equation with centrifugal potential and loading changes from the pole tide are unclear. We might assume that the relevant feedbacks are largely removed as a processing step in rendering Level 1b and Level 2 GRACE data products. We keep, however, rotational feedback effects of an interannual nature in one set of monthly solutions, and another set of solutions lack these effects. The users of these data should understand the differences, as those employing approaches to using the data to evaluate altimetric time series of order 10 years in length will certainly be interested in using the rotational feedback version for the analysis of interannual trends and variability adjacent to Greenland, for example Müller et al. (2019), whereas users focusing on seasonal timescale fingerprints are recommended to employ those coefficients that lack the rotational feedback, as the altimetry and space gravimetry products employed likely have the seasurface height and gravity effects of the annual polar motion, Chandler wobble, and associated pole tides removed.
It is also worthwhile to note that on timescales of decades the mantle primarily behaves elastically, perhaps with the exception at places where the tectonic history has brought heat, volatiles, and changes in mineral structure, such as water or reduced grain size, into the region, thus reducing the effective viscosity to values below about 5×10^{18} Pa s (e.g., Lange et al., 2014; Mitrovica et al., 2018; Nield et al., 2018). At such low values of viscosity in the upper mantle, stress relaxation can reduce both the effective influence of gravitational loading and the amplitude of fingerprints. While we acknowledge that this effect is quite difficult to quantify, it should be a secondorder effect.
We presently store data in a public repository hosted by Harvard Dataverse (https://doi.org/10.7910/DVN/8UC8IR; Adhikari et al., 2019).
The first set of data we supply are SH coefficients of global field variables. The zip file “SLFsh_coefficients.zip” contains a total of 1780 data files: 133×4 for GFZ and 156×4 each for CSR and JPL. For the given data center, four files are provided for a particular GRACE month: with and without Earth's rotational feedback included while solving for sealevel fingerprints in both CM and CF reference frames. File names follow the GRACE naming convention. Solutions that correspond to the GSM file “GSM2_20020952002120_GRAC_UTCSR_BA01_0600”, for example, are stored in four files named “SLF2_20020952002120_GRAC_UTCSR_BA01_0600” under appropriate directories; we simply replace “GSM” by “SLF” to denote “sealevel fingerprints”. The time stamp (in YYYYDoYYYYYDoY format) and the corresponding data center (fivecharacter string containing CSR, GFZ, or JPL) also appear in the file name. The example file considered above contains sealevel fingerprint solutions for the period 95–120 d of year 2002 based on the Stokes coefficients provided by CSR. Header lines 5–7 in each file further clarify which data center the solutions correspond to, which reference frame is considered, and whether or not Earth's rotational feedback is accounted for. Each data file consists of a total of 18 columns: SH degree l, SH order m, and SH coefficients for model forcing function F_{lm} (four columns), relative sea level S_{lm} (four columns), geoid height change [X_{lm}+Y_{lm}] (four columns), and vertical bedrock motion $[{P}_{lm}+{Q}_{lm}]$ (four columns). For each field, the first (last) two columns store cosine (sine) coefficients for our predicted mean and 1σ uncertainty, respectively. Users should note that the finite degree 0 order 0 harmonic in the monthly SLF files represents the finite mass changes for the global oceans.
The second set of data we supply are gridded maps of global field variables. We provide a total of 12 NetCDF files: four each for CSR, GFZ, and JPL. The file “SLFgrids_GFZOP_CF_WITHrotation.nc”, for example, stores solutions based on GFZ Stokes coefficients that are computed in the CF reference frame with the rotational feedback accounted for.
In this paper we describe a data product that emerges from the Release06 GRACE Level2 Stokes coefficients, provided by CSR, GFZ, and JPL, which contain the basic information necessary to create monthly sealevel fingerprints, and these are general enough that they may be employed in reconstructions of vertical bedrock motion, perturbed relative sea surface, and geoid height change. We provide SH coefficients of each field and uncertainty therein, computed in both CM and CF reference frames with and without rotational feedback included. For user convenience, we also provide spatial maps at $\mathrm{0.5}{}^{\circ}\times \mathrm{0.5}{}^{\circ}$ resolution.
A future space altimetry mission (Surface Water and Ocean Topography, or SWOT) is aimed at providing realtime twodimensional imaging of the seasurface height without the necessity of having to patch together onedimensional profiles (e.g., Gaultier et al., 2016). In addition to providing higher resolution, this will allow improved accuracy. When coupled to GRACEFO mapping of gravity changes, we should begin to see the emergence of sealevel fingerprints. Perhaps more importantly, we may begin to more confidently remove a part of the ocean altimetry signal that should not be assimilated into dynamic ocean models: that which is associated with selfgravitation and loading. Here we present the effects of landbased mass transport and rotational effects, both together and separately. Recent appreciation of the effects of solid Earth elastic and viscoelastic response is now receiving increased scrutiny for the potential bias that may be introduced into the altimetry trend record when not accounting for these effects properly (e.g., Desai et al., 2015; Lickley et al., 2018). We have not treated the influences of GIA on rotational deformation and/or the associated axial displacement of the centrifugal potential, although we have employed the GIA model of Caron et al. (2018) to analyze GRACE Level2 for proper representation of the monthly water height equivalent masses. As a consequence, users of the data that we supply here should understand that folding the fingerprints into analyses of any geodetic data, including tide gauges and ocean altimetry, might want to carefully consider that the secular polar motion effects in the Release06 GRACE Level2 products (Wahr et al., 2015; Bettadpur, 2018) have been removed and that the seasurface height variability associated with polar drift, annual, and Chandler wobble effect is currently removed from ocean altimetry data in the manner described by Desai et al. (2015). This fact allows users to rather straightforwardly remove landmasschangerelated fingerprints from either GRACE, ocean altimetry, tide gauge, or GPSdetermined vertical land motion data from April 2002 to August 2016 using the monthly solutions we supply here as a data product.
The fundamental theoretical concept of the socalled sealevel fingerprint is summarized in Sect. 2. Here we provide explicit mathematical expressions for all of the terms appearing in Eq. (4).

The barystatic term is given by
$$\begin{array}{}\text{(A1)}& {\displaystyle}E\left(t\right)={\displaystyle \frac{\mathrm{1}}{{A}_{\mathrm{O}}}}\left[{a}^{\mathrm{2}}\int H(\mathit{\theta},\mathit{\varphi},t)\left\{\mathrm{1}\mathcal{O}(\mathit{\theta},\mathit{\varphi})\right\}\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}\right],\end{array}$$where
$$\begin{array}{}\text{(A2)}& {\displaystyle}{A}_{\mathrm{O}}={a}^{\mathrm{2}}\int \mathcal{O}(\mathit{\theta},\mathit{\varphi})\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}\end{array}$$is the ocean surface area, a is already defined in Eq. (10), and 𝕊 is the surface domain of a unit sphere. The term enclosed by brackets in Eq. (A1) yields the net change in continental water volume.

Changes in gravitational potential, Φ_{g}, and associated changes in Earth's surface displacement, U_{g}, are obtained by convolving the surface loading function (Eq. 3) with respective Green's functions, 𝒢_{Φ} and 𝒢_{U}, as follows:
$$\begin{array}{}\text{(A3)}& {\displaystyle}\left\{\begin{array}{l}{\mathrm{\Phi}}_{\mathrm{g}}(\mathit{\theta},\mathit{\varphi},t)\\ {U}_{\mathrm{g}}(\mathit{\theta},\mathit{\varphi},t)\end{array}\right\}={a}^{\mathrm{2}}{\mathit{\rho}}_{\mathrm{w}}\int \left\{\begin{array}{l}{\mathcal{G}}_{\mathrm{\Phi}}\left(\mathit{\alpha}\right)\\ {\mathcal{G}}_{U}\left(\mathit{\alpha}\right)\end{array}\right\}L({\mathit{\theta}}^{\prime},{\mathit{\varphi}}^{\prime},t)\phantom{\rule{0.125em}{0ex}}\mathrm{d}{\mathbb{S}}^{\prime},\end{array}$$where $({\mathit{\theta}}^{\prime},{\mathit{\varphi}}^{\prime})$ are the variable coordinates. These variable coordinates at which the loading function is defined are related to (θ,ϕ) at which Φ_{g} and U_{g} are evaluated via the greatcircle distance, α, as follows: $\mathrm{cos}\mathit{\alpha}=\mathrm{cos}\mathit{\theta}\mathrm{cos}{\mathit{\theta}}^{\prime}+\mathrm{sin}\mathit{\theta}\mathrm{sin}{\mathit{\theta}}^{\prime}\mathrm{cos}({\mathit{\varphi}}^{\prime}\mathit{\varphi})$. Green's functions are represented in the Legendre transform domain as follows:
$$\begin{array}{}\text{(A4)}& {\displaystyle}\left\{\begin{array}{l}{\mathcal{G}}_{\mathrm{\Phi}}\left(\mathit{\alpha}\right)\\ {\mathcal{G}}_{U}\left(\mathit{\alpha}\right)\end{array}\right\}={\displaystyle \frac{\mathrm{3}}{\mathrm{4}\mathit{\pi}{a}^{\mathrm{2}}{\mathit{\rho}}_{e}}}\sum _{l=\mathrm{0}}^{\mathrm{\infty}}\left\{\begin{array}{c}g\left(\mathrm{1}+{k}_{l}^{\prime}\right)\\ {h}_{l}^{\prime}\end{array}\right\}{\mathcal{P}}_{l}\left(\mathrm{cos}\mathit{\alpha}\right),\end{array}$$where 𝒫_{l} are Legendre polynomials (Eq. B3), and ${k}_{l}^{\prime}$ and ${h}_{l}^{\prime}$ are the load Love numbers.

Changes in rotational potential, Φ_{r}, and associated changes in Earth's surface displacement, U_{r}, follow from the Eulerian theory of rotation (Lambeck, 1980):
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}\left\{\begin{array}{c}{\mathrm{\Phi}}_{\mathrm{r}}(\mathit{\theta},\mathit{\varphi},t)\\ {U}_{\mathrm{r}}(\mathit{\theta},\mathit{\varphi},t)\end{array}\right\}=\left\{\begin{array}{c}\mathrm{1}\\ \mathrm{0}\end{array}\right\}{\mathrm{\Lambda}}_{\mathrm{00}}\left(t\right){\mathcal{Y}}_{\mathrm{00}}(\mathit{\theta},\mathit{\varphi})\\ \text{(A5)}& {\displaystyle}& {\displaystyle}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}+\sum _{m=\mathrm{2}}^{\mathrm{2}}\left\{\begin{array}{c}\left(\mathrm{1}+{k}_{\mathrm{2}}\right)\\ {h}_{\mathrm{2}}/g\end{array}\right\}{\mathrm{\Lambda}}_{\mathrm{2}m}\left(t\right){\mathcal{Y}}_{\mathrm{2}m}(\mathit{\theta},\mathit{\varphi}),\end{array}$$where 𝒴_{lm} are degree l order m spherical harmonics (Eq. B2), Λ_{lm} are SH coefficients of perturbation in rotational potential, and k_{2} and h_{2} are degree 2 tidal Love numbers. We may express changes in rotational potential in terms of changes in Earth's rotation parameters, moment of inertia, and hence surface loading function. Considering leadingorder terms only, we get the following nonzero coefficients (Milne and Mitrovica, 1998):
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}\left\{\begin{array}{l}{\mathrm{\Lambda}}_{\mathrm{21}}\left(t\right)\\ {\mathrm{\Lambda}}_{\mathrm{2}\mathrm{1}}\left(t\right)\end{array}\right\}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{15}}}}{a}^{\mathrm{2}}{\mathrm{\Omega}}^{\mathrm{2}}\\ \text{(A6)}& {\displaystyle}& {\displaystyle}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\left({\displaystyle \frac{\mathrm{\Omega}(\mathrm{1}+{k}_{\mathrm{2}}^{\prime})}{\mathbb{A}{\mathit{\sigma}}_{\mathrm{0}}}}\left[{\displaystyle \frac{\mathrm{4}\mathit{\pi}}{\sqrt{\mathrm{15}}}}{\mathit{\rho}}_{\mathrm{w}}{a}^{\mathrm{4}}\left\{\begin{array}{c}{L}_{\mathrm{21}}\left(t\right)\\ {L}_{\mathrm{2}\mathrm{1}}\left(t\right)\end{array}\right\}\right]\right),\end{array}$$and
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}\left\{\begin{array}{l}{\mathrm{\Lambda}}_{\mathrm{00}}\left(t\right)\\ {\mathrm{\Lambda}}_{\mathrm{20}}\left(t\right)\end{array}\right\}=\left\{\begin{array}{c}\mathrm{2}/\mathrm{3}\\ \mathrm{2}/\left(\mathrm{3}\sqrt{\mathrm{5}}\right)\end{array}\right\}{a}^{\mathrm{2}}{\mathrm{\Omega}}^{\mathrm{2}}\\ \text{(A7)}& {\displaystyle}& {\displaystyle}\phantom{\rule{0.25em}{0ex}}\left({\displaystyle \frac{\mathrm{1}+{k}_{\mathrm{2}}^{\prime}}{\mathbb{C}}}\left[{\displaystyle \frac{\mathrm{8}\mathit{\pi}}{\mathrm{3}}}{\mathit{\rho}}_{\mathrm{w}}{a}^{\mathrm{4}}\left\{{L}_{\mathrm{00}}\left(t\right){\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{5}}}}{L}_{\mathrm{20}}\left(t\right)\right\}\right]\right),\end{array}$$where Ω is the Earth's mean rotational velocity, 𝔸 and ℂ are the mean equatorial and polar moment of inertia, respectively, σ_{0} is the socalled Chandler wobble frequency, and L_{lm} are SH coefficients of the surface loading function (Eq. B11). Note that the terms inside brackets represent changes in Earth's moment of inertia: ΔI_{11} and ΔI_{22} (Eq. A6) and ΔI_{33} (Eq. A7). Similarly, the terms enclosed by outer parentheses represent Earth's rotation parameters: polar motion (m_{1},m_{2}) (Eq. A6) and change in length of day m_{3} (Eq. A7).

The oceanaveraged term in Eq. (4), denoted by 〈*〉, may be written as follows:
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}\u2329{\displaystyle \frac{\mathrm{1}}{g}}\mathrm{\Phi}(\mathit{\theta},\mathit{\varphi},t)U(\mathit{\theta},\mathit{\varphi},t)\u232a={\displaystyle \frac{\mathrm{1}}{{A}_{\mathrm{O}}}}[{a}^{\mathrm{2}}\int \mathit{\{}{\displaystyle \frac{\mathrm{1}}{g}}\mathrm{\Phi}(\mathit{\theta},\mathit{\varphi},t)\\ \text{(A8)}& {\displaystyle}& {\displaystyle}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}U(\mathit{\theta},\mathit{\varphi},t)\mathit{\left\}}\mathcal{O}\right(\mathit{\theta},\mathit{\varphi}\left)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}\right].\end{array}$$When rotational feedback is excluded, Φ and U should be replaced by Φ_{g} and U_{g}, respectively.
B1 Primer

Spherical harmonics. For brevity, we define $\mathit{\omega}=(\mathit{\theta},\mathit{\varphi})$ and drop explicit dependence of a function on time so that $f(\mathit{\theta},\mathit{\varphi},t)\equiv f\left(\mathit{\omega}\right)$. Any squareintegrable function f(ω) can be expanded as the infinite sum of SHs as follows:
$$\begin{array}{}\text{(B1)}& {\displaystyle}f\left(\mathit{\omega}\right)=\sum _{l=\mathrm{0}}^{\mathrm{\infty}}\sum _{m=l}^{l}{f}_{lm}{\mathcal{Y}}_{lm}\left(\mathit{\omega}\right)\equiv \sum _{lm}{f}_{lm}{\mathcal{Y}}_{lm}\left(\mathit{\omega}\right),\end{array}$$where f_{lm} are SH coefficients, and 𝒴_{lm}(ω) are (real) normalized SHs of degree l and order m. These SHs may be expressed in terms of associated Legendre polynomials, ${\mathcal{P}}_{l\leftm\right}$, as follows:
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}{\mathcal{Y}}_{lm}\left(\mathit{\omega}\right)=\sqrt{(\mathrm{2}{\mathit{\delta}}_{\mathrm{0}m})(\mathrm{2}l+\mathrm{1}){\displaystyle \frac{(lm\left\right)\mathrm{!}}{(l+m\left\right)\mathrm{!}}}}\\ \text{(B2)}& {\displaystyle}& {\displaystyle}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{\mathcal{P}}_{l\leftm\right}\left(\mathrm{cos}\mathit{\theta}\right)\left\{\begin{array}{l}\mathrm{cos}\left(m\mathit{\varphi}\right)\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathrm{if}\phantom{\rule{0.25em}{0ex}}m\ge \mathrm{0}\\ \mathrm{sin}\left(\rightm\left\mathit{\varphi}\right)\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathrm{if}\phantom{\rule{0.25em}{0ex}}m\mathrm{0}\end{array}\right),\end{array}$$where δ_{0m} is the Kronecker delta. For $x\in \left[\mathrm{1},\mathrm{1}\right]$ and m≥0, polynomials 𝒫_{lm}(x) are given by
$${\mathcal{P}}_{lm}\left(x\right)=(\mathrm{1}{x}^{\mathrm{2}}{)}^{m/\mathrm{2}}{\displaystyle \frac{{\mathrm{d}}^{m}{\mathcal{P}}_{l}\left(x\right)}{\mathrm{d}{x}^{m}}},$$where
$$\begin{array}{}\text{(B3)}& {\displaystyle}{\mathcal{P}}_{l}\left(x\right)={\displaystyle \frac{\mathrm{1}}{{\mathrm{2}}^{l}\phantom{\rule{0.25em}{0ex}}l\mathrm{!}}}{\displaystyle \frac{{\mathrm{d}}^{l}({x}^{\mathrm{2}}\mathrm{1}{)}^{l}}{\mathrm{d}{x}^{l}}}\end{array}$$are the Legendre polynomials. This definition of SHs and their normalization are consistent with those employed for GRACE data generation and processing (Bettadpur, 2018) and can be evaluated straightforwardly, for example, using MATLAB's associated Legendre functions (https://www.mathworks.com/help/matlab/ref/legendre.html, last access: 7 May 2019).

SH addition theorem. It is useful to note here that the following relationship holds:
$$\begin{array}{}\text{(B4)}& {\displaystyle}{\mathcal{P}}_{l}\left(\mathrm{cos}\mathit{\alpha}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}l+\mathrm{1}}}\sum _{m=l}^{l}{\mathcal{Y}}_{lm}\left(\mathit{\omega}\right){\mathcal{Y}}_{lm}\left({\mathit{\omega}}^{\prime}\right),\end{array}$$where α once again is the greatcircle distance between coordinates ω and ω^{′}.

Evaluation of SH coefficients. For the chosen normalization, SHs obey the following orthogonality relationship
$$\begin{array}{}\text{(B5)}& {\displaystyle}\int {\mathcal{Y}}_{lm}\left(\mathit{\omega}\right){\mathcal{Y}}_{{l}^{\prime}{m}^{\prime}}\left(\mathit{\omega}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}=\mathrm{4}\mathit{\pi}{\mathit{\delta}}_{l{l}^{\prime}}{\mathit{\delta}}_{m{m}^{\prime}},\end{array}$$where ${\mathit{\delta}}_{l{l}^{\prime}}$ and ${\mathit{\delta}}_{m{m}^{\prime}}$ are Kronecker deltas. Using this property, SH coefficients of f(ω) are obtained as follows:
$$\begin{array}{}\text{(B6)}& {\displaystyle}{f}_{lm}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}}\int f\left(\mathit{\omega}\right){\mathcal{Y}}_{lm}\left(\mathit{\omega}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}.\end{array}$$ 
Evaluation of surface integrals on a unit sphere. We discretize the surface of a unit sphere using the socalled icosahedral pixelization method (Tegmark, 1996). It yields uniformly distributed quadrature points with equal pixel area. This makes numerical integration fairly straightforward as follows:
$$\begin{array}{}\text{(B7)}& {\displaystyle}\int f\left(\mathit{\omega}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}={\displaystyle \frac{\mathrm{4}\mathit{\pi}}{{N}_{\mathrm{T}}}}\phantom{\rule{0.25em}{0ex}}\sum _{j=\mathrm{1}}^{{N}_{\mathrm{T}}}f\left({\mathit{\omega}}_{j}\right),\end{array}$$where ω_{j} is the centroid of the jth pixel, and N_{T} is the total number of pixels. Note that the factor 4π∕N_{T} represents the area of each pixel on the surface of a unit sphere.
B2 SH coefficients of some basic functions

Ocean function. By definition, the ocean function is given by
$$\mathcal{O}\left(\mathit{\omega}\right)=\left\{\begin{array}{ll}\mathrm{1}& \mathrm{if}\phantom{\rule{0.25em}{0ex}}\mathit{\omega}\in {\mathbb{S}}_{\mathrm{O}}\\ \mathrm{0}& \mathrm{otherwise}\end{array}\right),$$where 𝕊_{O} is the ocean surface domain on a unit sphere. As in Eq. (B1), we may write
$$\mathcal{O}\left(\mathit{\omega}\right)=\sum _{lm}{\mathcal{O}}_{lm}{\mathcal{Y}}_{lm}\left(\mathit{\omega}\right).$$Following Eq. (B6) and using the definition of ocean function, we get
$$\begin{array}{}\text{(B8)}& {\displaystyle}{\mathcal{O}}_{lm}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}}\underset{{\mathbb{S}}_{\mathrm{O}}}{\int}{\mathcal{Y}}_{lm}\left(\mathit{\omega}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S},\end{array}$$where integration is performed only within the ocean surface domain. Following Eq. (B7), we obtain
$$\begin{array}{}\text{(B9)}& {\displaystyle}{\mathcal{O}}_{lm}={\displaystyle \frac{\mathrm{1}}{{N}_{\mathrm{T}}}}\sum _{j\in \phantom{\rule{0.125em}{0ex}}{\mathbb{S}}_{\mathrm{O}}}{\mathcal{Y}}_{lm}\left({\mathit{\omega}}_{j}\right).\end{array}$$ 
Model forcing function. By definition, $F\left(\mathit{\omega}\right)=H\left(\mathit{\omega}\right)\left[\mathrm{1}\mathcal{O}\left(\mathit{\omega}\right)\right]$. We may write
$${F}_{lm}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}}\underset{{\mathbb{S}}_{\mathrm{C}}}{\int}H\left(\mathit{\omega}\right){\mathcal{Y}}_{lm}\left(\mathit{\omega}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S},$$where 𝕊_{C} is the continental domain on a unit sphere. We derive H(ω) from the GRACE Stokes coefficients as detailed in Sect. 3. Following Eq. (B7), we get
$$\begin{array}{}\text{(B10)}& {\displaystyle}{F}_{lm}={\displaystyle \frac{\mathrm{1}}{{N}_{\mathrm{T}}}}\sum _{j\in \phantom{\rule{0.125em}{0ex}}{\mathbb{S}}_{\mathrm{C}}}H\left({\mathit{\omega}}_{j}\right){\mathcal{Y}}_{lm}\left({\mathit{\omega}}_{j}\right).\end{array}$$ 
Global surface loading function. Since $L=F+\widehat{S}$, we may write SH coefficients of L (Eq. 3) as follows:
$$\begin{array}{}\text{(B11)}& {\displaystyle}{L}_{lm}={F}_{lm}+{\widehat{S}}_{lm}.\end{array}$$
B3 Some useful integrals and barystatic sea level

Ocean surface area on a unit sphere. Since 𝒴_{00}(w)=1 (see Appendix B1), the SH coefficient of the ocean function (Eq. B9) for l=0 and m=0 is given by
$$\begin{array}{}\text{(B12)}& {\displaystyle}{\mathcal{O}}_{\mathrm{00}}={\displaystyle \frac{{N}_{\mathrm{O}}}{{N}_{\mathrm{T}}}}\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}}\left[\left({\displaystyle \frac{\mathrm{4}\mathit{\pi}}{{N}_{\mathrm{T}}}}\right){N}_{\mathrm{O}}\right],\end{array}$$where N_{O} is the number of pixels in the ocean surface domain 𝕊_{O}. Since the term enclosed by parentheses represents the area of each pixel, the total area of ocean surface (i.e., N_{O} times the pixel area) on a unit sphere is given by
$$\begin{array}{}\text{(B13)}& {\displaystyle}\int \mathcal{O}\left(\mathit{\omega}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}=\mathrm{4}\mathit{\pi}\phantom{\rule{0.25em}{0ex}}{\mathcal{O}}_{\mathrm{00}}.\end{array}$$ 
Continental water volume on a unit sphere. The SH coefficient of the forcing function (B10) for $l=m=\mathrm{0}$ is given by
$$\begin{array}{ll}{\displaystyle}{F}_{\mathrm{00}}& {\displaystyle}={\displaystyle \frac{\mathrm{1}}{{N}_{\mathrm{T}}}}\sum _{j\in \phantom{\rule{0.125em}{0ex}}{\mathbb{S}}_{\mathrm{C}}}H\left({\mathit{\omega}}_{j}\right)\\ \text{(B14)}& {\displaystyle}& {\displaystyle}\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}}\sum _{j\in \phantom{\rule{0.125em}{0ex}}{\mathbb{S}}_{\mathrm{C}}}\left[\left({\displaystyle \frac{\mathrm{4}\mathit{\pi}}{{N}_{\mathrm{T}}}}\right)H\left({\mathit{\omega}}_{j}\right)\right].\end{array}$$Since the term enclosed by parentheses represents the area of each pixel, the sum of the bracketed term over 𝕊_{C} essentially yields the total continental water volume. Consequently, we may write
$$\begin{array}{}\text{(B15)}& {\displaystyle}\int H\left(\mathit{\omega}\right)\left[\mathrm{1}\mathcal{O}\left(\mathit{\omega}\right)\right]\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}=\mathrm{4}\mathit{\pi}\phantom{\rule{0.25em}{0ex}}{F}_{\mathrm{00}}.\end{array}$$ 
Barystatic sea level. Using Eqs. (A1), (A2), (B13), and (B15), we get
$$\begin{array}{}\text{(B16)}& {\displaystyle}E={\displaystyle \frac{{F}_{\mathrm{00}}}{{\mathcal{O}}_{\mathrm{00}}}}.\end{array}$$
B4 SH coefficients appearing in the sealevel equation

The coefficient ${\widehat{E}}_{lm}$. This coefficient is used as the first guess solution of ${\widehat{S}}_{lm}$ (Eq. 8) and remains unchanged during the recursive process. Recalling that $\widehat{E}=\mathcal{O}E$ and that E is a spatial invariant, we may write
$${\widehat{E}}_{lm}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}}E\int \mathcal{O}\left(\mathit{\omega}\right){\mathcal{Y}}_{lm}\left(\mathit{\omega}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}.$$Noting that the integral is equivalent to 4π𝒪_{lm} (see Eq. B6), and using Eq. (B16), we get
$$\begin{array}{}\text{(B17)}& {\displaystyle}{\widehat{E}}_{lm}={\displaystyle \frac{{F}_{\mathrm{00}}}{{\mathcal{O}}_{\mathrm{00}}}}{\mathcal{O}}_{lm}.\end{array}$$ 
Other hatted coefficients. All of the coefficients appearing in Eq. (8) may be evaluated in a similar manner. Consider ${\widehat{X}}_{lm}$, for example. Recalling the definition that $\widehat{X}=\mathcal{O}X$, and following Eq. (B6), we may write
$${\widehat{X}}_{lm}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}}\int \mathcal{O}\left(\mathit{\omega}\right)X\left(\mathit{\omega}\right){\mathcal{Y}}_{lm}\left(\mathit{\omega}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}.$$Using the definition of ocean function, and expanding X(ω) as in Eq. (B1), we get
$${\widehat{X}}_{lm}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}}\underset{{\mathbb{S}}_{\mathrm{O}}}{\int}\left[\sum _{pq}{X}_{pq}{\mathcal{Y}}_{pq}\left(\mathit{\omega}\right)\right]{\mathcal{Y}}_{lm}\left(\mathit{\omega}\right)\phantom{\rule{0.25em}{0ex}}\mathrm{d}\mathit{\omega}.$$Following Eq. (B7), we evaluate the integral as follows:
$$\begin{array}{}\text{(B18)}& {\displaystyle}{\widehat{X}}_{lm}={\displaystyle \frac{\mathrm{1}}{{N}_{\mathrm{T}}}}\sum _{j\in \phantom{\rule{0.125em}{0ex}}{\mathbb{S}}_{\mathrm{O}}}\left[\sum _{pq}{X}_{pq}{\mathcal{Y}}_{pq}\left({\mathit{\omega}}_{j}\right)\right]{\mathcal{Y}}_{lm}\left({\mathit{\omega}}_{j}\right).\end{array}$$ 
The coefficient X_{pq}. By definition, $X={\mathrm{\Phi}}_{\mathrm{g}}/g$. Using Eqs. (A3) and (A4), we may write
$$X\left(\mathit{\omega}\right)={\displaystyle \frac{\mathrm{3}{\mathit{\rho}}_{\mathrm{w}}}{\mathrm{4}\mathit{\pi}{\mathit{\rho}}_{e}}}\int \sum _{p=\mathrm{0}}^{\mathrm{\infty}}\left(\mathrm{1}+{k}_{p}^{\prime}\right){\mathcal{P}}_{p}\left(\mathrm{cos}\mathit{\alpha}\right)\phantom{\rule{0.125em}{0ex}}L\left({\mathit{\omega}}^{\prime}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}{\mathbb{S}}^{\prime}.$$Using the SH addition theorem (Eq. B4), and expanding L(ω^{′}) as in Eq. (B1), we get
$$\begin{array}{ll}{\displaystyle}X\left(\mathit{\omega}\right)=& {\displaystyle \frac{\mathrm{3}{\mathit{\rho}}_{\mathrm{w}}}{\mathrm{4}\mathit{\pi}{\mathit{\rho}}_{e}}}\int \sum _{pq}\left\{{\displaystyle \frac{\mathrm{1}+{k}_{p}^{\prime}}{\mathrm{2}p+\mathrm{1}}}{\mathcal{Y}}_{pq}\left(\mathit{\omega}\right){\mathcal{Y}}_{pq}\left({\mathit{\omega}}^{\prime}\right)\right\}\\ {\displaystyle}& {\displaystyle}\sum _{{p}^{\prime}{q}^{\prime}}{L}_{{p}^{\prime}{q}^{\prime}}{\mathcal{Y}}_{{p}^{\prime}{q}^{\prime}}\left({\mathit{\omega}}^{\prime}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}{\mathbb{S}}^{\prime}\end{array}$$or
$$\begin{array}{ll}{\displaystyle}X\left(\mathit{\omega}\right)=& {\displaystyle \frac{\mathrm{3}{\mathit{\rho}}_{\mathrm{w}}}{\mathrm{4}\mathit{\pi}{\mathit{\rho}}_{e}}}\sum _{pq}\sum _{{p}^{\prime}{q}^{\prime}}{\displaystyle \frac{\mathrm{1}+{k}_{p}^{\prime}}{\mathrm{2}p+\mathrm{1}}}{L}_{{p}^{\prime}{q}^{\prime}}{\mathcal{Y}}_{pq}\left(\mathit{\omega}\right)\\ {\displaystyle}& {\displaystyle}\int {\mathcal{Y}}_{pq}\left({\mathit{\omega}}^{\prime}\right){\mathcal{Y}}_{{p}^{\prime}{q}^{\prime}}\left({\mathit{\omega}}^{\prime}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}{\mathbb{S}}^{\prime}.\end{array}$$Using the SH orthogonality relationship (Eq. B5), we get
$$\begin{array}{}\text{(B19)}& {\displaystyle}X\left(\mathit{\omega}\right)={\displaystyle \frac{\mathrm{3}{\mathit{\rho}}_{\mathrm{w}}}{{\mathit{\rho}}_{e}}}\sum _{pq}{\displaystyle \frac{\mathrm{1}+{k}_{p}^{\prime}}{\mathrm{2}p+\mathrm{1}}}{L}_{pq}{\mathcal{Y}}_{pq}\left(\mathit{\omega}\right).\end{array}$$Using Eqs. (B6) and (B19), SH coefficients X_{pq} are given by
$${X}_{pq}={\displaystyle \frac{\mathrm{3}{\mathit{\rho}}_{\mathrm{w}}}{\mathrm{4}\mathit{\pi}{\mathit{\rho}}_{e}}}\int \sum _{{p}^{\prime}{q}^{\prime}}\left\{{\displaystyle \frac{\mathrm{1}+{k}_{{p}^{\prime}}^{\prime}}{\mathrm{2}{p}^{\prime}+\mathrm{1}}}{L}_{{p}^{\prime}{q}^{\prime}}{\mathcal{Y}}_{{p}^{\prime}{q}^{\prime}}\left(\mathit{\omega}\right)\right\}{\mathcal{Y}}_{pq}\left(\mathit{\omega}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}.$$Rearranging terms and applying the orthogonality relationship (Eq. B5), we obtain
$$\begin{array}{}\text{(B20)}& {\displaystyle}{X}_{pq}={\displaystyle \frac{\mathrm{3}{\mathit{\rho}}_{\mathrm{w}}}{{\mathit{\rho}}_{e}}}{\displaystyle \frac{\mathrm{1}+{k}_{p}^{\prime}}{\mathrm{2}p+\mathrm{1}}}{L}_{pq}.\end{array}$$ 
The coefficient Y_{pq}. By definition, $Y={\mathrm{\Phi}}_{\mathrm{r}}/g$. Using Eqs. (A5) and (B6), SH coefficients Y_{pq} are given by
$$\begin{array}{ll}{\displaystyle}{Y}_{pq}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}g}}\int \left[{\mathrm{\Lambda}}_{\mathrm{00}}{\mathcal{Y}}_{\mathrm{00}}\left(\mathit{\omega}\right)+\sum _{\mathrm{2}{q}^{\prime}}(\mathrm{1}+{k}_{\mathrm{2}}){\mathrm{\Lambda}}_{\mathrm{2}{q}^{\prime}}{\mathcal{Y}}_{\mathrm{2}{q}^{\prime}}\left(\mathit{\omega}\right)\right]\\ {\displaystyle}& {\displaystyle}{\mathcal{Y}}_{pq}\left(\mathit{\omega}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}.\end{array}$$Rearranging terms and applying the orthogonality relationship (Eq. B5), we get
$$\begin{array}{}\text{(B21)}& {\displaystyle}{Y}_{pq}={\displaystyle \frac{{\mathrm{\Lambda}}_{\mathrm{00}}}{g}}{\mathit{\delta}}_{p\mathrm{0}}{\mathit{\delta}}_{q\mathrm{0}}+{\displaystyle \frac{(\mathrm{1}+{k}_{\mathrm{2}}){\mathrm{\Lambda}}_{\mathrm{2}q}}{g}}{\mathit{\delta}}_{p\mathrm{2}}.\end{array}$$ 
The coefficient P_{pq}. By definition, $P={U}_{\mathrm{g}}$. Using Eqs. (A3) and (A4), and following the procedure to derive X_{pq} (Eq. B20), we obtain
$$\begin{array}{}\text{(B22)}& {\displaystyle}{P}_{pq}={\displaystyle \frac{\mathrm{3}{\mathit{\rho}}_{\mathrm{w}}}{{\mathit{\rho}}_{e}}}{\displaystyle \frac{{h}_{p}^{\prime}}{\mathrm{2}p+\mathrm{1}}}{L}_{pq}.\end{array}$$ 
The coefficient Q_{pq}. By definition, $Q={U}_{\mathrm{r}}$. Using Eq. (A5), and following the procedure to derive Y_{pq} (Eq. B21), we get
$$\begin{array}{}\text{(B23)}& {\displaystyle}{Q}_{pq}={\displaystyle \frac{{h}_{\mathrm{2}}{\mathrm{\Lambda}}_{\mathrm{2}q}}{g}}{\mathit{\delta}}_{p\mathrm{2}}.\end{array}$$ 
The coefficient C_{pq}. By definition, $C=E\u2329\mathrm{\Phi}/gU\u232a\equiv E\u2329X+P+Y+Q\u232a$. Using Eqs. (A5) and (B19) and similar equations for P, we may write
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}C=E{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}{\mathcal{O}}_{\mathrm{00}}}}\underset{{\mathbb{S}}_{\mathrm{O}}}{\int}\left(\right)open="[">{\displaystyle \frac{\mathrm{3}{\mathit{\rho}}_{\mathrm{w}}}{{\mathit{\rho}}_{e}}}\sum _{pq}{\displaystyle \frac{\mathrm{1}+{k}_{p}^{\prime}{h}_{p}^{\prime}}{\mathrm{2}p+\mathrm{1}}}{L}_{pq}{\mathcal{Y}}_{pq}\left(\mathit{\omega}\right)\end{array}{\displaystyle}& {\displaystyle}\left(\right)close="]">+{\displaystyle \frac{\mathrm{1}}{g}}\left\{{\mathrm{\Lambda}}_{\mathrm{00}}{\mathcal{Y}}_{\mathrm{00}}\left(\mathit{\omega}\right)+\sum _{\mathrm{2}q}\left(\mathrm{1}+{k}_{\mathrm{2}}{h}_{\mathrm{2}}\right){\mathrm{\Lambda}}_{\mathrm{2}q}{\mathcal{Y}}_{\mathrm{2}q}\left(\mathit{\omega}\right)\right\}\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}.$$Using Eqs. (B8) and (B16), we get
$$\begin{array}{ll}{\displaystyle}C=& {\displaystyle}{\displaystyle \frac{{F}_{\mathrm{00}}}{{\mathcal{O}}_{\mathrm{00}}}}{\displaystyle \frac{\mathrm{3}{\mathit{\rho}}_{\mathrm{w}}}{{\mathit{\rho}}_{e}{\mathcal{O}}_{\mathrm{00}}}}\sum _{pq}{\displaystyle \frac{\mathrm{1}+{k}_{p}^{\prime}{h}_{p}^{\prime}}{\mathrm{2}p+\mathrm{1}}}{L}_{pq}{\mathcal{O}}_{pq}\\ \text{(B24)}& {\displaystyle}& {\displaystyle}{\displaystyle \frac{{\mathrm{\Lambda}}_{\mathrm{00}}}{g}}{\displaystyle \frac{\mathrm{1}}{g{\mathcal{O}}_{\mathrm{00}}}}\sum _{\mathrm{2}q}(\mathrm{1}+{k}_{\mathrm{2}}{h}_{\mathrm{2}}){\mathrm{\Lambda}}_{\mathrm{2}q}{\mathcal{O}}_{\mathrm{2}q}.\end{array}$$Note that C and all of the righthand side terms are spatially invariant. Using Eqs. (B6) and (B24), we get
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}{C}_{pq}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi}}}\left(\right)open="[">{\displaystyle \frac{{F}_{\mathrm{00}}}{{\mathcal{O}}_{\mathrm{00}}}}+{\displaystyle \frac{\mathrm{3}{\mathit{\rho}}_{\mathrm{w}}}{{\mathit{\rho}}_{e}{\mathcal{O}}_{\mathrm{00}}}}\sum _{{p}^{\prime}{q}^{\prime}}{\displaystyle \frac{\mathrm{1}+{k}_{{p}^{\prime}}^{\prime}{h}_{{p}^{\prime}}^{\prime}}{\mathrm{2}{p}^{\prime}+\mathrm{1}}}{L}_{{p}^{\prime}{q}^{\prime}}{\mathcal{O}}_{{p}^{\prime}{q}^{\prime}}\end{array}{\displaystyle}& {\displaystyle}\left(+{\displaystyle \frac{{\mathrm{\Lambda}}_{\mathrm{00}}}{g}}+{\displaystyle \frac{\mathrm{1}}{g{\mathcal{O}}_{\mathrm{00}}}}\sum _{\mathrm{2}{q}^{\prime}}(\mathrm{1}+{k}_{\mathrm{2}}{h}_{\mathrm{2}}){\mathrm{\Lambda}}_{\mathrm{2}{q}^{\prime}}{\mathcal{O}}_{\mathrm{2}{q}^{\prime}}\right]\\ {\displaystyle}& {\displaystyle}\int {\mathcal{Y}}_{pq}\left(\mathit{\omega}\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathbb{S}.$$Since 𝒴_{00}(ω)=1, we introduce a virtual expression 𝒴_{00}(ω) inside the integral and use Eq. (B5) to find
$$\begin{array}{ll}{\displaystyle}& {\displaystyle}{C}_{pq}={\mathit{\delta}}_{p\mathrm{0}}{\mathit{\delta}}_{{q}_{\mathrm{0}}}\left[{\displaystyle \frac{{F}_{\mathrm{00}}}{{\mathcal{O}}_{\mathrm{00}}}}+{\displaystyle \frac{\mathrm{3}{\mathit{\rho}}_{\mathrm{w}}}{{\mathit{\rho}}_{e}{\mathcal{O}}_{\mathrm{00}}}}\sum _{{p}^{\prime}{q}^{\prime}}{\displaystyle \frac{\mathrm{1}+{k}_{{p}^{\prime}}^{\prime}{h}_{{p}^{\prime}}^{\prime}}{\mathrm{2}{p}^{\prime}+\mathrm{1}}}{L}_{{p}^{\prime}{q}^{\prime}}{\mathcal{O}}_{{p}^{\prime}{q}^{\prime}}\right)\\ \text{(B25)}& {\displaystyle}& {\displaystyle}\left(\right)close="]">+\left\{{\displaystyle \frac{{\mathrm{\Lambda}}_{\mathrm{00}}}{g}}+{\displaystyle \frac{\mathrm{1}}{g{\mathcal{O}}_{\mathrm{00}}}}\sum _{\mathrm{2}{q}^{\prime}}(\mathrm{1}+{k}_{\mathrm{2}}{h}_{\mathrm{2}}){\mathrm{\Lambda}}_{\mathrm{2}{q}^{\prime}}{\mathcal{O}}_{\mathrm{2}{q}^{\prime}}\right\}.\end{array}$$The terms inside the braces vanish when rotational feedback is not included. Note that the first term on the righthand side of Eq. (B21) and the first term inside the braces above cancel out while solving Eq. (9). It is, however, important to consider these terms explicitly for a clean isolation of SH coefficients of the desired fields.
B5 Summary
Here we briefly outline the workflow of our computation.

Given the onland change in water equivalent height $H(\mathit{\theta},\mathit{\varphi},t)$ over the time period t, we compute SH coefficients of the model forcing function F_{lm} using Eq. (B10).

We compute SH coefficients of the “global loading function” L_{lm} (Eq. B11) by initializing ${\widehat{S}}_{lm}$, such that ${\widehat{S}}_{lm}\equiv {\widehat{S}}_{lm}^{\mathrm{0}}={\widehat{E}}_{lm}$, where ${\widehat{E}}_{lm}$ is given by Eq. (B17).

Once ${\widehat{S}}_{lm}$, and hence L_{lm}, are initialized, we solve the recursion equation for ${\widehat{S}}_{lm}^{r}$ (Eq. 8) until the solution is converged. The hatted SH coefficients appearing in Eq. (8) are expressed in terms of their nonhatted counterparts as in, for example, Eq. (B18). For the chosen solid Earth model and the reference frame origin, we compute these nonhatted coefficients using Eqs. (B20)–(B23) and (B25).

The choice of the solid Earth model determines the load Love numbers ${k}_{p}^{\prime}$ and ${h}_{p}^{\prime}$ and the tidal Love numbers k_{2} and h_{2}. In this study, we consider the Preliminary Reference Earth Model (Dziewonski and Anderson, 1981).

The choice of the reference frame origin determines the degree 1 load Love numbers. In this study, we take the values from Blewitt (2003) for the CM and CF reference frames.

Rotational feedback is accounted for via Eqs. (B21), (B23), and (B25); SH coefficients of perturbation in rotational potential Λ_{pq} appearing in these equations are given by Eqs. (A6) and (A7). To deactivate this feedback, we set Y_{pq}=0 (Eq. B21) and Q_{pq}=0 (Eq. B23) and remove the terms enclosed by braces in Eq. (B25).

As for the convergence criterion, we track the relative change in L2 norm after each recursion and call the solution converged when it is less than 0.001 % of the L2 norm of the solution itself. This level of solution convergence is typically achieved after four to six iterations.

Unless stated otherwise, constants and parameters used in this study are taken directly from Table 1 of Adhikari et al. (2016).


Once the solution is converged for ${\widehat{S}}_{lm}$, we appropriately combine corresponding nonhatted coefficients (i.e., Eqs. B20–B23 and B25) in order to retrieve SH coefficients of relative sea level S_{pq} (Eq. 8), geoid height change [X_{pq}+Y_{pq}], and vertical bedrock motion $[{P}_{pq}{Q}_{pq}]$. We supply these coefficients along with corresponding spatial maps.
SA and ERI conceived the research and wrote the first draft of the paper. SA formulated the sealevel solver, with the help of ERI, and led the calculations, with the help of FWL (in estimating degree 1 Stokes coefficients) and TF and LC (in uncertainty quantification). All authors contributed to the analysis of the results and to the writing and editing of the paper.
The authors declare that they have no conflict of interest.
This research was carried out at the Jet Propulsion Laboratory (JPL), California Institute of Technology, under a contract with National Aeronautics and Space Administration (NASA). Constructive comments from the reviewers greatly improved the quality of the paper.
This research was supported by the JPL Research, Technology & Development Program (grant no. 01STCRR.17.235.118; 2017–2019) and the NASA SeaLevel Change Science Team (grant no. 16SLCT160015; 2018–2020), as well as the NASA MEaSUREs program (grant no. 17MEASURES0031; 2019–2023).
This paper was edited by Giuseppe M. R. Manzella and reviewed by Don Chambers, Xuebin Zhang, and Makan A. Karegar.
Adhikari, S. and Ivins, E. R.: Climatedriven polar motion: 2003–2015, Sci. Adv., 2, e1501693, https://doi.org/10.1126/sciadv.1501693, 2016. a, b, c
Adhikari, S., Ivins, E. R., and Larour, E.: ISSMSESAW v1.0: meshbased computation of gravitationally consistent sealevel and geodetic signatures caused by cryosphere and climate driven mass change, Geosci. Model Dev., 9, 1087–1109, https://doi.org/10.5194/gmd910872016, 2016. a, b, c, d
Adhikari, S., Ivins, E. R., Frederikse, T., Landerer, F. W., and Caron, L.: Changes in relative sealevel, geoid height, and bedrock displacement derived from the Release06 GRACE Level2 monthly Stokes coefficients, Harvard Dataverse, V2, https://doi.org/10.7910/DVN/8UC8IR, 2019. a, b
Bettadpur, S.: Gravity Recovery and Climate Experiment Level – 2 Gravity Field Product User Handbook, GRACE 327734 (CSRG0301) (Rev 4.0, 25 April 2018), available at: ftp://podaac.jpl.nasa.gov/allData/grace/docs/L2UserHandbook_v4.0.pdf (last access: 7 May 2019), 2018. a, b, c
Bilbao, R. A., Gregory, J. M., and Bouttes, N.: Analysis of the regional pattern of sea level change due to ocean dynamics and density change for 1993–2099 in observations and CMIP5 AOGCMs, Clim. Dynam., 45, 2647–2666, https://doi.org/10.1007/s003820152499z, 2015. a
Blewitt, G.: Selfconsistency in reference frames, geocenter definition, and surface loading of the solid Earth, J. Geophys. Res., 108, B22103, https://doi.org/10.1029/2002JB002082, 2003. a
Caron, L., Ivins, E. R., Larour, E., Adhikari, S., Nilsson, J., and Blewitt, G.: GIA model statistics for GRACE hydrology, cryosphere and ocean science, Geophys. Res. Lett., 45, 2203–2212, https://doi.org/10.1002/2017GL076644, 2018. a, b, c, d
Chambers, D. P.: Evaluation of new GRACE timevariable gravity data over the ocean, Geophys. Res. Lett., 33, L17603, https://doi.org/10.1029/2006GL027296, 2006. a
Chambers, D. P. and Willis, J. K.: Analysis of largescale ocean bottom pressure variability in the North Pacific, J. Geophys. Res., 113, C11003, https://doi.org/10.1029/2008JC004930, 2008. a
Cheng, M. K. and Ries, J. C.: The unexpected signal in GRACE estimates of C20, J. Geod., 91, 897–914, https://doi.org/10.1007/s0019001609955, 2017. a
Cheng, M. K., Ries, J. C., and Tapley, B. D.: Variations of the Earth's figure axis from satellite laser ranging and GRACE, J. Geophys. Res., 116, B01409, https://doi.org/10.1029/2010JB000850, 2011. a
Clark, J. A., Haidle, P. E., and Cunningham, N. L.: Comparison of satellite altimetry to tide gauge measurement of sea level: Predictions of glacioisostatic adjustment, J. Climate, 15, 22, 3291–3300, https://doi.org/10.1175/15200442(2002)015<3291:COSATT>2.0.CO;2, 2002. a
Davis, J. L. and Vinogradova, N. T.: Causes of accelerating sealevel on the East Coast of North America, Geophys. Res. Lett., 44, 5133–5141, https://doi.org/10.1002/2017GL072845, 2017. a
Desai, S. D.: Observing the pole tide with satellite altimetry, J. Geophys. Res., 107, 3186, https://doi.org/10.1029/2001JC001224, 2002. a
Desai, S. D. and Yuan, D.N.: Application of the convolution formalism to the ocean tide potential: Results from the Gravity Recovery and Climate Experiment (GRACE), J. Geophys. Res., 111, C06023, https://doi.org/10.1029/2005JC003361, 2006. a
Desai, S. D., Wahr, J., and Beckley, B.: Revisiting the pole tide for and from satellite altimetry, J. Geod., 89, 1233–1243, https://doi.org/10.1007/s0019001508487, 2015. a, b
Dieng, H. B., Cazenave, A., von Schuckmann, K., Ablain, M., and Meyssignac, B.: Sea level budget over 2005–2013: missing contributions and data errors, Ocean Sci., 11, 789–802, https://doi.org/10.5194/os117892015, 2015. a
Dobslaw, H., BergmannWolf, I., Dill, R., Poropat, L., Thomas, M., Dahle, C., Esselborn, S., Konig, R., and Flechtner, F.: A new highresolution model of nontidal atmosphere and ocean mass variability for dealiasing of satellite gravity observations: AOD1B RL06, Geophys. J. Int., 211, 263–269, https://doi.org/10.1093/gji/ggx302, 2017. a
Dziewonski, A. M. and Anderson, D. L.: Preliminary reference Earth model, Earth Planet. Inter., 25, 297–356, https://doi.org/10.1016/00319201(81)900467, 1981. a
Farrell, W. E. and Clark, J. A.: On postglacial sea level, Geophys. J. Roy. Astr. S., 46, 647–667, https://doi.org/10.1111/j.1365246X.1976.tb01252.x, 1976. a
Fasullo, J. T. and Nerem, R. S.: Altimeterera emergence of the patterns of forced sealevel rise in climate models and implications for the future, P. Natl. Acad. Sci. USA, 115, 12944–12949, https://doi.org/10.1073/pnas.1813233115, 2018. a
Frederikse, T., Jevrejeva, S., Riva, R. E. M., and Dangendorf, S.: A consistent sealevel reconstruction and its budget on basin and global scales over 1958–2014, J. Climate, 31, 1267–1280, https://doi.org/10.1175/JCLID170502.1, 2018. a
Galassi, G. and Spada, G.: Tide gauge observations in Antarctica (1958–2014) and recent ice loss, Antarct. Sci., 29, 369–381, https://doi.org/10.1017/S0954102016000729, 2017. a
Gaultier, L., Ubelmann, C. and Fu, L.L.: The challenge of using future SWOT data for oceanic field reconstruction, J. Atmos. Ocean. Tech., 33, 119–126, https://doi.org/10.1175/JTECHD150160.1, 2016. a
Han, S.C., Sauber, J., and Pollitz, F.: Coseismic compression/dilatation and viscoelastic uplift/subsidence following the 2012 Indian Ocean earthquakes quantified from satellite gravity observations, Geophys. Res. Lett., 42, 3764–3722, https://doi.org/10.1002/2015GL063819, 2015. a
Han, S.C., Sauber, J., and Pollitz, F.: Postseismic gravity change after the 2006–2007 great earthquake doublet and constraints on the asthenosphere structure in the central Kuril Islands, Geophys. Res. Lett., 43, 3169–3177, https://doi.org/10.1002/2016GL068167, 2016. a
Haubrich, R. and Munk, W.: The pole tide, J. Geophys. Res., 64, 2373–2388, https://doi.org/10.1029/JZ064i012p02373, 1959. a
Hirschi, M. and Seneviratne, S. I.: Basinscale waterbalance dataset (BSWB): an update, Earth Syst. Sci. Data, 9, 251–258, https://doi.org/10.5194/essd92512017, 2017. a
Hsu, C.W. and Velicogna, I.: Detection of sea level fingerprints derived from GRACE gravity data, Geophys. Res. Lett., 44, 8953–8961, https://doi.org/10.1002/2017GL074070, 2017. a
Ivins, E. R., Watkins, M. M., Yuan, D.N., Dietrich, R., Casassa, G., and Rülke, A.: Onland ice loss and glacial isostatic adjustment at the Drake Passage: 2003–2009, J. Geophys. Res., 116, B02403, https://doi.org/10.1029/2010JB007607, 2011. a
Ivins, E. R., James, T. S., Wahr, J., Schrama, E. J. O., Landerer, F. W., and Simon, K. M.: Antarctic contribution to sealevel rise observed by GRACE with improved GIA correction, J. Geophys. Res., 118, 3126–3141, https://doi.org/10.1002/jgrb.50208, 2013. a
Jacob, T., Wahr, J., Pfeffer, W. T., and Swenson, S.: Recent contributions of glaciers and ice caps to sea level rise, Nature, 482, 514–518, https://doi.org/10.1038/nature10847, 2012. a
Johnson, G. F. and Chambers, D. P.: Ocean bottom pressure seasonal cycles and decadal trends from GRACE release05: Ocean circulation implications, J. Geophys. Res., 118, 4228–4240, https://doi.org/10.1002/jgrc.20307, 2013. a
King, M. A. and SantamaríaGómez, A.: Ongoing deformation of Antarctica following recent Great Earthquakes, Geophys. Res. Lett., 43, 1918–1927, https://doi.org/10.1002/2016GL067773, 2016. a
Lambeck, K.: The Earth's Variable Rotation: Geophysical Causes and Consequnces, Cambridge Univ. Press, Cambridge, UK, 449 pp., 1981. a, b
Lambert, A. and Beaumont, C.: Nanovariations in gravity due to seasonal groundwater movements: implications for gravitational detection of tectonic movements, J. Geophys. Res., 82, 297–306, https://doi.org/10.1029/JB082i002p00297, 1977. a
Landerer, F. W. and Swenson, S. C.: Accuracy of scaled GRACE terrestrial water storage estimates, Water Resour. Res., 48, W04531, https://doi.org/10.1029/2011WR011453, 2012. a
Landerer, F. W., Wiese, D. N., Bentel, K., Böning, C., and Watkins, M. M.: North Atlantic meridional overturning circulation variations from GRACE ocean bottom pressure anomalies, Geophys. Res. Lett., 42, 8114–8121, https://doi.org/10.1002/2015GL065730, 2015. a
Lange, H., Casassa, G., Ivins, E. R., Schröder, L., Fritsche, M., Richter, A., Groh, A., and Dietrich, R.: Observed crustal uplift near the Southern Patagonian Icefield constrains improved viscoelastic Earth models, Geophys. Res. Lett., 41, 805–812, https://doi.org/10.1002/2013GL058419, 2014. a
Larour, E., Ivins, E. R., and Adhikari, S.: Should coastal planners have concern over where land ice is melting?, Sci. Adv., 3, e1700537, https://doi.org/10.1126/sciadv.1700537, 2017. a
Lickley, M. J., Hay, C. C., Tamisiea, M. E., and Mitrovica, J. X.: Bias in estimates of global mean sea level change inferred from satellite altimetry, J. Climate, 31, 5263–5271, https://doi.org/10.1175/JCLID180024.1, 2018. a
Luthcke, S. B., Sabaka, T. J., Loomis, B. D., Arendt, A. A., McCarthy, J. J., and Camp, J.: Antarctica, Greenland and Gulf of Alaska landice evolution from an iterated GRACE global mascon solution, J. Glaciol., 59, 613–631, https://doi.org/10.3189/2013JoG12J147, 2013. a
Mazloff, M. R. and Böning, C.: Rapid variability of Antarctic Bottom Water transport into the Pacific Ocean inferred from GRACE, Geophys. Res. Lett., 43, 3822–3829, https://doi.org/10.1002/2016gl068474, 2016. a
Milne, G. A. and Mitrovica, J.: Postglacial sealevel change on a rotating Earth, Geophys. J. Int., 133, 1–19, https://doi.org/10.1046/j.1365246X.1998.1331455.x, 1998. a, b
Mitrovica, J. X. and Peltier, W. R.: On postglacial geoid subsidence over the equatorial ocean, J. Geophys. Res., 96, 20053–20071, https://doi.org/10.1029/91JB01284, 1991. a, b
Mitrovica, J. X., Tamisiea, M. E., Davis, J. L., and Milne, G. A.: Recent mass balance of polar ice sheets inferred from patterns of global sealevel change, Nature, 409, 1026–1029, https://doi.org/10.1038/35059054, 2001. a
Mitrovica, J. X., Hay, C. C., Kopp, R. E., Harig, C., and Latychev, K.: Quantifying the sensitivity of sea level change in coastal localities to the geometry of polar ice mass flux, J. Climate, 31, 3701–3709, https://doi.org/10.1175/JCLID170465.1, 2018. a, b
Müller, F. L., Wekerle, C., Dettmering, D., Passaro, M., Bosch, W., and Seitz, F.: Dynamic ocean topography of the northern Nordic seas: a comparison between satellite altimetry and ocean modeling, The Cryosphere, 13, 611–626, https://doi.org/10.5194/tc136112019, 2019. a
Nield, G. A., Whitehouse, P. L., van der Wal, W., Blank, B., O'Donnell, J. P., and Stuart, G. W.: The impact of lateral variations in lithospheric thickness on glacial isostatic adjustment in West Antarctica, Geophys. J. Int., 214, 811–824, https://doi.org/10.1093/gji/ggy158, 2018. a
Peltier, W. R.: Dynamics of the Ice Age Earth, Adv. Geophys., 24, 1–146, https://doi.org/10.1016/S00652687(08)605191, 1982. a
Piecuch, C. G., Thompson, P. R., and Donohue, K. A.: Air pressure effects on sea level changes during the twentieth century, J. Geophys. Res., 121, 7917–7930, https://doi.org/10.1002/2016JC012131, 2016. a
Rapp, R. H.: The determination of geoid undulations and gravity anomalies from SEASAT altimeter data, J. Geophys. Res., 88, 1552–1562, https://doi.org/10.1029/JC088iC03p01552, 1983. a
Rietbroek, R., Brunnabend, S.E., Kusche, J., and Schröter, J.: Resolving sea level contributions by identifying fingerprints in timevariable gravity and altimetry, J. Geodynam., 59–60, 72–81, https://doi.org/10.1016/j.jog.2011.06.007, 2012. a
Rignot, E., Velicogna, I., van den Broeke, M. R., Monaghan, A., and Lenaerts, J.: Acceleration of the contribution of the Greenland and Antarctic ice sheets to sea level rise, Geophys. Res. Lett., 38, L05503, https://doi.org/10.1029/2011GL046583, 2011. a
Riva, R. E. M., Bamber, J. L., Lavallée, D. A., and Wouters, B.: Sealevel fingerprint of continental water and ice mass change from GRACE, Geophys. Res. Lett., 37, L19605, https://doi.org/10.1029/2010GL044770, 2010. a
Rodell, M., Beaudoing, H. K., L'Ecuyer, T. S., Olson, W. S., Famiglietti, J. S., Houser, P. R., Adler, R., Bosilovich, M. G., Clayson, C. A., Chambers, D., Clark, E., Fetzer, E. J., Gao, X., Gu, G., Hilburn, K., Huffman, G. J., Lettenmaier, D. P., Liu, W. T., Robertson, F. R., Schlosser, C. A., Sheffield, J., and Wood, E. F.: The observed state of the water cycle in the early TwentyFirst Century, J. Climate, 28, 8289–8317, https://doi.org/10.1175/JCLID1400555.1, 2015. a
Save, H., Bettadpur, S., and Tapley, B. D.: Highresolution CSR GRACE RL05 mascons, J. Geophys. Res., 121, 7547–7569, https://doi.org/10.1002/2016JB013007, 2016. a
Saynisch, J., BergmannWolf, I., and Thomas, M.: Assimilation of GRACEderived oceanic mass distributions with a global ocean circulation model, J. Geodesy, 89, 121–139, https://doi.org/10.1007/s0019001407660, 2015. a
Schrama, E. J. O., Wouters, B., and Rietbroek, R.: A mascon approach to assess ice sheet and glacier mass balances and their uncertainties from GRACE data, J. Geophys. Res., 119, 6048–6066, https://doi.org/10.1002/2013JB010923, 2014. a, b
Shepherd A. and the IMBIE2 Team: Mass balance of the Antarctic Ice Sheet from 1992 to 2017, Nature, 558, 219–222, https://doi.org/10.1038/s415860180179y, 2018. a
Simon, K. M., Riva, R. E. M., Kleinherenbrink, M., and Frederikse, T.: The glacial isostatic adjustment signal at present day in northern Europe and the British Isles estimated from geodetic observations and geophysical models, Solid Earth, 9, 777–795, https://doi.org/10.5194/se97772018, 2018. a
Spada, G.: Glacial isostatic adjustment and contemporary sea level rise: An overview, Surv. Geophys., 38, 153–185, https://doi.org/10.1007/s107120169379x, 2017. a
Spada, G. and Galassi, G.: Spectral analysis of sea level during the altimetry era, and evidence for GIA and glacial melting fingerprints, Global Planet. Change, 143, 34–49, https://doi.org/10.1016/j.gloplacha.2016.05.006, 2016. a
Stammer, D., Ray, R. D., Andersen, O. B., Arbic, B. K., Bosch, W., Carrère, L., Cheng, Y., Chinn, D. S., Dushaw, B. D., Egbert, G. D., Erofeeva, S. Y., Fok, H. S., Green, J. A. M., Griffiths, S., King, M. A., Lapin, V., Lemoine, F. G., Luthcke, S. B., Lyard, F., Morison, J., Müller, M., Padman, L., Richman, J. G., Shriver, J. F., Shum, C. K., Taguchi, E., and Yi, Y.: Accuracy assessment of global barotropic ocean tide models, Rev. Geophys., 52, 243–282, https://doi.org/10.1002/2014RG000450, 2014. a
Sterenborg,, M. G., Morrow, E., and Mitrovica, J. X.: Bias in GRACE estimates of ice mass change due to accompanying sealevel change, J. Geodesy, 87, 387–392, https://doi.org/10.1007/s001900120608x, 2013. a
Stocchi, P. and Spada, G.: Influence of glacial isostatic adjustment upon current sea level variations in the Mediterranean, Tectonophysics, 474, 56–68, https://doi.org/10.1016/j.tecto.2009.01.003, 2009. a
Swenson, S., Chambers, D., and Wahr, J.: Estimating geocenter variations from a combination of GRACE and ocean model output, J. Geophys. Res., 113, B08410, https://doi.org/10.1029/2007JB005338, 2008. a
Tamisiea, M. E.: Ongoing glacial isostatic contributions to observations of sea level change, Geophys. J. Int., 186, 1036–1044, https://doi.org/10.1111/j.1365246X.2011.05116.x, 2011. a, b
Tamisiea, M. E., Hughes, C. W., Williams, S. D. P., and Bingley, R. M.: Sea level: measuring the bounding surfaces of the ocean, Philos. T. Roy. Soc. A, 372, 20130336, https://doi.org/10.1098/rsta.2013.0336, 2014. a
Tapley, B. D., Watkins, M. M., Flechtner, F., Reigber, C., Bettadpur, S., Rodell, M., Sasgen, I., Famiglietti, J. S., Landerer, F. W., Chambers, D. P., Reager, J. T., Gardner, A. S., Save, H., Ivins, E. R., Swenson, S. C., Boening, C., Dahle, C., Wiese, D. N., Dobslaw, H., Tamisiea, M. E., and Velicogna, I.: Contributions of GRACE to understanding climate change, Nat. Clim. Change, 9, 358–369, https://doi.org/10.1038/s4155801904562, 2019. a
Tegmark, M.: An icosahedronbased method for pixelizing the celestial sphere, Astrophysical J. Lett., 470, L81–L84, https://doi.org/10.1086/310310, 1996. a
Velicogna, I.: Increasing rates of ice mass loss from the Greenland and Antarctic ice sheets revealed by GRACE, Geophys. Res. Lett., 36, L19503, https://doi.org/10.1029/2009GL040222, 2009. a
Wahr, J., Molenaar, M., and Bryan, F.: Time variability of the Earth's gravity field: Hydrological and oceanic effects and their possible detection using GRACE, J. Geophys. Res., 103, 30205–30230, https://doi.org/10.1029/98JB02844, 1998. a
Wahr, J., Nerem, R. S., and Bettadpur, S. V.: The pole tide and its effect on GRACE timevariable gravity measurements: Implications for estimates of surface mass variations, J. Geophys. Res.Sol. Ea., 120, 4597–4615, https://doi.org/10.1002/2015JB011986, 2015. a
Wahr, J. M.: Deformation induced by polar motion, J. Geophys. Res., 90, 9363–9368, https://doi.org/10.1029/JB090iB11p09363, 1985. a
Watkins, M. M., Wiese, D. N., Yuan, D.N., Böning, C., and Landerer, F. W.: Improved methods for observing Earth's time variable mass distribution with GRACE, J. Geophys. Res.Sol. Ea., 120, 2648–2671, https://doi.org/10.1002/2014JB011547, 2015. a, b, c
WCRP Global Sea Level Budget Group: Global sealevel budget 1993–present, Earth Syst. Sci. Data, 10, 1551–1590, https://doi.org/10.5194/essd1015512018, 2018. a
Wiese, D. N., Landerer, F. W., and Watkins, M. M.: Quantifying and reducing leakage errors in the JPL RL05M GRACE mascon solution, Water Resour. Res., 52, 7490–7502, https://doi.org/10.1002/2016WR019344, 2016. a, b
 Abstract
 Copyright statement
 Introduction
 Key variables and deliverables
 The sealevel equation
 GRACE and sealevel fingerprints
 Discussion
 Data availability
 Conclusions
 Appendix A: Theory of the sealevel fingerprint
 Appendix B: Spectral methods for sealevel equation
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Copyright statement
 Introduction
 Key variables and deliverables
 The sealevel equation
 GRACE and sealevel fingerprints
 Discussion
 Data availability
 Conclusions
 Appendix A: Theory of the sealevel fingerprint
 Appendix B: Spectral methods for sealevel equation
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References