A novel spatial filter to reduce north–south striping noise in GRACE spherical harmonic coefficients

Prevalent north–south striping (NSS) noise in the spherical harmonic coefficient products of the satellite missions gravity recovery and climate experiment greatly impedes the interpretation of signals. The overwhelming NSS noise always leads to excessive smoothing of the data, allowing a large room for improvement in the spatial resolution if this particular NSS noise can be mitigated beforehand. Here, we put forward a new spatial filter that can effectively remove NSS noise while remaining orthogonal to physical signals. This new approach overcomes the limitations of the previous method proposed by Swenson and Wahr (2006), where signal distortion was large and high-order coefficients were uncorrectable. The filter is based on autocorrelation in the longitude direction and cross-correlation in the latitude direction. The NSS-type noise identified by our method is mainly located in coefficients of spherical harmonic order larger than about 20 and degree beyond 30, spatially between latitudes ± 60°. After removing the dominating NSS noise with our method, a weaker filter than before is added to handle the residual noise. Thereby, the spatial resolution can be increased and the amplitude damping can be reduced. Our method can coincidentally reduce outliers in time series without significant trend bias, which underpins its effectiveness and reliability.


Introduction
GRACE's global gravity field products in the form of SHCs are known to suffer from overwhelming NSS noise (Fig. 1a), which is supposed to originate from its polar orbit, lack of cross-track sensitivity, and aliasing due to sampling frequency. The NSS noise in the raw observations is so strong that the physical signals are nearly indistinguishable, so how to remove/suppress the NSS noise is the priority in GRACE data post-processing. A number of miscellaneous filters have been developed for this purpose. Here, we provide a new spatial filter that makes the most of the geometric pattern of NSS noise.
Each element of SHCs is assumed to contain three components: signal, NSS noise (also termed correlated noise), and random noise. The idea behind GRACE data processing is to preserve the signal part and suppress the other noise parts. Generally, there are two strategies to achieve this goal: (i) applying the filter to the coefficient directly and (ii) proposing a pattern to disaggregate the NSS noise component and impose an additional filter to suppress residual noise.
The implementation of the first strategy can mainly be separated into spatial and temporal filters. The main difference between temporal and spatial filters is whether the filter can be applied to a single month of observations. Examples of spatial filters include the isotropic Gaussian filter (Wahr et al. 1998) and its non-isotropic variants (Han et al. 2005;Zhang et al. 2009), the Wiener filter (Klees et al. 2008;Sasgen et al. 2006), the DDK filter (Kusche 2007;Kusche et al. 2009), and the regularization filter (Devaraju and Sneeuw, 2017). These filters differ in the structure (e.g., whether they are homogeneous or inhomogeneous, isotropic or non-isotropic and for more details see (Devaraju and Sneeuw, 2017), whether the signal and noise covariance matrices are utilized, and how to derive prior covariance information). Temporal filters include the EOF filter Wouters and Schrama 2007), the stochastic filter (Wang et al. 2016), the SSA filter (Prevost et al. 2019), and the least square filter (Crowley and Huang 2020). The temporal filters regard the NSS noise as ) × DDK2 indicates that the difference between the results of DDK6 and DDK3 is further smoothed by DDK2 to emphasize the low-degree contribution. A lower index following the DDK filter indicates a stronger smoothing effect. G300/G500: a Gaussian filter with a 300km/500-km radius white/stochastic noise and suggest that this random signal can be identified as the time information grows. Although some temporal filters tried to exploit the prior information of correlation in the spherical harmonic order (e.g., the EOF and stochastic filters), they mostly depended on temporal information, rather than the particular spatial properties of the NSS noise. Both spatial and temporal filters have their own limitations: spatial filters, in order to ensure the quality of the monthly data, significantly underestimate the long-term trends; the performance of temporal filters depends heav-ily on the length and availability of the monthly GRACE products, so it is no surprise to find that the results of applying temporal filters are changing as the observation period expands. Hence, the temporal and spatial filters should not be regarded as complete substitutes for each other. Rather, they utilize information from different domains that can be combined to improve the data quality (e.g., Prevost et al. 2019). Nonetheless, this study only focuses on spatial type filters. To date, the de-striping method proposed by Swenson and Wahr (2006) (hereinafter SW) provides the only option for the second strategy. They found that NSS noise is associated with spurious correlations between specific coefficients. As a solution, SW's approach progressively took all coefficients in the same order as a series and fitted and removed polynomial functions for coefficients of odd or even sequences. Although an additional filter is usually required to suppress residual random noise and uncorrected NSS noise in the second strategy, the extra filter is theoretically weaker than the counterpart in the first strategy, so the quantity of preserved signals is hopefully larger. To corroborate this claim, we provide the example of the RL06 GSM product in December 2004 in Fig. 1. Two popular Gaussian and DDK filters are applied separately to suppress the NSS noise. The Gaussian filter with a 300-km radius (G300) can dampen most of the noise, but the residual NSS pattern indicates that a stronger filter such as G500 is needed. The same conclusion can be drawn when DDK6 and its stronger alternative DDK3 are applied. In order to demonstrate what signal is lost when the filter is strengthened, we calculate the difference between the results using DDK6 and DDK3 and smooth the difference by DDK2 to emphasize the signal damping in low degrees (Fig. 1b). The striking differences in Greenland, Amazon, Antarctica and Sumatra are definitely signals instead of noise. Therefore, in order to effectively suppress the overwhelming NSS noise, one has to use an excessively strong filter to simultaneously attenuate signals that can be kept if the NSS noise is smaller. This example highlights how the second strategy of filtering can potentially help improve the spatial resolution and reduce amplitude damping in the post-processing of GRACE data.
Although widely used and further developed (Chambers 2006;Chen et al. 2006;Duan et al. 2009), SW's method has two main limitations: (1) its potential distortion of physical signals is not well understood because of manipulating coefficients in the spectral domain and (2) it becomes inapplicable when the coefficient order approaches the maximum degree due to insufficient coefficients for the polynomial fitting. As shown later, the first limitation is nonnegligible, especially for coefficients of around degree 15 and for observations at high latitudes. However, the explicit rule of the distortions is unclear. The second limitation leaves coefficients of order over 50/86 for SHCs truncated at degree 60/96, where NSS noise abounds, unprocessed.
The user-friendly mascon products Watkins et al. 2015;Loomis et al. 2019) utilized other sophisticated techniques to suppress the noise and recover the signal intensity. The inversion of mascons makes use of a priori information on the parameter uncertainties that can be derived from ancillary geophysical observations/models, geographic boundaries, and temporal correlation. However, mascons products have two limitations compared to SHCs. First, mascons stem from the assumption that signals are located at the earth's surface, which is generally not applica-ble in the field of solid earth (Chao 2016), such as earthquakes (Han et al. 2013), crustal thickening (Yi et al. 2016), and post-glacial rebound (Peltier et al. 2015). Second, as the signal-to-noise ratio is different from place to place, the global regularization used in the mascons inversion does not always guarantee a regionally optimal solution. For example, mascon products are apt to remove small-sized geophysical signals that can be found in SHCs (Zhang et al. 2019). As a price for ease of use, mascons lack the flexibility to adjust the balance between signal and noise by users.
Here, we provide a new alternative that overcomes the limitations in SW's method for filters of the second strategy. The new spatial filtering method has many similarities with the multi-channel singular spectrum analysis (SSA), but implemented in the spatial domain instead of the temporal domain; hence, our method is named SSAS. Note that it is fundamentally different from another SSA study of ours (Yi and Sneeuw 2021), which was implemented in the temporal domain to fill data gaps in the GRACE products. We will first introduce the principles and procedures of the SSAS method in Sect. 2. The performance of our method is compared with SW's method P 3 M 10 , which means a third-degree polynomial fitting is carried out for coefficients of orders exceeding ± 9. The results of December 2004, in both spatial and spectral domains, are provided in Sect. 3, where we also compare the spatial resolution of results using various methods in the Himalayas and surroundings. In Sect. 4, we will discuss the potential signal distortion produced by our method.

Data and method
The following work is based on Release 06 GSM product provided by the Center for Space Research (CSR) at the University of Texas. The monthly data include measurements from both GRACE and GRACE-FO from April 2002 to August 2020, which are available, e.g., at https://isdc.gfzpotsdam.de. Conventionally, the low-degree terms are specially treated: the degree-1 terms are supplemented by ) and the C 2,0 term is replaced by the result of Satellite Laser Ranging (Cheng and Ries 2017). The results are anomalies relative to GGM05C (Ries et al. 2016).
Although NSS noise is well known, its definition is still vague. Here, we put forward an unambiguous definition of NSS noise: high-frequency periodic variations in the longitude direction that persist in the latitude direction. The definition describes a pattern that seldom exists in a real physical signal. Therefore, the goal of identifying NSS noise becomes the search for quasi-periodic signals along the longitude direction. The multi-channel SSA method is adept at extracting common variations among a batch of time series, so we utilize it to find shared periodic signals (i.e., NSS noise) in longitudinal series.
Since we define the NSS noise in the spatial domain, the measured time-varying global gravity field, expressed in SHCs, needs to be transformed into space observations, such as gravity disturbance: where indicates that only the time-varying part is concerned here; C nm and S nm are SHCs as a function of degree n and order m; G represents the gravitational constant; R and M are the mean radius and total mass of the earth, respectively; P nm represents the normalized associated Legendre functions; θ andφ are colatitude and longitude, respectively. Here, the beginning of the degree n 0 is set to 21, since we have known that noise mostly exists at high degrees. We tested the result with n 0 0 and found that the low degrees are mostly immune to our method even if they are included. However, the exclusion of low degrees can improve the performance at high degrees because the proportion of noise to be found is greater.
A maximum SHC degree of 60 corresponds with a sampling resolution of 3°on the sphere surface. Nonetheless, the sampling rate in the longitude direction is increased to 1°t o allow for a subtle shift in the series (shown later), while the sampling interval in the latitude direction remains at 3°. Therefore, g(φ, θ) is calculated to create a matrix of 360 by 60, which represents the number of observations along the longitude and latitude, respectively. The i-th row and jth column of the matrix are denoted by x j i , with1 ≤ i ≤ 360and1 ≤ j ≤ 60. Here, we use the superscript j to indicate the latitude index, not an exponent. For simplicity, we omit the index j and denote all observations at a latitude as x i , which compose an array X. The spatial sequence X is equivalent to the time series in the traditional SSA, but it differs in the construction of the trajectory matrix.
We first remove the mean value of the series. Due to the periodicity in the longitude direction, we know that x 1 and x 360 are spatially adjacent, so we can rearrange the series to [x 2 , . . . , x 360 , x 1 ] by shifting it forward by one element and moving the first one to the end. Likewise, the series can be shifted by an arbitrary number of elements without losing any information. We combine all these series with shifts less than M (in one direction) to create a 360 × M trajectory matrix Y : x 359 x 360 x 360 According to our experience, M should be a few times the period of the target signal. The period of NSS noise is around 10 degrees, so we empirically choose M 40 here. The trajectory matrix along the latitude θ j is denoted by Y j . We then align the matrices of 60 latitudes horizontally to build a new augmented trajectory matrix Z: where w j sin(θ j ) is a scalar weight to reduce the signal amplitude as the grid area decreases. The giant matrix Z is similar to the equivalent one in the multi-channel SSA method (Ghil et al. 2002). In the next step, Z is factorized by singular value decomposition (SVD): where is diagonal and U and V are orthonormal. Each column of U, U k , contains the spatial characteristics in the longitude direction and is sorted in descending order by the explained variance of all columns in Z. The matrix V can be reshaped to V m, j,k , where m (≤ M), j (≤ 60), and k (≤ 360) represent the number of shifts, latitude index, and mode index, respectively.
The Y matrix of the jth latitude, Y j , can be rebuilt by the summation of all modes S k : where S has the same size as Y ; the asterisk indicates that the entire dimension is used; K is the number of modes. Given a value of M greater than 6 (which likely happens), K will be 360. Note that we omit the j index from the S matrix, although S is latitude dependent. Therefore, we can formulate a new time series, termed reconstructed components (RCs), from the mean value of skew-diagonal elements: where (p,q) are the location index of S. The number of elements satisfying either of the conditions in the parentheses is M and all of them are used for the average. The idea is to use elements with the same location index in Eq. (2). The sum of all RCs is equal to the j-th latitude observation, but with the mean value removed: The original observations along the longitude have been decomposed into K modes, based on the self-correlation in the longitude direction and also cross-correlation in the latitude direction. According to our definition of NSS noise at the beginning of this section, its changing pattern repeatedly occurs in both directions, so it will be strengthened by the construction of the trajectory matrix Z (Eq. 3) and transformed into stable and dominant modes. Although strong non-noisy signals can also contribute to significant modes, they are physically less likely to be periodic and repeated along longitude and latitude (two exception signals in Sumatra and Greenland will be discussed later), so we only need to identify the NSS noise by searching for periodic modes.
In the correlation analysis of time-lagged/space-lagged series, modes that represent periodic components always occur in pairs. It can be explained by assuming a series of signals with a period of T . The correlation between the shifted and original series will be close to ± 1 when the shift is an even multiple of T /4, but nearly zero when an odd multiple of T /4. Therefore, the cluster of shifted periodic signals cannot be reconstructed by one mode, but by two similar modes in a comparable magnitude with a phase lag of T /4. Those series that are shifted by non-integral multiples of T /4 are explained by both modes. This is analogous to the notion that two trigonometric functions sin(x) and cos(x) are necessary to form arbitrary phase signals of the same period.
Therefore, this phenomenon of pairwise occurrence provides a way to identify periodic components. We can directly check the correlation between the columns of the matrix U (Eq. 4). Since the coupled U i and U j should have similar magnitudes (they come from the same signal source but with different phase shifts), they should be close in the column sequence of the matrix. Hence, we only need to compare adjacent modes. Considering that the mode strength could be quite similar, the coupled modes may occasionally be separated, so the modes satisfying i − j ≤ 2 are checked. We shift U j by the number of elements within [−t, t] (t is an integer, here it is set to 3, equal to a maximum period of 12 degrees; negative means a backward shift) and search for its extreme correlation with U i . If the maximum and minimum correlation values can exceed ± 0.9 (empirically determined), these two modes are regarded as coupled. Finally, all coupled modes are selected to reconstruct the periodic variation, i.e., the NSS noise, of the j-th latitude. Thus, Eq. (7) is updated: The search for coupled modes is carried out only in the first 20 columns of U, because the preceding columns are relatively more reliable and NSS noise is supposed to account for the main part after excluding coefficients of degrees ≤ 20.
We checked the RMS of the change in equivalent water height after applying the SASS filter (Fig. 2a). The impacts are strong in four regions, but caused by two different mechanisms. The effects are strong in Greenland, Antarctic, and Amazon because their signals are extremely strong, producing a remarkable net value even though the relative changes are negligible. Specifically, the impact on the long-term trend in Greenland and Antarctic is around 1%, and the case in Amazon will be discussed in Fig. 9d. In contrast, the large anomaly in Sumatra is due to the high geometric resemblance of its seismic signals to NSS noise ( Figure S1), so our patternbased method will partially remove the physical signals. The 2011 Tohoku and 2010 Chile earthquakes are not influenced by this problem ( Figure S1). In order to protect the seismic signals in Sumatra, we omit observations within the circle region to create a gap area. Although the relative impact in Greenland and Antarctic is trivial, we still suppress their ice melting signals since the polar regions are not the focus of our method (Fig. 8). Besides, these two regions have extremely high signal-to-noise ratios, so it is best to keep signals there intact.
Therefore, we add two extra amendments to improve the identification of NSS noise. First, we reduce the signals in Greenland and Antarctica by 90%. The attenuation of the known strong signals can (1) reduce the possibility of its mixture with intense NSS noise in middle and low latitudes and (2) increase the fraction of NSS noise in the whole observations. Second, we omit observations in the Sumatra region to protect the seismic signals. Hence, the SVD technique in Eq. (4) is no longer applicable, and we need a gap-filling method that can conjecture the NSS noise in the gap area. This is feasible according to the repeatability of NSS noise in both latitude and longitude directions. For this purpose, we adopt the iterative gap-filling method (Kondrashov and Ghil 2006) that was originally designed for the SSA technique, which performs two loops that iteratively update the missing data.
To sum it up, the SSAS method has the following steps: 1. Getting global gridded gravity disturbance based on Eq.
(1), where n 0 21; 2. Removing the average of each longitude series, suppressing polar signals by the weight mask shown in Fig. 2b, and setting grids in the gap region to NaN values; 3. Constructing the augmented trajectory matrix based on Eq.
(3), where M is suggested to be 40; 4. Decomposing the trajectory matrix by an SSA gap-filling method; 5. Finding and extracting the coupled modes in the first 20 columns of U; the residuals are returned to step-3; the loop from step-3 to step-5 exits when no more coupled modes are found. Fig. 2 a RMS of the change in EWH imposed by SASS. Due to the geometric similarity of the seismic signals in Sumatra (the circle) to the NSS noise, we create a gap mask within the circle to protect the signals there. We also apply two extra weights to reduce the drastic signals of ice mass changes in polar regions (marked by the ellipses) and b the final mask and weight map 6. Reconstructing the longitude series by using all coupled components (Eq. 8) and the series of all latitudes compose the NSS noise; 7. Getting global gravity grids based on Eq. (1), where n 0 0, and subtracting the extracted NSS noise from it to get the NSS-noise-removed gravity. 8. If necessary, the NSS-noise-removed gravity can be expanded into SHCs and applied with an additional filter to suppress the random noise (e.g., DDK6).
Omitting the gap mask will only cause the signals in Sumatra to be underestimated, so this step can be skipped if this region is not focused, then the SSA gap-filling method in step-4 can be simplified to the SVD method (Eq. 4). Our method is also applicable to regional studies, where only the necessary latitude matrices are used to build the Z matrix in Eq. (3). According to our experience, increasing the window size M and the value of t in the search for coupled modes will increase the amount of NSS noise to be extracted. Hence, the results will be cleaner and an extra filter is not always necessary, but the chance of misidentification will rise as well. Given that the majority of the NSS noise can be easily removed with insensitivity to the parameter combinations, we recommend identifying NSS noise conservatively and adopting an additional filter for the residual noise. Even so, the final results are quite insensitive to the choice of parameters ( Figure S2). In this study, we usually choose DDK6 as the supplementary filter, but any other filter is theoretically applicable.

Result in the spatial domain
We test our method in the product of December 2004 (as in Fig. 1), and the results are shown in Fig. 3. By using our SSAS method, we effectively extracted most of the NSS pattern variations, which are treated as NSS noise and removed. The residuals are attributed to physical signals plus random and uncorrected NSS noise (hereinafter we simply call it the signal component), which are much cleaner than before even without applying an extra filter. Although moderate NSS-type signals remain in oceans, the signal component smoothed by DDK6 is suitable for investigating signals in land or averaged in the whole ocean area. The NSS noise Fig. 3 Comparison of results using our method (SSAS) and P 3 M 10 in December 2004. The original signal is separated into signal and noise components by SSAS or P 3 M 10 , which are shown separately with and without the additional DDK filter. The DDK2 filter is applied to emphasize the impact of the method on coefficients of low degrees. Note the last row has a color scale of shorter range component with DDK6 only displays anomalies in the NSS pattern, implying what is improved by our method compared to the result in Fig. 1c. The extracted NSS noise is filled with noise from high degrees of coefficients, so it is possible that signals from medium degrees are spuriously extracted and hidden in them, but we cannot recognize it. The NSS noise with DDK2 highlights the contribution from medium degrees, and the result shows a tiny NSS pattern. It confirms that the NSS noise mainly comes from coefficients of high degrees and our method imposes few changes on coefficients of moderate degrees.
The results using the decorrelation filter P 3 M 10 are compared here. Note that the filtered result may be slightly different based on the choice of degree of the polynomial, the unchanged proportion of low-degree coefficients, and whether the polynomial fits the entire coefficients or in the windows (Duan et al. 2009). Coefficients with orders larger than 50 are untreated due to insufficient data for the polynomial fitting. As a result, these intact coefficients greatly contaminate the observations using P 3 M 10 if no other filters are applied. Therefore, an additional filter is always necessary for SW's method. Compared to the signal result using SSAS + DDK6, using P 3 M 10 + DDK6 has fewer residual NSS-pattern signals at high latitudes, but more at middle and low latitudes (e.g., in Australia and the Amazon). Observations at high latitudes have a much higher quality than these at lower latitudes, as a result of denser orbit coverage. So, a reasonable decorrelation filter should impose changes weakened with the increase in latitude. However, the noise component of P 3 M 10 with DDK6 shows an opposite tendency and the remarkable signal modification in Greenland is apparently a distortion. Besides, the noise with DDK2 emphasizes that this method greatly modifies coefficients around degree 15, which are not likely to contain that much NSS noise, so it shows another sign of signal distortion.

Result in the spectral domain
We have compared the results in the spatial domain. The comparison in the spectral domain can be informative and provides an insight into the source of the differences. Figure 4 shows the amplitude per coefficient for the December 2004 results. The original amplitude (a) shows unreasonably large values in high orders, which is the main cause of the NSS noise. As a result, the root mean square (RMS) of the amplitudes by the degree (termed amplitude per degree) of the original dataset shows a trend reversal after degree 30 (b). The SSAS filter can effectively reduce the strength in the high orders. The extracted NSS noise shows an order-dependent pattern (d), consistent with the previous understanding that NSS noise is reflected by correlation in odd/even coefficients of the same order (Swenson and Wahr 2006). Moreover, the amplitude pattern of the extracted NSS noise indicates that it exists mainly in orders larger than about 20. Coefficients of low orders have little contribution to the NSS pattern, so it is reasonable to implement P 3 M 10 with the polynomial fitting roughly beginning from order 10.
The amplitude per coefficient for the result using P 3 M 10 is mostly decreased relative to the original one in locations where the fit-and-remove procedure is applied. However, the near-sectoral coefficients (where n approximates m) are still over one order of magnitude stronger than the coefficients of the same order or degree. The difference to the original, given in (f), shows mostly random patterns throughout all the filtered coefficients. P 3 M 10 leads to a reduction in strength beginning from the degree of 15, which accounts for the large signal distortion in the last graph of Fig. 3. In contrast, the amplitude-per-degree curve of SSAS agrees with the original GRACE product in the first 30 degrees and an increasing suppression effect for coefficients of higher degrees. The phenomena shown here also exists in the results of other months. Therefore, compared with P 3 M 10 , SSAS tends to better preserve the signal energy at medium degrees and diminish noise at high degrees to a greater extent.

Regional research
To show our method's improvements in spatial resolution and signal preservation, we compare the results by using SSAS + DDK6, DDK6, and three mascon products: Goddard Space Flight Center (GSFC) mascons RL06 v1 (Loomis et al. 2019), CSR mascons RL06 v2 , Jet Propulsion Laboratory (JPL) mascons RL06 v2 (Watkins et al. 2015). The post-glacial rebound effect in SHCs is corrected by the ICE-6G_C model (Peltier et al. 2015). The study area is chosen as the Tibetan Plateau where the mass signals are complex but well-studied (Fig. 5): glacier and snow are abundant in May in the Nyainqentanglha Mountain (the circle sign) (Yi et al. 2020), and the groundwater has been rapidly depleting in northern India (plus) during the GRACE era (Rodell et al. 2009), and glacier mass in Pamir (triangle) is varying heterogeneously (Brun et al. 2017). As these signals occur near the earth's surface, the results are shown in the form of EWH. The anomalies from 2004 to 2016 exhibit the temporal evolution of the long-term trend and interannual variation. In terms of major signal anomalies, particularly the trends in glacier and water mass changes, the results of SSAS + DDK6 and DDK6 demonstrate a similar pattern, confirming that our method does not erase or fabricate signals. Nevertheless, the signals are rendered with different details and amplitudes by different results. The results with DDK6 hold some NSS-type signals (the red dashed curves) and the center of the Indian ground-  Fig. 6. Comparison of SSAS + DDK6 and DDK6 shows that the application of the SSAS filter only imposes trifling impact on the trend and variabilities, but it can significantly reduce spikes in the signals (e.g., the dashed circle in Fig. 6a; more results are given in Sect. 4.2). Conventionally, the DDK4 filter-stronger than DDK6-is preferable due to the severe NSS noise, but it inevitably underestimates the trend by 10-20% compared to DDK6. With the help of SSAS, it is possible to effectively suppress the NSS noise without significant loss of signals. Another merit of our method is its dynamic adjustability. Unlike filters that treat each monthly product indiscriminately, our method does not play a significant effect when the data quality is good, such as in May 2016 (Fig. 5). The comparison of maps and time series shows how our method can help improve spatial resolution and reduce signal damping.
The mascon product is dedicated to recover the signal strength, but their results are inconsistent in the magnitude and the location of regional extrema due to the difference in the regularization strategy. The GSFC mascons yield results quite similar to the SSAS + DDK6 ones in terms of signal patterns and magnitudes. The CSR mascons are apparently stronger than the other products (Fig. 6, Figure S3). The failure to remove some isolated outliers (the dashed circle in Fig. 6b) implies that the CSR mascons may have used a weaker constraint for noise suppression, and therefore, its signals are stronger. The relatively lower resolution of 3degree of the JPL mascons causes the groundwater depletion signal in North India to drift eastward, so its trend in Fig. 6a is exceptionally weak. Unlike the divergent results produced by different mascon solutions, SSAS is a mild filter and only introduces a tiny impact on the pattern and intensity of the original signals.

Test of simulated signals
We simulate a gravity field of January 2004 from the Global Land Data Assimilation System (GLDAS) Noah v2.1 model (Rodell et al. 2004) to test the distortion of our method on hydrological signals. The data are available at http://disc.sci. gsfc.nasa.gov. We also add five synthetic signals with regular geometries with sizes from 3°by 3°to 10°by 10°in the Pacific Ocean as testing benchmarks (Fig. 7a). The input observations in the form of EWH, composed of soil moisture, snow coverage, and canopy, are expanded into SHCs (Wahr et al. 1998). The values in Greenland could be extraordinarily large, so the values greater than 120 cm are truncated to 120 cm. We implement two different simulations without and with NSS noise extracted from the corresponding GRACE observations by our method, which are termed simulation-A and simulation-B, respectively. Since the NSS noise is derived from SSAS per se, we only apply P 3 M 10 to the noisefree simulation-A and apply SSAS to the simulation-B. The P 3 M 10 filter imposes a significant signal distortion mostly in the NSS pattern, which is an undesirable consequence of adjusting coefficients by the order. The remarkable deviation also corroborates our findings in Fig. 3, where coefficients of medium degrees and signals in high latitudes are vulnerable to distortions introduced by SW's method. We tested SSAS on simulation-A, but found no distortion at all (so not shown here), suggesting that our method has a good characteristic of orthogonality to physical signals. When NSS noise is added, our method can effectively separate the signal and noise components again, manifesting the robustness of the method.

Spatial distribution of the influence on real data
Here, using the CSR RL06 product, we calculate the RMS of the difference between the original and smoothed time series in each grid to check the spatial characteristic of the SSAS filter. The DDK3 and DDK6 filters are applied additionally to highlight the discrepancies in medium and high degrees, respectively. The impact of SSAS on the degree variance in each monthly product ( Figure S4) indicates that there is a significant difference in data quality between GRACE and GRACE Follow-on (GRACE-FO) products. Therefore, we show the results of GRACE (April 2002-June 2017) and GRACE FO (June 2018-August 2020) missions separately. The results shown in Fig. 8 demonstrate that the SSAS filter mostly changes grids between latitudes of ± 60°and seems to be mostly latitude-dependent. The absolute differences (left column) smoothed by DDK6 between ± 60°latitudes have a median value of 2.8 cm for GRACE and 1.3 cm for GRACE-FO, indicating the proportion of NSS noise is much reduced in GRACE-FO. The absolute differences with DDK3 are about 20% of these with DDK6, implying that most of the NSS noise found by our method is located in high-degree coefficients. The relative differences (right column) show that oceanic grids are mostly modified by about 60% and 20% in GRACE products with DDK6 and DDK3, respectively, while the land grids are apparently less adjusted, except in regions where the signals are always weak, like the Sahara Desert. The same conclusion can be drawn from the GRACE-FO results, but with a smaller degree of changes (35% and 10% for DDK6 and DDK3, respectively). Time series of four spots (A: Sahara; B: Australia; C: Pacific; D: Amazon) with their RMS values are shown in Fig. 9. Mass changes in the Sahara and Pacific are mostly tiny, so a smaller signal RMS is generally preferable. Mass changes in Australia are usually small as well, but El Nino-Southern Oscillation (ENSO) events may occasionally bring tremendous precipitation. The well-known 2011 La Nina event (Boening et al. 2012) and associated rainy seasons are better captured by the SSAS + DDK6 series, which shows a continuous progression from wet to recovery period in Australia in 2011. The adoption of the SSAS filter generally reduces the RMS by~4 cm and can effectively remove isolated outliers in the time series of A, B, and C. The time series of D has a tremendous seasonal variation, and our method preserves this feature well. The series of the differentials is scrutinized, and it shows random variations without a significant trend or seasonal variation. Therefore, our method exhibits reasonable behavior in terms of noise removal and signal maintenance.

Potential distortion of long-term trends
Lastly, we test the impact of our method on long-term trends. Note in step 2 of our method, we have protected signals in Sumatra and diminished these in Greenland and Antarctic. If not, the trend in Sumatra will sometimes be partly absorbed into the noise component and thus underestimated. This is the only region known to be vulnerable to the SASS method without applying the masks in Fig. 2. The long-term trend results for the periods 2003-2010 and 2011-June 2017 are presented in Fig. 10. Since the trend's signal-to-noise ratio is improved relative to the monthly product, a weaker DDK7 filter is used. The differences in trends show only weak or moderate NSS patterns, depending on the data quality, confirming that our method does not skew the long-term trends. This is a significant improvement over traditional spatial filters, which inevitably underestimate trends if the data quality at the monthly scale is guaranteed.
As the most reliable component in the time series, the annual variation is even less distorted ( Figure S5). The impact in global basins is quantified in Figures S6-8 and Table S1. Results at the basin scale are more consistent than those at Fig. 7 Test of distortion of SSAS and P 3 M 10 with simulated signals. Simulation-A (a) is noise-free. Simulation-B (b) is the sum of simulation-A and NSS noise extracted from GRACE data using SSAS. c Application of P 3 M 10 to simulation-A (noise-free) and the residuals are shown in (e). P 3 M 10 is not applied to simulation-B since its signal distortion is already significant. d Application of SSAS to simulation-B and the difference relative to simulation-A is shown in (f). SSAS is completely orthogonal to simulation-A, so the result is not shown. The simulated signal comes from the GLDAS land surface model, with five synthetic signal sources in regular shapes in the Pacific Ocean the grid scale because the average within the basin offsets the positive/negative differences to some extent. In terms of trends, the difference between SSAS and DDK7 is up to 0.25 cm/yr, with about 88% of basins having a difference less than 0.1 cm/yr. The relative RMS value between SSAS and DDK6 is smaller than 0.2 for more than half of the basins.

Conclusions
Here, we present a new SSAS spatial filtering method dedicated to removing NSS noise, which is defined as high-frequency periodic signals prevailing in the longitude direction and persistent in the latitude direction. The SSAS method is developed based on the multi-channel SSA technique but in the spatial domain, and incidentally benefit from the geographic periodicity, which avoids the truncation of embedded series due to the lag window. The pattern-based method may occasionally misidentify signals in Sumatra as NSS noise due to the geometrical similarity. Therefore, to overcome the known distortion, we leave out signals in Sumatra to create a gap area, which can be filled by the SSA gap-filling method. It demonstrates that our method is adjustable.
Results in December 2004 are shown in the spatial and spectral domains. The NSS noise extracted using SSAS has the following features: (1) it only holds the NSS pattern; (2) it lies mainly between latitudes of ± 60°; (3) it is mainly located Fig. 8 Global distribution of absolute (left) and relative (right) changes imposed by SASS. Two DDK filters are used to highlight the impacts on low/high degrees. The results of GRACE (upper four maps) and GRACE-FO (lower four) are shown separately. The green circles (labeled in graph-h) mark four spots with time series compared in Fig. 9 in coefficients of orders greater than 20 and degrees greater than 30; (4) it seems to be order dependent and its amplitude per degree increases as the degree increases. These are consistent with the understanding of NSS noise and add credibility to our method. Detailed examination in the Himalayan region suggests that results by using SSAS can improve the spatial resolution and signal preservation of GRACE data.
Apart from noise removal, an equally important criterion for a filter is potential signal distortion. We tested SSAS on simulated signals and scrutinized spatial distribution of changes imposed by SSAS and time series at four spots and checked the bias in the long-term trend. It can be concluded that: Fig. 9 Time-series comparison of the four spots marked in Fig. 8. The values represent the RMS of the individual series 1. real physical signals, especially hydrological ones, are generally immune to the SSAS filter; 2. NSS noise extracted by the method can be correctly found again if mixed with simulated signals; 3. the NSS noise with DDK6 is~2.8 cm for GRACE and~1.3 cm for GRACE-FO in regions between latitudes of ± 60°; 4. SSAS with DDK6 can reduce variations in oceanic signals by 60% and 35% for GRACE and GRACE-FO, respectively; 5. SSAS can apparently reduce isolated outliers in time series, although no temporal constraint is applied; 6. SSAS puts little impact on long-term trends.
A filter can be evaluated in two aspects: (1) the effectiveness of removing target noise and (2) the probability of mistaking signal for noise. As for the former, SW's method cannot handle high-order coefficients. Our method can effectively address this issue, but cannot find the NSS noise present at latitudes higher than ± 60°. We consider that this is not a big problem because the signal-noise-ratio is greatly improved there. Regarding the latter, SW's method easily modifies physical signals (starting from SHC degree 15). Our method shows better orthogonality to physical signals, but is prone to temper signals in Sumatra. We introduce the gap filling method to remedy this limitation.
As a new alternative for filters of the second strategy (that removes NSS noise before filtering), the SSAS method does not have the limitations of signal distortion and incomplete correction rooted in SW's method. We encourage GRACE data users to adopt this type of filter because it directly mitigates the strongest interference noise with little signal attenuation, relieving the burden of other spatial or temporal filters, thus increasing the chance of identifying more and subtler signals.