Home Greenland Antarctica Arctic Methods Code Publications/Press About

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 \begin{equation} 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. \end{equation} The Slepian function coefficients, $g_{\alpha,lm}$, are found by solving the eigenvalue equation \begin{equation} \sum_{l'=0}^{L} \sum_{m'=-l'}^{l'} D_{lm,l'm'} g_{l'm'} = \lambda g_{lm}. \end{equation} 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 \begin{equation} \int_R Y_{lm}Y_{l'm'} \,d\Omega = D_{lm,l'm'}. \end{equation} 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: \begin{equation} \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}. \end{equation} 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: \begin{equation} N = (L+1)^2 \frac{A}{4\pi}, \end{equation} The truncated sum lets us represent the data very sparsely, yet with very good reconstruction properties within the region [Simons, 2010]: \begin{equation} {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 . \end{equation} 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. 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. 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. 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.