## Introduction

Understanding of the global water cycle, the exchange of water mass among oceans, atmosphere and land, is essential to decipher Earth's climate change. The global water and ice mass redistribution occurs via precipitation, evapotranspiration, sublimation, surface/subsurface water flow and ice flow. Continuous hydrologic, meteorologic and cryospheric observations for theses phenomena are important to understand and predict climate changes and to manage corresponding natural hazards (IPCC, 2014). Despite much efforts for the observation, sparse in-situ stations hinder understanding the water and ice mass balance. This limitation is particularly true for the Antarctica and Greenland.

The Earth’s gravity monitoring from space can be an alternative observation of water cycle because the redistribution of water mass can cause changes in the gravity field (Cazenave and Chen, 2010). In particular, the Gravity Recovery and Climate Experiment (GRACE) (2002-2017) mission has provided unprecedented observations for changes in the gravity field. The GRACE mission consists of co-orbiting twin satellites apart from about 220 km at altitude of ~450 km (Tapley et al., 2004). Onboard K-band microwave interferometers measure inter-satellite distance at accuracy of ~μm (Dunn et al., 2003). The measurements are used to determine the gravity field and further estimate water mass variations with accuracy of less than 1 cm equivalent water height. This is about 2 orders of magnitude better than those of EGM96, the pre-GRACE gravity model (Wahr et al., 2004). Since its launch, the monthly gravity anomalies provided by GRACE have been widely used in hydrologic (Alsdorf et al., 2010;Eom et al., 2017a;Reager et al., 2014), atmospheric (Kim et al., 2016;Seo et al., 2014), oceanographic (Leuliette and Miller, 2009;Llovel et al., 2014;Makowski et al., 2015) and cryospheric (Eom et al., 2017b;Seo et al., 2015) studies with global to regional scale.

However, there are some limitations of the monthly GRACE solutions for a region with small spatial scale less than a few hundred kilometers (Longuevergne et al., 2010). Monthly gravity solutions are determined by integrated along-track observations in combination with a background gravity model during the observation epoch. During the monthly sampling, sub-monthly gravity fluctuations cause aliasing error including spurious longer period variabilities and peculiar meridional stripes (Wahr et al., 2004). To suppress the aliasing error, numerical models are incorporated into GRACE’s data processing in order to remove sub-monthly fluctuation of atmospheric surface pressure and ocean bottom pressure (Bettadpur, 2012). However, those models are insufficient to reduce such aliasing error: high-frequency gravity residuals remain in GRACE observation resulting in significant aliasing error.

To remedy the aliasing error, spatial filters are commonly applied. After filtering, the error and noise can be remarkably reduced, leading to reveal the pattern of water mass signal but showing broader spatial features with attenuation. To overcome such shortcomings, several approaches have been developed. For example, Velicogna et al. (2005) used an empirical scale factor to compensate the signal loss to nearby oceans when estimating Greenland ice mass loss. More recently, Chen et al. (2013) used a forward modeling method that iteratively adjusts mass fields until its smoothed mass fields are close to the GRACE’s observation after spatial filtering. These signal recovering methods are efficient to estimate basin-averaged mass changes, but limited to understand details of spatial pattern because the former method does not alter the signal pattern and the latter suffers from spurious signal patterns.

The empirical orthogonal functions (EOF) method can be an alternative for GRACE data to separate water mass signal from aliasing errors. The method decomposes a space- and time-variable dataset into spatial patterns and associated time coefficients describing each spatial pattern’s temporal variability (Hannachi et al., 2007). In atmospheric science, therefore, EOF analysis has been frequently used to find some dominant climate modes such as AO and MJO in climate data (Hurrell and Van Loon, 1997;Madden and Julian, 1971). Based on this feature, EOF method can be also used as a filter by exclusions of non-significant modes. In this study, we suggest the way to use EOF method as a filter for GRACE monthly solutions. We examine the EOF method in the Western Antarctic Ice Sheet (WAIS) and the Antarctic Peninsula (AP) where the conventional GRACE filtering is particularly problematic to recover accurate ice mass loss pattern.

## Data and Methods

### GRACE data

The monthly GRACE solutions are freely available from the ftp server at the NASA PODAAC (ftp:// podaac.jpl.nasa.gov). Conventionally, a global gravity model is represented by spherical harmonic (SH) function because the gravitational potential field outside Earth’s surface satisfies Laplace’s equation. Therefore, the geoid height observed by GRACE is also described by coefficients of SH function, known as Stoke’s coefficients:(1)

where Δ*N* is geoid height variations at colatitude, *θ*, and longitude, *ϕ* at a given month *t*. *R _{e}* is the Earth's radius and ${\tilde{P}}_{lm}$ is normalized associated Legendre functions. The subscripts

*l*and

*m*represent the degree and order of SH coefficients, respectively. Δ

*C*and Δ

_{lmt}*S*are Stokes coefficients after removing mean values representing static geoid. The solutions used in this study are provided by the Center for Space Research (CSR) at University of Texas at Austin, one of the three official GRACE data processing centers. Recently, CSR released a new version of gravity solutions, RL06. A solution consists of SH coefficients up to degree and order 60. In this study, 149 monthly solutions from January 2003 to August 2016 are used and linearly interpolated 15 omitting months completing full-time series for more than 13 years.

_{lmt}There are several procedures to derive surface mass distribution from GRACE geoid height data. Table 1 outlines the post-processing flow-chart for GRACE SH solutions, which are commonly adopted by GRACE studies (Velicogna and Wahr, 2013). First, there are four geophysical data reductions as follows. (i) GRACE’s range-rate observation using on-board interferometers is not sensitive to the center of mass changes, which are represented by degree 1 SH coefficients, ΔC_{10}, ΔC_{11} and ΔS_{11}. Commonly these three coefficients are synthesized with the other highorder SH coefficients assuming that ocean mass changes are passively determined by terrestrial water and ice mass variations (Swenson et al., 2008). It is necessary to add the estimated monthly values of the degree 1 coefficients to GRACE SH solutions. (ii) GRACE ΔC_{20} shows significant aliasing errors due to S_{2} and K_{2} tides (Chen and Wilson, 2008), and are replaced with Satellite Laser Ranging (SLR) observations (Cheng and Tapley, 2004). (iii) The effect of viscoelastic mantle responses to the last glaciation, glacial isostatic adjustment (GIA), is corrected by the model based on ICE-5G ice history (*A* et al., 2013). (iv) Time-varying Stokes coefficients (Δ*C _{lmt}* and Δ

*S*) representing geoid changes are converted to SH coefficients for surface mass density ($\Delta {\tilde{C}}_{lmt}$ and $\Delta {\tilde{S}}_{lmt}$) considering elastic loading response of solid Earth (Wahr et al., 1998).

_{lmt}Second, because the reduced SH coefficients for surface mass density suffer from the aliasing error and random noise, two successive filtering procedures are usually practiced. Swenson and Wahr (2006) found that the aliasing error evidently exhibits correlated patterns and hence developed a de-correlation filter to suppress the error. In addition, random noise is suppressed by applying a kernel function for spatial average (Jekeli, 1981). The shape of the spatial kernel is a Gaussian distribution, and the kernel’s radius, *r*, was defined at which the kernel’s amplitude drops to a half. The radius range of about 300-500 km is empirically proven to effectively reduce most residual noise after de-correlation filter (Velicogna and Wahr, 2013).

Therefore, there are three different reduced SHs for surface mass density; ones without any spatial filtering, the other with the de-correlation filtering only and another with the de-correlation filtering followed by 300-km Gaussian filtering (Table 1). Finally, these three SHs representing global surface density fields are localized on our study region including WAIS and AP (Harig and Simons, 2012). We named these datasets ‘no-filtered’, ‘de-correlated’ and ‘two-step’ (represented by gray-, blue- and greencolored boxes in Table 1, respectively) for simplicity.

### 2.2. Empirical Orthogonal Functions

EOF is an appropriate method for multivariate geophysical data field with space-time domain. The method is also known as principal component analysis or loading factor analysis. Although all methods have different nomenclatures and approaches of data interpretation, they all share a common mathematical technique, singular value decomposition (SVD). SVD decomposes a given dataset into orthogonal matrices; it is useful for EOF to separate signal and errors from GRACE observations.

The monthly GRACE solutions can be arranged into a data matrix with size of *n*×*p*. For our GRACE data, *n* is the number of months (164 in this study), and *p* is the number of SH coefficients equals to 3721 when the maximum degree is 60. For a given month (*t*), SH mass anomaly coefficients ($\Delta {\tilde{C}}_{lmt}$ and $\Delta {\tilde{S}}_{lmt}$) can be rearranged into a column vector:

To construct a data matrix, all column vectors for the study period are merged into an *n*×*p* matrix(3):

where ** D** is a data matrix consisting of monthly solutions according to Eq. (2).

A data matrix is normally weighted by the number of coefficients at each degree prior to SVD. This procedure is equivalent to the latitude weighting for equal-spaced grid data, which is non-uniformly distributed on the Earth’s surface (i.e., denser poleward). Similarly, the degree weighting for the data matrix consisting of SH coefficients attains uniform contribution of each coefficient into the variance of a decomposed mode. Since the number of coefficients for the degree *l* is 2*l*+1, the diagonal weighting matrix ** W** is(4)

Then the weighted data matrix is

In the EOF analysis, spatial patterns and associated time series of dataset are decomposed by SVD technique:

where both ** U** and

**are orthogonal matrices. Each column vector in**

*V***and**

*V***represents spatial pattern and associated time variability, respectively.**

*U***is a**

*S**n*×

*p*rectangular matrix whose diagonal elements represent singular values of

**, which equal the explained variance of each mode (spatial pattern and associated temporal variation). There are many highlevel computing languages to perform SVD calculation such as**

*D**Fortran*,

*Python*and MATLAB

^{®}. For example, we can easily obtain

**,**

*U***and**

*S***by typing ‘[U,S,V] =svd(D)’ in MATLAB**

*V*^{®}command line. The simplicity of SVD computation makes the EOF filter a powerful filtering strategy for GRACE SH solutions.

## Results

### Long-term variations in GRACE data

Firstly, we examine long-term variations in ice mass change using the three GRACE datasets (‘no-filtered’, ‘de-correlated’ and ‘two-step’). Since the correlated striping errors are random in time, they do not significantly affect the trends over regions where mass changes steadily increase or decrease. This implies that in our study region, WAIS and AP, the no-filtered dataset can provide accurate long-term variations. On the other hand, the long-term variations from the filtered date sets such as ‘de-correlated’ and ‘twostep’, are likely underestimated due to spatial filtering. This examination highlights signal smoothing effects when filters are applied on GRACE data.

There are many previous studies reporting that WAIS and AP have experienced rapid decreases in ice mass with acceleration in recent decades (IMBIE, 2018). Our three datasets also show similar results, but there are minor differences depending on filtering strategies. Figure 1 shows the annual rate of ice loss (a-c) and its acceleration (d-f) in each GRACE dataset. In Fig. 1a, large ice loss is observed in the Southern AP and the Amundsen sector of WAIS. A positive anomaly may represent ice mass gain in the center of WAIS, which is likely associated with the stagnation of ice flow in Kamb ice stream (Retzlaff and Bentley, 1993). The acceleration term in the nofiltered dataset (Fig. 1d) shows similar pattern with Fig. 1a in the Amundsen sector, but AP have a dipole structure with ice loss deceleration in the Southern region and acceleration in the Northern region.

Figures 1b and 1e show results of de-correlation filter on long-term variations in GRACE estimates, which are close to those of ‘no-filter’ cases (Figs. 1a and 1d). However, there are some differences in spatial pattern and its amplitude of ice mass changes. The positive/negative structure in WA region is not clearly separated, compared to those of no-filtered dataset. In AP region, both of trend and acceleration are slightly reduced compared to Figs. 1a and 1d, respectively. These differences in spatial pattern and signal intensity arise because de-correlation filtering removes north-south stripes without distinguishing between signals and aliasing errors. Therefore, the signals with spatial pattern elongated in north-south direction are likely suppressed during de-correlation filtering.

Figures 1c and 1f show the results when a Gaussian averaging filter is additionally applied to the GRACE data after de-correlation filtering. Amplitudes of the estimated trend and acceleration (Fig. 1c and 1f) are much more reduced than those shown in the previous results. The maximum values of ice loss rate and its acceleration in this case are 201.4 Gton/year and -0.5 Gton/year^{2}, respectively. The recovered estimates by the two-step filter are about 50 % smaller than those of no-filtered dataset in WAIS. This reduction agrees well with numerical results of *Velicogna* and *Wahr* (2013). Through this analysis, we conclude that additional caution is necessary when GRACE data is reduced by both de-correlation and Gaussian filtering.

### EOF decomposition and reconstruction

Before applying EOF filter, we remove trend and acceleration components from the no-filtered dataset. As described above, the long-term variations are likely recovered from the no-filtered dataset with high accuracy. Indeed, the ice loss rate and its acceleration shown in Figs. 1a and 1d generally agreed with the estimates from recent study (IMBIE, 2018). During GRACE era, there has been a consensus about the long-term ice mass variations over the Antarctic continent (IMBIE, 2018). However, there are still large uncertainties about seasonal ice mass variations with a sub-basin scale (Velicogna et al., 2014;Xiaolei et al., 2014). To examine this, we apply EOF filtering method to the de-trended dataset.

As explained in Section 2, the EOF filtering is divided into two procedures: SVD decomposition of a given data matrix and data reconstruction with selected modes including signals. We therefore decompose the de-trended and no-filtered dataset using SVD. Figure 2 shows time series of the leading four modes (the *k*th column of ** U** in Eq. (6), hereafter referred as to PC

*k*). The PC

*1*and PC

*2*represent apparent inter-annual variations, which are likely associated with ice mass variations. In contrast, PC

*3*and PC

*4*show random variations resulting from observational noise and aliasing error. Higher modes also show the similar white noise (not presented here). This indicates that only the leading two modes from SVD decomposition of the dataset represent ice mass signals.

A similar conclusion also be led by visual inspections of spatial modes. Figure 3 shows results for the first 4 EOF spatial modes (hereafter the *k*th column of ** V** referred as to EOF

*k*). Before presenting spatial modes,

**is restored by the effect of degree weighting in Eq. (5) (i.e.,**

*V*

*W*^{−1}

**). Figure 3a (EOF**

*V**1*) reveals a spatially coherent structure of ice mass loss. The first mode explains 89% of the total variance. In addition, EOF

*1*displays the similar spatial patterns shown in Fig. 1a and 1d. This indicates that a large fraction of ice mass fluctuations arises in a region where ice loss and its acceleration are dominant. EOF

*2*shows a dipole pattern separating between AP and WAIS, implying that different phases of ice mass variations in AP and WAIS. Although the second mode accounts only for 4.5 % of total variance, larger than three times of the explained variance of 3rd mode, it plays an important role in spatial and temporal variability of ice mass change in the study area.

The EOF patterns in higher modes shows characteristics of the aliasing errors (Figs. 3c and 3d). In contrast to the first two EOFs, these represent distinctive patterns of meridional stripes in the association with noise and aliasing error. It is reasonable to conclude that these noisy modes should be excluded before examining the ice mass changes in the area. Note that third and higher modes account for 5.6 % of the total variance, which is larger than that of the second mode. This implies that amplitudes of some important ice mass variability in the area are smaller than those of aliasing errors so that the signal and error are hardly separated from each other using conventional filters.

Based on the inspections of temporal and spatial modes, we conclude that the first two modes include most of ice mass signals, and the signal can be reconstructed with only the two modes. The reconstruction is achieved by multiplying the three matrices with the third and higher columns to be zero in the second term of Eq. (6). This simple algebra produces the EOF-filtered dataset (red box in Table 1).

### Month-to-month mass changes in GRACE data

We calculate monthly difference of ice mass changes using the four GRACE datasets (‘no-filtered’, ‘de-correlated’, ‘two-step’ and ‘EOF reconstructed’). These represent month-to-month ice mass variability at sub-basin scale except the contribution of the trend and acceleration components. Figure 4a shows the monthly differences in ice mass from the no-filtered dataset. The map is dominated by distinctive northsouth stripes over whole study region, which are similar with EOF3 and EOF4 shown in Fig. 3. Note that the scale bars in Fig. 4 differ for each plot, which indicates that aliasing errors are significantly amplified when ice mass changes are subtracted between successive months. It further implies, during the considered time period, signal variations are relatively small compared to the aliasing error and noise.

Figure 4b shows the result from the de-correlated dataset. Most of stripes disappear due to the effect of de-correlation filter, but pattern of ice mass changes much differs from what we expected in the region. Small scale positive and negative anomalies are likely artifacts. Figure 4c is a smoothed version of Fig. 4b with the Gaussian spatial averaging. Although the amplitude is significantly reduced, the spatial pattern of ice mass change is still similar with Fig. 4b. Finally, Figs. 4d shows month-to-month variations from the EOF-filtered dataset. The spatial pattern is very similar to those of Figs. 1a and 1d, which represents regions arising ice mass changes. If we presume that regions experinecing steady ice loss for decades also exhibit fluctuations in seasonal ice mass changes, we could conclude that the similarity in spatial pattern between Figs. 1a and 4d provides that EOF filter is most effective in suppressing aliasing errors without signal loss.

## Summary

GRACE SH solutions are contaminated by aliasing errors due to sub-monthly atmospheric and oceanic fluctuations, which cannot be removed by numerical models during pre-processing procedures. A conventional filtering strategy to reduce the errors is a two-step procedure consisting of decorrelation and Gaussian filters. Although it is effective, the conventional method may also cause signal loss and hinder examination over a sub-basin scale. The limitation is particularly true in WAIS and AP because ice mass signal has a similar pattern of aliasing errors.

In this study, we suggest the EOF filtering method to overcome the limitations of conventional filtering methods. The feature of EOF analysis decomposing a given dataset into two orthogonal matrices can be applicable to geophysical multi-variable data such as GRACE gravity solutions. Using SVD decomposition, signals and noise are separated into different modes, and we simply reconstruct a filtered-data with signal dominant modes. The comparison with conventional filters over the WAIS and AP regions shows the excellence of EOF filter to reduce aliasing errors.

In May 22, 2018, the GRACE Follow-On (FO), a successor to GRACE, was successfully launched at Vandenberg, California. GRACE FO is the first satellite-to-satellite gravity mission equipped with the laser ranging interferometer instrument, and the accuracy is expected to be increased by about ten times compared with the GRACE mission. The gravity solutions provided by GRACE FO are expected to continue its wide applications for hydrologic, oceanographic, atmospheric and cryospheric studies. Therefore, EOF analysis on the SH solutions will also be useful for GRACE FO.