Journal of Geodesy (2022) 96:23 https://doi.org/10.1007/s00190-022-01614-z ORIG INAL ART ICLE A novel spatial filter to reduce north–south striping noise in GRACE spherical harmonic coefficients Shuang Yi1 · Nico Sneeuw1 Received: 28 July 2021 / Accepted: 8 March 2022 / Published online: 4 April 2022 © The Author(s) 2022 Abstract 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. Keywords GRACE · Striping noise · Filter · Spatial filtering · Singular spectra analysis 1 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 fre- quency. 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 spa- tial filter that makes the most of the geometric pattern of NSS noise. Each element of SHCs is assumed to contain three compo- nents: 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. B Shuang Yi shuang.yi@gis.uni-stuttgart.de 1 Institute of Geodesy, University of Stuttgart, 70174 Stuttgart, Germany Generally, there are two strategies to achieve this goal: (i) applying the filter to the coefficient directly and (ii) propos- ing a pattern to disaggregate the NSS noise component and impose an additional filter to suppress residual noise. The implementation of the first strategy canmainly be sep- arated 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 DDKfilter (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 homoge- neous 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 (Schrama et al. 2007; 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 123 http://crossmark.crossref.org/dialog/?doi=10.1007/s00190-022-01614-z&domain=pdf http://orcid.org/0000-0003-2976-2351 http://orcid.org/0000-0003-1796-0131 23 Page 2 of 17 S. Yi, N. Sneeuw Fig. 1 GRACE observations of global gravity disturbance (µGal) in December 2004 without (a) and with different filters (c-f). (DDK6- DDK3)×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 300- km/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 infor- mation, 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 apply- ing 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 com- bined 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 123 A novel spatial filter to reduce north–south striping noise in GRACE spherical harmonic… Page 3 of 17 23 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 poly- nomial 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 sec- ond 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 pro- vide 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 andSumatra are definitely signals instead of noise. Therefore, in order to effectively suppress the overwhelming NSSnoise, one has to use an excessively strongfilter to simul- taneously 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 resolu- tion 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 coef- ficients 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 (Save et al. 2016; Watkins et al. 2015; Loomis et al. 2019) utilized other sophis- ticated 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 strat- egy. The new spatial filtering method has many similarities with themulti-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.Wewill first introduce the principles and procedures of the SSAS method in Sect. 2. The performance of our method is compared with SW’s method P3M10, which means a third-degree polyno- mial fitting is carried out for coefficients of orders exceed- ing±9. The results of December 2004, in both spatial and spectral domains, are provided in Sect. 3, wherewe also com- pare 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. 2 Data andmethod The following work is based on Release 06 GSM product provided by the Center for Space Research (CSR) at the Uni- versity 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.gfz- potsdam.de. Conventionally, the low-degree terms are spe- cially treated: the degree-1 terms are supplemented by (Sun et al. 2016) and the C2,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 lon- gitude 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 lon- gitude direction. The multi-channel SSA method is adept at extracting common variations among a batch of time series, sowe utilize it to find shared periodic signals (i.e., NSSnoise) in longitudinal series. 123 https://isdc.gfz-potsdam.de 23 Page 4 of 17 S. Yi, N. Sneeuw 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: �g(φ, θ) � GM R2 N∑ n�n0 n∑ m�0 (n + 1)Pnm(cosθ) ( �Cnmcosmφ + �Snmsinmφ ) , (1) where � indicates that only the time-varying part is con- cerned here; Cnm and Snm are SHCs as a function of degree n and orderm;G represents the gravitational constant; R and M are the mean radius and total mass of the earth, respec- tively; Pnm represents the normalized associated Legendre functions; θandφ are colatitude and longitude, respectively. Here, the beginning of the degree n0 is set to 21, sincewe have known that noise mostly exists at high degrees. We tested the result with n0 � 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 sam- pling resolution of 3° on the sphere surface. Nonetheless, the sampling rate in the longitude direction is increased to 1° to 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 j- th column of the matrix are denoted by x j i ,with1 ≤ i ≤ 360and1 ≤ j ≤ 60. Here, we use the superscript j to indi- cate the latitude index, not an exponent. For simplicity, we omit the index j and denote all observations at a latitude as xi , which compose an array X. The spatial sequence X is equiv- alent 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 x1 and x360 are spatially adjacent, so we can rearrange the series to [x2, . . . , x360, x1] 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 elementswithout losing any information.We combine all these series with shifts less than M (in one direction) to create a 360 × M trajectory matrix Y : Y360×M � ⎡ ⎢⎢⎢⎢⎢⎣ x1 x2 x2 x3 · · · xM xM+1 ... . . . ... x359 x360 x360 x1 · · · xM−2 xM−1 ⎤ ⎥⎥⎥⎥⎥⎦ . (2) 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 Yj. We then align the matrices of 60 latitudes horizontally to build a new augmented trajectory matrix Z: Z360×60M � [ w1Y 1 · · · w60Y 60 ] , (3) where w j � sin(θ j ) is a scalar weight to reduce the sig- nal 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): Z360×60M � U360×360�360×360V ′ 60M×360, (4) where � is diagonal and U and V are orthonormal. Each column of U, Uk , 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 Vm, 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 Sk : Y j � K∑ k�1 Sk � K∑ k�1 U∗,k × σk,k × V ′∗, j,k, (5) 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: (6) RCk i � mean ( Skp,q ) , (p + q � i + 1orp + q � i + 360; 1 ≤ i ≤ 360), where (p,q) are the location index of S. The number of ele- ments 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: X j � 1 w j K∑ k�1 RCk . (7) 123 A novel spatial filter to reduce north–south striping noise in GRACE spherical harmonic… Page 5 of 17 23 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 lat- itude 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 trans- formed 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 Suma- tra 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 sig- nals 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 can- not 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 neces- sary to form arbitrary phase signals of the same period. Therefore, this phenomenon of pairwise occurrence pro- vides a way to identify periodic components. We can directly check the correlation between the columns of the matrix U (Eq. 4). Since the coupled Ui and Uj 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 sepa- rated, so themodes satisfying i− j ≤ 2 are checked.We shift Uj 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 Ui . If the maximum and minimum cor- relation values can exceed±0.9 (empirically determined), these twomodes 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: X̂ j � 1 w j ∑ coupled RCk . (8) The search for coupled modes is carried out only in the first 20 columns of U, because the preceding columns are relativelymore reliable andNSSnoise 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 mech- anisms. The effects are strong in Greenland, Antarctic, and Amazon because their signals are extremely strong, produc- ing 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 toNSSnoise (Figure S1), so our pattern- based 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 obser- vations. 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 n0 � 21; 2. Removing the average of each longitude series, suppress- ing 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 trajectorymatrix 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. 123 23 Page 6 of 17 S. Yi, N. Sneeuw 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 n0 � 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 gapmaskwill only cause the signals inSuma- tra 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 sizeM 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 param- eters (Figure S2). In this study, we usually choose DDK6 as the supplementary filter, but any other filter is theoretically applicable. 3 Result 3.1 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 pat- tern 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 123 A novel spatial filter to reduce north–south striping noise in GRACE spherical harmonic… Page 7 of 17 23 Fig. 3 Comparison of results using our method (SSAS) and P3M10 in December 2004. The original signal is separated into signal and noise components by SSAS or P3M10, which are shown separately with and without the additional DDK filter. The DDK2 filter is applied to empha- size the impact of the method on coefficients of low degrees. Note the last row has a color scale of shorter range 123 23 Page 8 of 17 S. Yi, N. Sneeuw 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 P3M10 are com- pared 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 poly- nomial fitting. As a result, these intact coefficients greatly contaminate the observations using P3M10 if no other filters are applied. Therefore, an additional filter is always neces- sary for SW’s method. Compared to the signal result using SSAS + DDK6, using P3M10 + DDK6 has fewer residual NSS-pattern signals at high latitudes, but more at middle and low latitudes (e.g., in Australia and the Amazon). Observa- tions at high latitudes have amuch higher quality than these at lower latitudes, as a result of denser orbit coverage. So, a rea- sonable decorrelation filter should impose changesweakened with the increase in latitude. However, the noise component of P3M10 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. 3.2 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. Fig- ure 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 SSASfilter 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 P3M10 with the polynomial fitting roughly beginning from order 10. The amplitude per coefficient for the result using P3M10 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), showsmostly random patterns throughout all the filtered coefficients. P3M10 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 origi- nal GRACE product in the first 30 degrees and an increasing suppression effect for coefficients of higher degrees. The phe- nomena shown here also exists in the results of other months. Therefore, compared with P3M10, SSAS tends to better pre- serve the signal energy atmediumdegrees and diminish noise at high degrees to a greater extent. 3.3 Regional research To showourmethod’s improvements in spatial resolution and signal preservation, we compare the results by using SSAS + DDK6, DDK6, and three mascon products: Goddard Space FlightCenter (GSFC)masconsRL06v1 (Loomis et al. 2019), CSRmasconsRL06v2 (Save et al. 2016), Jet PropulsionLab- oratory (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 hetero- geneously (Brun et al. 2017). As these signals occur near the earth’s surface, the results are shown in the formofEWH.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 andDDK6 demonstrate a similar pattern, confirming that our method does not erase or fabricate signals. Nevertheless, the signals are renderedwith different details and amplitudes bydifferent results. The results with DDK6 hold some NSS-type signals (the red dashed curves) and the center of the Indian ground- 123 A novel spatial filter to reduce north–south striping noise in GRACE spherical harmonic… Page 9 of 17 23 Fig. 4 Comparison of amplitude per degree (b) and amplitude per coefficient (others) for the December 2004 results. Shaded areas indicate coefficients that are kept unchanged by the method. Note d and f are differences of logarithm values, which represent differences in orders of magnitude water signal in 2004 is thus shifted westward. The time series at three grid points located in North India, Nyainqentanglha Mountain and Pamir are examined in Fig. 6. Comparison of SSAS + DDK6 and DDK6 shows that the application of the SSASfilter only imposes trifling impact on the trend and vari- abilities, 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 indiscrimi- nately, 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 damp- ing. 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 123 23 Page 10 of 17 S. Yi, N. Sneeuw Fig. 5 Mass changes using two different filtering strategies and three mascon products in the Himalayas and surroundings in May of years of 2004, 2010, and 2016. Nyainqentanglha Mountain (circle) North India (plus), and Pamir (triangle) are marked. The results are anomalies rela- tive to the respective average for the period 2003–2010. The red dashed curves indicate likely NSS noise remaining in the DDK6 results 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 fail- ure 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 3- degree 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. 123 A novel spatial filter to reduce north–south striping noise in GRACE spherical harmonic… Page 11 of 17 23 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 Time (year) -200 -150 -100 -50 0 50 W at er h ei gh t ( cm ) (a) North India (76oE, 29oN) GSFC mascon: -6.29 ± 0.10 CSR mascon: -7.40 ± 0.20 JPL mascon: -4.00 ± 0.08 No filter: -6.93 ± 0.25 DDK6: -6.20 ± 0.13 SSAS+DDK6: -5.99 ± 0.12 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 Time (year) -150 -100 -50 0 50 W at er h ei gh t ( cm ) (b) Nyainqentanglha Mountain (95oE,30oN) GSFC mascon: -3.01 ± 0.06 CSR mascon: -4.68 ± 0.21 JPL mascon: -3.51 ± 0.14 No filter: -4.41 ± 0.25 DDK6: -2.95 ± 0.10 SSAS+DDK6: -2.91 ± 0.09 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 Time (year) -80 -60 -40 -20 0 20 40 W at er h ei gh t ( cm ) (c) Pamir (71oE,38oN) GSFC mascon: -1.35 ± 0.07 CSR mascon: -0.94 ± 0.10 JPL mascon: -1.58 ± 0.14 No filter: -1.93 ± 0.21 DDK6: -1.42 ± 0.11 SSAS+DDK6: -1.37 ± 0.10 Fig. 6 Time series of mass variations in North India (a), Nyainqentanglha Mountain (b) and Pamir (c) using various methods. The series show anomalies relative to the average in 2003. The trend (cm/yr) over April 2002–August 2020 is annotated 123 23 Page 12 of 17 S. Yi, N. Sneeuw 4 Discussion on signal distortion 4.1 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 reg- ular 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 soilmoisture, snow coverage, and canopy, are expanded into SHCs (Wahr et al. 1998). The values in Greenland could be extraordinar- ily large, so the values greater than 120 cm are truncated to 120 cm.We implement two different simulationswithout 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 P3M10 to the noise- free simulation-A and apply SSAS to the simulation-B. The P3M10 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 vulner- able 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 charac- teristic 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. 4.2 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 inGRACE-FO.The absolute differenceswithDDK3 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 GRACEproductswithDDK6 andDDK3, 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 Aus- tralia in 2011. The adoption of the SSAS filter generally reduces the RMS by~4 cm and can effectively remove iso- lated outliers in the time series of A, B, andC. 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 sig- nificant trend or seasonal variation. Therefore, our method exhibits reasonable behavior in terms of noise removal and signal maintenance. 4.3 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, con- firming that our method does not skew the long-term trends. This is a significant improvement over traditional spatial fil- ters, 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 123 http://disc.sci.gsfc.nasa.gov A novel spatial filter to reduce north–south striping noise in GRACE spherical harmonic… Page 13 of 17 23 Fig. 7 Test of distortion of SSAS and P3M10 with simulated sig- nals. 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 P3M10 to simulation-A (noise-free) and the residuals are shown in (e). P3M10 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 andDDK6 is smaller than 0.2 formore than half of the basins. 5 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 tech- nique 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 theNSS pattern; (2) it liesmainly between latitudes of±60°; (3) it ismainly located 123 23 Page 14 of 17 S. Yi, N. Sneeuw 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 ampli- tude per degree increases as the degree increases. These are consistent with the understanding ofNSS noise and add cred- ibility 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: 123 A novel spatial filter to reduce north–south striping noise in GRACE spherical harmonic… Page 15 of 17 23 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. NSSnoise extracted by themethod 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 forGRACE-FO in regions between latitudes of±60°; 4. SSAS with DDK6 can reduce variations in oceanic sig- nals 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 effec- tiveness 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 effec- tively 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 fromSHCdegree 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 123 23 Page 16 of 17 S. Yi, N. Sneeuw Fig. 10 Impact of SSAS on long-term trends for the periods 2003–2010 and 2011–June 2017 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. Supplementary Information The online version contains supplemen- tary material available at https://doi.org/10.1007/s00190-022-01614-z. Acknowledgements S.Y. was supported by the Alexander von Hum- boldt Foundation. Authors’ contributions S.Y. conceived, designed, and programmed the study. S.Y. and N.S. wrote the manuscript. Funding Open Access funding enabled and organized by Projekt DEAL. Data availability The GRACE data that support this study are available at the International Center for Global Earth Models (http://icgem.gfz- potsdam.de). The MATLAB code implementing the method and the CSR RL06 product smoothed by SSAS are deposited at https://github. com/shuang-yi/SSAS-GRACE-filter. Declarations Conflict of interest The authors declare no competing financial inter- ests. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adap- tation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indi- cate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, youwill need to obtain permission directly from the copy- right holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. 123 https://doi.org/10.1007/s00190-022-01614-z http://icgem.gfz-potsdam.de https://github.com/shuang-yi/SSAS-GRACE-filter http://creativecommons.org/licenses/by/4.0/ http://creativecommons.org/licenses/by/4.0/ A novel spatial filter to reduce north–south striping noise in GRACE spherical harmonic… Page 17 of 17 23 References Boening C, Willis J, Landerer F, Nerem R, Fasullo J (2012) The 2011 La Niña: so strong, the oceans fell. Geophys Res Lett. https://doi. org/10.1029/2012GL053055 Brun F, Berthier E, Wagnon P, Kääb A, Treichler D (2017) A spatially resolved estimate of highmountain asia glaciermass balances from 2000 to 2016. Nat Geosci 10(9):668–673. https://doi.org/10.1038/ ngeo2999 Chambers DP (2006) Evaluation of new GRACE time-variable gravity data over the ocean. Geophys Res Lett. https://doi.org/10.1029/ 2006GL027296 Chao BF (2016) Caveats on the equivalent water thickness and sur- facemascon solutions derived from theGRACESatellite-observed time-variable gravity. J Geodesy 90:807–813 Chen JL, Wilson CR, Seo K-W (2006) Optimized smoothing of gravity recovery and climate experiment (GRACE) time-variable gravity observations. J Geophys Res Solid Earth. https://doi.org/10.1029/ 2005JB004064 Cheng M, Ries J (2017) The unexpected signal in GRACE estimates of C20. J Geodesy 91:897–914 Crowley JW, Huang J (2020) A least-squares method for estimat- ing the correlated error of GRACE models. Geophys J Int 221(3):1736–1749. https://doi.org/10.1093/gji/ggaa104 Devaraju B, Sneeuw N (2017) The polar form of the spherical har- monic spectrum: implications for filtering grace data. J Geodesy 91(12):1475–1489. https://doi.org/10.1007/s00190-017-1037-7 Duan XJ, Guo JY, Shum CK, van der Wal W (2009) On the postpro- cessing removal of correlated errors in GRACE temporal gravity field solutions. J Geodesy 83(11):1095–1106. https://doi.org/10. 1007/s00190-009-0327-0 GhilM,AllenMR,DettingerMD, IdeK,KondrashovD,MannME et al (2002) Advanced spectral methods for climatic time series. Rev Geophys 40(1):3-1-3–41. https://doi.org/10.1029/2001RG000092 Han S-C, Shum CK, Jekeli C, Kuo C-Y, Wilson C, Seo K-W (2005) Non-isotropic filtering of GRACE temporal gravity for geophysi- cal signal enhancement. Geophys J Int 163(1):18–25. https://doi. org/10.1111/j.1365-246X.2005.02756.x Han S-C, Riva R, Sauber J, Okal E (2013) Source parameter inversion for recent great earthquakes from a decade-long observation of global gravity fields. J Geophys Res Solid Earth 118:1240–1267. https://doi.org/10.1002/jgrb.50116 KleesR, Revtova EA,Gunter BC,Ditmar P,OudmanE,WinsemiusHC, Savenije HHG (2008) The design of an optimal filter for monthly GRACE gravity models. Geophys J Int 175(2):417–432. https:// doi.org/10.1111/j.1365-246X.2008.03922.x Kondrashov D, Ghil M (2006) Spatio-temporal filling of missing points in geophysical data sets. Nonlinear Process Geophys 13(2):151–159. https://doi.org/10.5194/npg-13-151-2006 Kusche J (2007) Approximate decorrelation and non-isotropic smooth- ing of time-variable GRACE-type gravity field models. J Geodesy 81(11):733–749. https://doi.org/10.1007/s00190-007-0143-3 Kusche J, Schmidt R, Petrovic S, Rietbroek R (2009) Decorrelated GRACE time-variable gravity solutions by GFZ, and their valida- tion using a hydrological model. J Geodesy 83(10):903–913 Loomis BD, Luthcke SB, Sabaka TJ (2019) Regularization and error characterization of GRACE Mascons. J Geodesy 93:1381–1398 Peltier WR, Argus DF, Drummond R (2015) Space geodesy con- strains ice age terminal deglaciation: the global ICE-6G_C (VM5a) model. J Geophys Res Solid Earth 120:450–487 Prevost P, Chanard K, Fleitout L, Calais E, Walwer D, van Dam T, Ghil M (2019) Data-adaptive spatio-temporal filtering of GRACE data. Geophys J Int 219(3):2034–2055. https://doi.org/10.1093/ gji/ggz409 Ries J, Bettadpur S, Eanes R, Kang Z, Ko U, McCullough C, Nagel P, Pie N, Poole S, Richter T, Save H, Tapley B (2016) The combined gravity model GGM05C. GFZ Data Services. https://doi.org/10. 5880/ICGEM.2016.002 Rodell M, Houser P, Jambor U, Gottschalck J, Mitchell K, Meng C et al (2004) The global land data assimilation system. Bull AmMeteor Soc 85(3):381–394. https://doi.org/10.1175/bams-85-3-381 Rodell M, Velicogna I, Famiglietti JS (2009) Satellite-based estimates of groundwater depletion in India. Nature 460(7258):999–1002. https://doi.org/10.1038/nature08238 Sasgen I, Martinec Z, Fleming K (2006) Wiener optimal filtering of GRACE data. Stud Geophys Geod 50(4):499–508. https://doi.org/ 10.1007/s11200-006-0031-y Save H, Bettadpur S, Tapley BD (2016) High-resolution CSR GRACE RL05 Mascons. J Geophys Res Solid Earth 121:7547–7569 Schrama EJO,Wouters B, Lavallée DA (2007) Signal and noise in grav- ity recovery and climate experiment (GRACE) observed surface mass variations. J Geophys Res 112(B8):B08407. https://doi.org/ 10.1029/2006JB004882 Sun Y, Riva R, Ditmar P (2016) Optimizing estimates of annual varia- tions and trends in geocenter motion and J2 from a combination of GRACE data and geophysical models. J Geophys Res Solid Earth 121:8352–8370 Swenson S, Wahr J (2006) Post-processing removal of correlated errors in GRACE data. Geophys Res Lett. https://doi.org/10.1029/ 2005gl025285 Wahr J, Molenaar M, Bryan F (1998) Time variability of the earth’s gravity field: hydrological and oceanic effects and their possible detection using GRACE. J Geophys Res-Solid Earth 103(B12):30205–30229. https://doi.org/10.1029/98jb02844 Wang L, Davis JL, Hill EM, Tamisiea ME (2016) Stochastic filtering for determining gravity variations for decade-long time series of GRACE gravity. J Geophys Res Solid Earth 121(4):2915–2931. https://doi.org/10.1002/2015JB012650 Watkins MM et al (2015) Improved methods for observing earth’s time variable mass distribution with GRACE using spherical cap mas- cons: improved gravity observations fromGRACE. JGeophys Res Solid Earth 120:2648–2671 Wouters B, Schrama EJO (2007) Improved accuracy of GRACE gravity solutions through empirical orthogonal function filtering of spher- ical harmonics: EOF filtering of grace coefficients. Geophys Res Lett. https://doi.org/10.1029/2007GL032098 Yi S, Freymueller JT, Sun W (2016) How fast is the middle-lower crust flowing in Eastern Tibet? A constraint from geodetic observations: flow in Eastern Tibet. J Geophys Res Solid Earth 121:6903–6915 Yi S, Sneeuw N (2021) Filling the data gaps within GRACE mis- sions using singular spectrum analysis. J Geophys Res Solid Earth. https://doi.org/10.1029/2020JB021227 Yi S, Song C, Heki K, Kang S, Wang Q, Chang L (2020) Satellite- observed monthly glacier and snow mass changes in Southeast Tibet: implication for substantial meltwater contribution to the Brahmaputra. Cryosphere 14(7):2267–2281. https://doi.org/10. 5194/tc-14-2267-2020 Zhang L et al (2019) Evaluation of GRACE Mascon solutions for small spatial scales and localized mass sources. Geophys J Int 218(2):1307–1327 Zhang Z-Z, Chao BF, Lu Y, Hsu H-T (2009) An effective filtering for GRACE time-variable gravity: fan filter. Geophys Res Lett 36(17):L17311. https://doi.org/10.1029/2009GL039459 123 https://doi.org/10.1029/2012GL053055 https://doi.org/10.1038/ngeo2999 https://doi.org/10.1029/2006GL027296 https://doi.org/10.1029/2005JB004064 https://doi.org/10.1093/gji/ggaa104 https://doi.org/10.1007/s00190-017-1037-7 https://doi.org/10.1007/s00190-009-0327-0 https://doi.org/10.1029/2001RG000092 https://doi.org/10.1111/j.1365-246X.2005.02756.x https://doi.org/10.1002/jgrb.50116 https://doi.org/10.1111/j.1365-246X.2008.03922.x https://doi.org/10.5194/npg-13-151-2006 https://doi.org/10.1007/s00190-007-0143-3 https://doi.org/10.1093/gji/ggz409 https://doi.org/10.5880/ICGEM.2016.002 https://doi.org/10.1175/bams-85-3-381 https://doi.org/10.1038/nature08238 https://doi.org/10.1007/s11200-006-0031-y https://doi.org/10.1029/2006JB004882 https://doi.org/10.1029/2005gl025285 https://doi.org/10.1029/98jb02844 https://doi.org/10.1002/2015JB012650 https://doi.org/10.1029/2007GL032098 https://doi.org/10.1029/2020JB021227 https://doi.org/10.5194/tc-14-2267-2020 https://doi.org/10.1029/2009GL039459 A novel spatial filter to reduce north–south striping noise in GRACE spherical harmonic coefficients Abstract 1 Introduction 2 Data and method 3 Result 3.1 Result in the spatial domain 3.2 Result in the spectral domain 3.3 Regional research 4 Discussion on signal distortion 4.1 Test of simulated signals 4.2 Spatial distribution of the influence on real data 4.3 Potential distortion of long-term trends 5 Conclusions Acknowledgements References