The methods we use in this analysis are described in the Supporting Information from our recent paper Harig and Simons [2012]. Here we present a summary of the technique and outline the general principle of Slepian localization.

In geophysics, it is quite common to be interested in some continuous data (such as gravity) in some specific region of interest (such as Greenland). A common way to represent and analyze this data is by using spherical harmonics. Since spherical harmonics spread their energy over the entire globe, they are not usually amenable to regional use without some sort of numerical regularization in the signal inversion. This is due to the fact that spherical harmonics are orthogonal over the globe, but not the region.

Instead, we desire an orthogonal set of functions on the sphere that is both optimally concentrated in our spatial region of interest and bandlimited to a chosen degree. These functions will have most of their energy within our region, and will have minimal energy outside. For this purpose we use spherical Slepian functions, which are analogs to the classic Slepian concentration problem [Slepian, 1983; Wieczorek and Simons, 2005; Simons et al., 2006; Simons and Dahlen, 2006]. The new functions are defined in a similar way to spherical harmonics, as $$g_\alpha(\mathbf{\hat{r}}) = \sum_{l=0}^L \sum_{m=-l}^l g_{\alpha , lm} Y_{lm}(\mathbf{\hat{r}}), \quad\quad g_{\alpha,lm} = \int_\Omega g_\alpha(\mathbf{\hat{r}}) Y_{lm}(\mathbf{\hat{r}}) \,d \Omega.$$ The Slepian function coefficients, $g_{\alpha,lm}$, are found by solving the eigenvalue equation $$\sum_{l'=0}^{L} \sum_{m'=-l'}^{l'} D_{lm,l'm'} g_{l'm'} = \lambda g_{lm}.$$ If we were considering the entire sphere, then the integral of the products of spherical harmonics would be the delta function, illustrating that they are orthogonal. When we just consider the region $R$, the integrals are no longer delta functions, and instead are $$\int_R Y_{lm}Y_{l'm'} \,d\Omega = D_{lm,l'm'}.$$ Here $\mathbf{D}$ is a matrix whose elements $D_{lm,l'm'}$ are products of spherical harmonics integrated over the region $R$. By performing an eigenvalue decomposition of this matrix, we find the functions that have the most energy concentrated in our region: $$\lambda = \frac{\displaystyle \int_R g_\alpha^2(\mathbf{\hat{r}})\,d\Omega} {\displaystyle \int_\Omega g_\alpha^2(\mathbf{\hat{r}})\,d\Omega} = \mbox{maximum}.$$ We can then use a truncated sum of just these most concentrated functions as a basis for our region. Commonly we truncate up to the Shannon number, $N$, which is a function both of the region and of the bandwidth of our data: $$N = (L+1)^2 \frac{A}{4\pi},$$ The truncated sum lets us represent the data very sparsely, yet with very good reconstruction properties within the region [Simons, 2010]: $${d}(\hat{\mathbf{r}})\approx \sum_{\alpha=1}^{N} d_{\alpha}\,g_{\alpha}(\hat{\mathbf{r}}) \quad\mbox{for}\quad \hat{\mathbf{r}}\in R .$$ In this manner, the process is analogous to a singular-value decomposition in the pixel basis depending on inverse eigenvalues of $D$, truncated when they get large. These properties make the Slepian basis is an ideal tool to conduct estimation problems that are linear or quadratic in the data [Simons and Dahlen, 2006; Dahlen and Simons, 2008].

Below we show an example of a set of Slepian functions that are concentrated in a region around Greenland, for a certain spherical harmonic bandwidth up to degree and order 60.
Figure 1: Slepian eigenfunctions $g_1, g_2, ..., g_{12}$ that are optimally concentrated within a region buffering Greenland by 0.5$^\circ$. Dashed line indicates the region of concentration. Functions are bandlimited to $L = 60$ and are scaled to unit magnitude. The parameter$\alpha$ denotes which eigenfunction is shown. The parameter $\lambda$ is the corresponding eigenvalue for each function, indicating the amount of concentration. Magnitude values whose absolute values are smaller than 0.01 are left white.
The GRACE level-2 data are released as spherical harmonic coefficients up to degree and order 60. This corresponds to $(L + 1)^2 = 3721$ coefficients. We transform the data into a Slepian basis for Greenland, and after truncation are instead left with just 20 functions and their coefficients. This transformation increases the signal-to-noise properties and allows us to fit functions to the coefficients as they change through time. Below we show what some of these fits look like for some of the functions.
Figure 2: Time series of various ($\alpha$ = 1, 5, 6, 7, 11, 20) Slepian coefficients and their best-fit polynomial (blue lines). Each coefficient is fit by an annual periodic and linear function, as well as quadratic and cubic polynomial terms if those terms pass an $F$-test for variance reduction. Shown here are the coefficient and fitted function values with the annual periodic function subtracted from both. The mean is removed from each time series.
Below we show what these fits look like in space. Each Slepian function is expanded in space and multiplied by its fitted coefficient averaged for a year. We can see that the magnitude of mass loss is mostly represented by the first five functions. The remaining functions, while they have a low mass integral change, are important for representing the signal in space and showing the full pattern of mass change.
Figure 3: Predicted GRACE annual mass change in the Slepian basis for each of the first twelve eigenfunctions. Each eigenfunction, denoted by its index $\alpha$, is scaled by the total change in that coefficient from 1/2003-11/2011 divided by the time span (yrs), expressed as the cm/yr water equivalent of surface density. Thus, this represents the mass change for an average year during this time span. The inset variable "Int" displays the integral of each function in the concentration region within the dashed line expressed as the mass change rate of gigatons per year. Surface-density change of absolute value smaller than 0.1 cm/year is left white.
Finally, we can examine the time series for the three most-contributing Slepian functions (below), which we express as the integral of the product of the expansion coefficient and the function. It is easy to see that the mass signal trends can be well estimated relative to the variance seen from month to month. The Slepian functions significantly enhance signal-to-noise within the region of interest compared to traditional spherical harmonics.
Figure 4: Mass change in gigatons (Gt) for the three most significant Slepian-function terms ($\alpha$ = 1, 3, 11), which contribute more than 70% of total mass change over the data time span. Monthly data are drawn as the solid black lines while the $2\sigma$ uncertainty envelopes are drawn in grey. Each function has a mean of zero but has been offset from zero for clarity.

To aid other researchers and follow the tenets of reproducible research, we have posted online the Matlab code we used to perform the GRACE research presented in this website. Currently the code is hosted on Frederik's software page http://geoweb.princeton.edu/people/simons/software.html#grace. Below we've listed the code available with descriptions, and outlined a general set of instructions by which to use it. Since it depends heavily on other code Frederik has written and posted, we recommend downloading all of the functions on the software site, and adding that directory to your matlab path. Some of these functions also require the Matlab mapping toolbox.

Making the most of GRACE
Function Description
Clmlmp2Crrp Given a spectral covariance matrix, evaluates it in space
Clmlmp2Cab Given a spectral covariance matrix, turns it into a Slepian covariance matrix
cov2plm Given a spectral covariance matrix, generates spherical harmonics realizations
grace2plmt Turns monthly GRACE data files into a single matrix for time-dependent analysis. This function does corrections for C2,0 from the SLR values of Cheng and Tapley, [2004] and degree 1 values from Swenson, [2008].
grace2slept Transform the result of GRACE2PLMT into a Slepian basis
grs Computes parameters for a certain geodetic reference system
integratebasis Integrates Slepian eigenfunctions given as spherical harmonics expansions
kernelcp Finds the kernel whose eigenfunctions are optimally concentrated (parallel version of KERNELC)
lovenums Returns elastic Love numbers for a certain Earth model
periodfit Find and fit periodic cycles through a data set
plm2avg Integrates and averages spherical harmonic expansions
plm2avgp Integrates and averages spherical harmonic expansions (parallel version of PLM2AVG)
pmlt2diff Turns monthly GRACE data matrix into a month-to-month difference map
plmt2resid Turns monthly GRACE data matrix into residuals after fitting analysis in the spherical harmonic basis
plmresid2cov Turns GRACE residual time series into a spherical-harmonic spectral covariance matrix
slepresid2cov Turns GRACE residual time series into a Slepian covariance matrix
resid2plot Plots GRACE residual time series
rotateGp Rotates a matrix returned by GLMALPHA (parallel)
slept2resid Turns monthly GRACE data matrix into residuals after fitting in the Slepian basis
timeseriesfit Fits polynomial functions to time series with an F-test criterion

For the specific results here on the website and in our paper Harig and Simons [2010], we have posted the functions that make the figures on Frederik's software page. If you want to study other regions or signals, here is the general outline with which to do that.

1. Set up your data and directory structure, and then with GRACE2PLMT read in the GRACE data files into a matrix for use in Matlab. This function does corrections for C2,0 from the SLR values of Cheng and Tapley, [2004] and degree 1 values from Swenson, [2008].

2. Decide on your choice of basis, depending on your region of interest and the bandwidth you want. Using GRACE2SLEPT, project the results of GRACE2PLMT into your chosen basis. We recommend that your basis is chosen based on a set of synthetic experiments which estimate the leakage/recovery tradeoffs.

3. Next run SLEPT2RESID to fit a choice of functions (e.g. lines, quadratics, etc.) to the Slepian coefficients. Note: if you want to remove a model of GIA then you should do that before this step, using a function like CORRECT4PGR.

4. If you want to examine the total mass trend in a region, this information is calculated and returned as a result from SLEPT2RESID. To summarize, each Slepian coefficient up to the Shannon truncation is multiplied by the integral (found using INTEGRATEBASIS) of the corresponding function over the region. This total data is then fit with TIMESERIESFIT which allows fitting with a custom data variance.

5. If you want the total map of mass change, multiply the difference between estimated signal coefficients by the corresponding Slepian function, and add them up. Remember the appropriate units. This sum can then be expanded to space using PLM2XYZ By limiting the set of coefficients you use, you can instead make this for any date span, such as for several years or a single year.

References

Harig, C. and F. J. Simons. Mapping Greeenland's mass loss in space and time. In Press, PNAS. 2012.

Dahlen, F. A & Simons, F. J. (2008) Spectral estimation on a sphere in geophysics and cosmology. Geophys. J. Int. 174, 774-807.

Slepian, D. (1983) Some comments on Fourier analysis, uncertainty, and modeling. SIAM Review 25, 379-393.

Simons, F. J, Dahlen, F. A, & Wieczorek, M. A. (2006) Spatiospectral concentration on a sphere. SIAM Review 48, 504-536.

Simons, F. J & Dahlen, F. A. (2006) Spherical Slepian functions and the polar gap in geodesy. Geophys. J. Int. 166, 1039-1061.

Simons, F. J. (2010) in Handbook of Geomathematics, eds. Freeden, W, Nashed, M. Z, & Sonar, T. (Springer-Verlag), pp. 891-923.

Wieczorek, M. A & Simons, F. J. (2005) Localized spectral analysis on the sphere. Geo136 phys. J. Int. 162, 655-675.