Sea-level fingerprints emergent from GRACE mission data

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 climate-related changes on decadal timescales, such as ice mass loss and sea-level rise. As the changes in the latter are significant over the past two decades, there is a concomitant self-attraction 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 post-seismic viscoelastic mantle flow, or GIA, occurring beneath the seafloor. Here we perform a systematic calculation of sea-level fingerprints of on-land water mass changes using monthly Release-06 GRACE Level-2 Stokes coefficients for the span April 2002 to August 2016, which result in a set of solutions for the time-varying geoid, sea-surface 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 center-of-mass and center-of-figure reference frames. These data may be applied for either study of the fields themselves or as fundamental filter components for the analysis of ocean-circulation- and earthquake-related fields or for improving ocean tide models.

1 Introduction

Geodesists have long understood that the ocean mean sea surface follows the shape of the Earth's geoid (Rapp, 1983) and that changes in on-land water storage are a source of time-varying gravity (Lambert and Beaumont, 1977). The fundamental relationship of changes in land ice and water, solid Earth, and sea-surface 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 supra-seasonal 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 sea-level variability associated with a self-attraction 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 sea-level change: the prediction of “sea-level fingerprints”. Sea-level 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 land-mass-driven 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 on-land hydrology, two-way land-to-ocean 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 ocean-bottom pressure data, in conjunction with space geodetic data, to claim that fingerprints associated with decadal-scale on-land 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 mid-1990s. Galassi and Spada (2017) have noted that the influence of a land-mass-induced 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 sea-level 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 sea-level 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 sea-level fingerprints (e.g., Rietbroek et al., 2012; Frederikse et al., 2018). Consideration of land-ice- and water-driven 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 sea-level 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ía-Gó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 Release-06 GRACE Level-2 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 long-standing 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 GRACE-based estimates of global surface mass transport.

2 Key variables and deliverables

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:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{\Delta }S(\mathit{\theta },\mathit{\varphi },\mathrm{\Delta }t)=\mathrm{\Delta }N(\mathit{\theta },\mathit{\varphi },\mathrm{\Delta }t)-\mathrm{\Delta }U(\mathit{\theta },\mathit{\varphi },\mathrm{\Delta }t),\end{array}$$

where ΔN and ΔU are corresponding changes in sea-surface 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 center-of-mass (CM) of the total Earth system or the center-of-figure (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 sea-surface height, we write

$$\begin{array}{ll}{\displaystyle }\mathrm{\Delta }S(\mathit{\theta },\mathit{\varphi },\mathrm{\Delta }t)=& {\displaystyle \frac{\mathrm{1}}{g}}\mathrm{\Delta }\mathrm{\Phi }(\mathit{\theta },\mathit{\varphi },\mathrm{\Delta }t)-\mathrm{\Delta }U(\mathit{\theta },\mathit{\varphi },\mathrm{\Delta }t)\\ \text{(2)}& {\displaystyle }& {\displaystyle }+\mathrm{\Delta }C\left(\mathrm{\Delta }t\right),\end{array}$$

where ΔΦ is the net perturbation in Earth's surface potential, ΔC is a spatial invariant that explains the discrepancy between the sea-surface 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. Space-based gravity missions, such as GRACE and GRACE Follow-On (GRACE-FO), provide direct measurements of the non-rotational 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 on-land 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 sea-level fingerprint of monthly changes in on-land 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 sea-surface 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.

3 The sea-level equation

Here we briefly summarize the fundamental concept and a numerical technique of solving the so-called “sea-level 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, mass-conserving load function so that

$$\begin{array}{ll}{\displaystyle }L(\mathit{\theta },\mathit{\varphi },t)=& {\displaystyle }H(\mathit{\theta },\mathit{\varphi },t)\left[\mathrm{1}-\mathcal{O}(\mathit{\theta },\mathit{\varphi })\right]\\ \text{(3)}& {\displaystyle }& {\displaystyle }+S(\mathit{\theta },\mathit{\varphi },t)\mathcal{O}(\mathit{\theta },\mathit{\varphi }),\end{array}$$

where $H(\mathit{\theta },\mathit{\varphi },t)$ is the on-land 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 on-land (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 sea-level fingerprint: a unique spatial pattern of relative sea level that is consistent with fundamental physical features of a realistic Earth. For a self-gravitating elastically compressible rotating Earth, we compute the sea-level fingerprint by satisfying the following sea-level equation:

$$\begin{array}{ll}{\displaystyle }S(\mathit{\theta },\mathit{\varphi },t)=& {\displaystyle }E\left(t\right)+{\displaystyle \frac{\mathrm{1}}{g}}\mathrm{\Phi }(\mathit{\theta },\mathit{\varphi },t)-U(\mathit{\theta },\mathit{\varphi },t)\\ \text{(4)}& {\displaystyle }& {\displaystyle }-〈{\displaystyle \frac{\mathrm{1}}{g}}\mathrm{\Phi }(\mathit{\theta },\mathit{\varphi },t)-U(\mathit{\theta },\mathit{\varphi },t)〉.\end{array}$$

The physical interpretation of the right-hand 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 sea-level 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 ocean-averaged value of $(\mathrm{\Phi }/g-U)$. This spatial invariant is essential to ensure that the global mean relative sea-level change is the same as the barystatic term.

To solve for the sea-level 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:

$$\begin{array}{ll}{\displaystyle }S(\mathit{\theta },\mathit{\varphi },t)=& {\displaystyle }X(\mathit{\theta },\mathit{\varphi },t)+Y(\mathit{\theta },\mathit{\varphi },t)+P(\mathit{\theta },\mathit{\varphi },t)\\ \text{(6)}& {\displaystyle }& {\displaystyle }+Q(\mathit{\theta },\mathit{\varphi },t)+C\left(t\right),\end{array}$$

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 }/g-U>$. 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:

$$\begin{array}{ll}{\displaystyle }\widehat{S}(\mathit{\theta },\mathit{\varphi },t)=& {\displaystyle }\widehat{X}(\mathit{\theta },\mathit{\varphi },t)+\widehat{Y}(\mathit{\theta },\mathit{\varphi },t)+\widehat{P}(\mathit{\theta },\mathit{\varphi },t)\\ \text{(7)}& {\displaystyle }& {\displaystyle }+\widehat{Q}(\mathit{\theta },\mathit{\varphi },t)+\widehat{C}\left(t\right),\end{array}$$

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 right-hand 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:

$$\begin{array}{ll}{\displaystyle }{\widehat{S}}_{lm}^r=& {\displaystyle }{\widehat{X}}_{lm}\left({\widehat{S}}_{lm}^{r-\mathrm{1}}\right)+{\widehat{Y}}_{lm}\left({\widehat{S}}_{lm}^{r-\mathrm{1}}\right)+{\widehat{P}}_{lm}\left({\widehat{S}}_{lm}^{r-\mathrm{1}}\right)+{\widehat{Q}}_{lm}\left({\widehat{S}}_{lm}^{r-\mathrm{1}}\right)\\ \text{(8)}& {\displaystyle }& {\displaystyle }+{\widehat{C}}_{lm}\left({\widehat{S}}_{lm}^{r-\mathrm{1}}\right),\end{array}$$

where $r=\mathrm{1},\mathrm{2},\mathrm{\dots }{r}_{\max }$ is the recursion counter, and r_max is the value of r for which the desired convergence is attained. Note that dependence of right-hand 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 non-hatted 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 non-hatted (degree p, order q) coefficients are obtained as well. These non-hatted coefficients automatically satisfy the sea-level equation itself (Eq. 6) in the SH domain; i.e.,

$$\begin{array}{ll}{\displaystyle }{S}_{pq}=& {\displaystyle }{X}_{pq}\left({\widehat{S}}_{lm}\right)+Y_{pq}\left({\widehat{S}}_{lm}\right)+P_{pq}\left({\widehat{S}}_{lm}\right)+Q_{pq}\left({\widehat{S}}_{lm}\right)\\ \text{(9)}& {\displaystyle }& {\displaystyle }+C_{pq}\left({\widehat{S}}_{lm}\right).\end{array}$$

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 self-gravitating Earth are given by X_pq(t) and those for a self-gravitating rotating Earth by [X_pq(t)+Y_pq(t)]. Similarly, SH coefficients of bedrock elevation change are given by −P_pq(t) for a self-gravitating Earth and by $-\left[P_{pq}\left(t\right)+Q_{pq}\left(t\right)\right]$ for a self-gravitating rotating Earth.

4 GRACE and sea-level fingerprints

Here we give a brief summary of the steps undertaken to develop sea-level 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 Level-2 Release-06 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 Release-06 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 post-processed 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 sea-level rise. To do this correctly, we derive degree 1 coefficients from JPL Release-06 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 Release-06 data products (available at ftp://podaac.jpl.nasa.gov/allData/grace/docs/TN-11_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/solid-earth/gia/ (last access: 7 May 2019). Finally, for all coefficients, we remove corresponding 11-year (January 2003–December 2013) mean values to retrieve Strokes coefficient anomalies.

By combining GSM Stokes coefficient anomalies with GIA and low-degree coefficients as noted above, we may derive corresponding coefficients for land water storage anomalies, H_lm(t), as follows (Wahr et al., 1998):

$$\begin{array}{ll}{\displaystyle }{H}_{lm}\left(t\right)=& {\displaystyle \frac{a{\mathit{\rho }}_e}{\mathrm{3}{\mathit{\rho }}_{\mathrm{w}}}}{\displaystyle \frac{\mathrm{2}l+\mathrm{1}}{\mathrm{1}+k_l^{\prime }}}\left\{\exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\ln \left(\mathrm{2}\right)}}{\left({\displaystyle \frac{lr}{a}}\right)}^{\mathrm{2}}\right]\right\}\\ \text{(10)}& {\displaystyle }& {\displaystyle }\left[C_{lm}^{\mathrm{gsm}*}\left(t\right)-C_{lm}^{\mathrm{gia}}\left(t\right)\right],\end{array}$$

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 so-called 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 low-degree Stokes coefficients as noted above.

Figure 1 Barystatic sea-level time series. Our estimates of trends and seasonal amplitudes for all three data centers are compared to JPL mascon solutions (Watkins et al., 2015). Results are plotted relative to the time means over the period January 2003 through December 2013. Trend values are provided for the period January 2005 through December 2015, except for GFZ solutions (January 2005 to December 2013), for a comparison to the sum of individual mass components (during January 2005 to December 2016) listed in Table 13 of the WCRP (2018) report: 1.65±0.23 mm yr⁻¹. As an additional point of comparison, Dieng et al. (2015) find GRACE-determined mass changes for the barystatic sea-level trend at 2.04±0.08 mm yr⁻¹ for January 2005–December 2013 from the mean of CSR, GFZ, and JPL Level-2 Release-05 spherical harmonic solutions.

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 non-overlapping sub-domains 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 non-cryospheric 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 (GLDAS-Noah), 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).

Figure 2 Effects of scaling on the select spherical harmonic coefficients. (a) Scaling effects on the average rate of change in land water storage, relative sea level and geoid height, during April 2002 through March 2016. These solutions are based on JPL Stokes coefficients and are computed in the CM reference frame with the rotational feedback included. (b) Distribution of energy, computed for a given degree l as a sum of the squares of corresponding orders m, for unscaled and scaled solutions. Note the log scale on the x axis of the right panel. Unlike the geoid height change, the relative sea-level change has nonzero energy at l=0 with magnitudes of 1.16 and 2.37 mm² yr⁻² for unscaled and scaled solutions, respectively. Also note that solutions for geoid height change employ a different scale (factor of 4) for appropriate visualization.

Figure 3 Land load function, sea-level fingerprint, and uncertainties therein. Average rate of water equivalent height change in land water storage (a, b) and associated change in relative sea level (c, d) for the period April 2002 to March 2016 are shown with their corresponding 1σ uncertainties. Maps for sea-level fingerprint and uncertainty are produced by assembling the corresponding spherical harmonic coefficients provided with this article. These solutions are based on JPL Stokes coefficients and are computed in the CM reference frame with the rotational feedback included. A zoomed-in map of the Mediterranean Sea (e) is meant to highlight the local variability in sea-level fingerprint. The fingerprint-predicted trends for tide gauges at TRIESTE (1.26±0.18 mm yr⁻¹) and CASCAIS (1.50±0.17 mm yr⁻¹) reflect differences that are comparable to those associated with interdecadal atmospheric pressure trends (0.2±0.2 mm yr⁻¹) and the GIA-related fingerprint (≈0.2 mm yr⁻¹) (see Piecuch et al., 2016, and Stocchi and Spada, 2009, respectively). This illustrates one example of the importance of contemporary sea-level fingerprints for tide gauge data analysis and interpretation. Panel (f) is meant to show that our solutions are comparable to those obtained from a well-validated and higher-resolution ISSM sea-level solver (Adhikari et al., 2016); note that the solution discrepancy is within the 1σ uncertainties (compare with panel d).

A detailed description of scaling may be found in Adhikari and Ivins (2016), who used the same recipe to post-process the CSR Release-05 GRACE Level-2 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 sea-level change (Fig. 3), vertical bedrock motion (not shown), and geoid height change (not shown) are consistent with the respective solutions computed using a well-validated sea-level solver (Adhikari et al., 2016) that operates on an unstructured global mesh of the Ice and Sea-level 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 sea-level fingerprints is generally within 1σ uncertainties (compare Figs. 3d and 4d). We show the origin of discrepancies by plotting the degree variance spectrum.

Figure 4 Comparison of data centers for select fields. JPL solutions are subtracted from CSR and GFZ solutions for trend in land water storage change (a, c) and relative sea-level change (b, d) during January 2003–December 2013. The difference in the spectrum of energy distribution is also shown (e, f). Results are computed in the CM reference frame with the rotational feedback accounted for.

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 sea-level 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 sea-surface height may be obtained by summing coefficients for bedrock motion and those for relative sea-level 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 Release-06 (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 sea-level 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.

Figure 5 Effects of (a) Earth's rotation and (b) reference frame on model solutions. Rotational feedback on the relative sea-level change exhibits a degree 2 spherical harmonic pattern and generally accounts for ≈10 % of the total signal (compare with Fig. 3c). Degree 1 load Love numbers depend on the choice of reference frame origin. The effect of reference frame – quantified here as the model solution computed in the CM frame minus the solution computed in the CF frame – is more pronounced in the vertical bedrock motion and the geoid height change (not shown). These results are based on JPL Stokes coefficients for the period April 2002 to March 2016.

5 Discussion

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 on-land mass changes from ocean altimetry, bottom-pressure, and tide gauge studies. Such removal is essential for future studies of the patterns of sea-level 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 sea-level 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 sea-level 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 1-b 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 sea-surface 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¹⁸ 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 second-order effect.

6 Data availability

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 sea-level fingerprints in both CM and CF reference frames. File names follow the GRACE naming convention. Solutions that correspond to the GSM file “GSM-2_2002095-2002120_GRAC_UTCSR_BA01_0600”, for example, are stored in four files named “SLF-2_2002095-2002120_GRAC_UTCSR_BA01_0600” under appropriate directories; we simply replace “GSM” by “SLF” to denote “sea-level fingerprints”. The time stamp (in YYYYDoY-YYYYDoY format) and the corresponding data center (five-character string containing CSR, GFZ, or JPL) also appear in the file name. The example file considered above contains sea-level 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.

7 Conclusions

In this paper we describe a data product that emerges from the Release-06 GRACE Level-2 Stokes coefficients, provided by CSR, GFZ, and JPL, which contain the basic information necessary to create monthly sea-level 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 real-time two-dimensional imaging of the sea-surface height without the necessity of having to patch together one-dimensional profiles (e.g., Gaultier et al., 2016). In addition to providing higher resolution, this will allow improved accuracy. When coupled to GRACE-FO mapping of gravity changes, we should begin to see the emergence of sea-level 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 self-gravitation and loading. Here we present the effects of land-based 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 Level-2 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 Release-06 GRACE Level-2 products (Wahr et al., 2015; Bettadpur, 2018) have been removed and that the sea-surface 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 land-mass-change-related fingerprints from either GRACE, ocean altimetry, tide gauge, or GPS-determined vertical land motion data from April 2002 to August 2016 using the monthly solutions we supply here as a data product.

Appendix A: Theory of the sea-level fingerprint

The fundamental theoretical concept of the so-called sea-level 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\}\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 })\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)\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 great-circle distance, α, as follows: $\cos \mathit{\alpha }=\cos \mathit{\theta }\cos {\mathit{\theta }}^{\prime }+\sin \mathit{\theta }\sin {\mathit{\theta }}^{\prime }\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(\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 }+\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₂ and h₂ 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 leading-order 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 }\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 }\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, σ₀ is the so-called 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₁₁ and ΔI₂₂ (Eq. A6) and ΔI₃₃ (Eq. A7). Similarly, the terms enclosed by outer parentheses represent Earth's rotation parameters: polar motion (m₁,m₂) (Eq. A6) and change in length of day m₃ (Eq. A7).

• The ocean-averaged term in Eq. (4), denoted by 〈*〉, may be written as follows:

  $$\begin{array}{ll}{\displaystyle }& {\displaystyle }〈{\displaystyle \frac{\mathrm{1}}{g}}\mathrm{\Phi }(\mathit{\theta },\mathit{\varphi },t)-U(\mathit{\theta },\mathit{\varphi },t)〉={\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 }-U(\mathit{\theta },\mathit{\varphi },t)\mathit{\left\}}\mathcal{O}\right(\mathit{\theta },\mathit{\varphi }\left)\mathrm{d}\mathbb{S}\right].\end{array}$$

  When rotational feedback is excluded, Φ and U should be replaced by Φ_g and U_g, respectively.

Appendix B: Spectral methods for sea-level equation

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 square-integrable 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\left|m\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{(l-|m\left|\right)\mathrm{!}}{(l+|m\left|\right)\mathrm{!}}}}\\ \text{(B2)}& {\displaystyle }& {\displaystyle }{\mathcal{P}}_{l\left|m\right|}\left(\cos \mathit{\theta }\right)\left\{\begin{array}{l}\cos \left(m\mathit{\varphi }\right)\mathrm{if}m\ge \mathrm{0}\\ \sin \left(\right|m\left|\mathit{\varphi }\right)\mathrm{if}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}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(\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 great-circle 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)\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)\mathrm{d}\mathbb{S}.\end{array}$$

• Evaluation of surface integrals on a unit sphere. We discretize the surface of a unit sphere using the so-called 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)\mathrm{d}\mathbb{S}={\displaystyle \frac{\mathrm{4}\mathit{\pi }}{N_{\mathrm{T}}}}\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}\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)\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 {\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)\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 {\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 𝒴₀₀(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)\mathrm{d}\mathbb{S}=\mathrm{4}\mathit{\pi }{\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 {\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 {\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]\mathrm{d}\mathbb{S}=\mathrm{4}\mathit{\pi }{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 sea-level 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)\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)\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)\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 {\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(\cos \mathit{\alpha }\right)L\left({\mathit{\omega }}^{\prime }\right)\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)\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)\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)\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)\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-〈\mathrm{\Phi }/g-U〉\equiv E-〈X+P+Y+Q〉$. 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[{\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)\right.\\ {\displaystyle }& {\displaystyle }\left.+{\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\}\right]\mathrm{d}\mathbb{S}.\end{array}$$

  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 right-hand 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[{\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.\\ {\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)\mathrm{d}\mathbb{S}.\end{array}$$

  Since 𝒴₀₀(ω)=1, we introduce a virtual expression 𝒴₀₀(ω) 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.+\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\}\right].\end{array}$$

  The terms inside the braces vanish when rotational feedback is not included. Note that the first term on the right-hand 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 on-land 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 non-hatted counterparts as in, for example, Eq. (B18). For the chosen solid Earth model and the reference frame origin, we compute these non-hatted 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₂ and h₂. 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 L-2 norm after each recursion and call the solution converged when it is less than 0.001 % of the L-2 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 non-hatted 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.