Laurent Longuevergne

Y Wada

M Khaki

E Forootan

M Kuhn

J Awange

Efficient Basin Scale Filtering of GRACE Satellite Products

Keywords: GRACE Filter, River basin, Leakage reduction, Hydrology, Hydrological fluxes

Efficient basin scale filtering of GRACE satellite products

Introduction

Since 2002, the Gravity Recovery And Climate Experiment (GRACE) satellite mission has been providing time-variable global gravity field solutions [START_REF] Tapley | GRACE measurements of mass variability in the Earth system[END_REF]. These variations are primarily caused by temporal changes in the gravity field due to changes in hydrology, ice masses of the cryosphere, or surface deformation, e.g., glacial isostatic adjustment (GIA). Within a temporal and spatial resolution of respectively one day to one month and a few hundred kilometers, GRACE products have proved to be very useful for various geophysical and hydrological studies (see, e.g., [START_REF] Kusche | Mass distribution and mass transport in the Earth system[END_REF]Wouters et al., 2014, for applications).

In particular, the so-called level 2 (L2) time-variable gravity fields are widely used to quantify global (e.g., [START_REF] Rodell | The global land data assimilation system[END_REF][START_REF] Eicker | Does GRACE see the terrestrial water cycle 'intensifying[END_REF][START_REF] Kusche | Mapping probabilities of extreme continental water storage changes from space gravimetry[END_REF] and regional (e.g., [START_REF] Chen | 2005 drought event in the Amazon River basin as measured by GRACE and estimated by climate models[END_REF][START_REF] Awange | Water storage changes and climate variability within the Nile Basin between 2002 and 2011[END_REF][START_REF] Munier | Combining data sets of satellite-retrieved products for basin-scale water balance study: 2. Evaluation on the Mississippi Basin and closure correction model[END_REF]Khaki et al., 2017a,b) terrestrial total water storage (TWS) changes, i.e., the sum of changes in surface and sub-surface water storage compartments. GRACE products are also applied to estimate changes of the terrestrial water cycle (e.g., [START_REF] Ogawa | Acceleration signal in GRACE time-variable gravity in relation to interannual hydrological changes[END_REF][START_REF] Eicker | Does GRACE see the terrestrial water cycle 'intensifying[END_REF] or to validate the water cycle in atmospheric reanalyses (e.g., [START_REF] Springer | New estimates of variations in water ux and storage overEurope based on regional (Re) analyses and multisensor observations[END_REF][START_REF] Kusche | Mapping probabilities of extreme continental water storage changes from space gravimetry[END_REF][START_REF] Forootan | Large-scale total water storage and water flux changes over the arid and semi-arid 38 parts of the Middle East from GRACE and reanalysis products[END_REF].

Combined with information observed from other monitoring techniques (e.g., GPS and satellite altimetry) or simulations by land surface models, L2 products are applied to estimate surface (e.g., lakes and rivers) and subsurface (e.g., soil moisture and groundwater) storage changes at (river) basin scales (e.g., [START_REF] Syed | Total basin discharge for the Amazon and Mississippi River basins from GRACE and a landatmosphere water balance[END_REF][START_REF] Longuevergne | GRACE Hydrological estimates for small basins: Evaluating processing approaches on the High Plains Aquifer, USA[END_REF][START_REF] Famiglietti | Water in the balance[END_REF]Forootan et al., 2014b). GRACE L2 products are provided in terms of potential spherical harmonic coefficients, e.g., up to degree and order 60 or 90, which mainly represent the large-to medium-scale (e.g., few hundred km) time-variable gravity changes. However, the L2 potential coefficients contain different types of errors. A part of these errors is related to colored/correlated noise due to the anisotropic spatial sampling of the mission, instrumental noise (K-band ranging system, GPS, and the accelerometer observations and star cameras), and temporal aliasing caused by the incomplete reduction of short-term mass variations by models [START_REF] Forootan | Comparisons of atmospheric data and reduction methods for the analysis of satellite gravimetry observations[END_REF](Forootan et al., , 2014a;;[START_REF] Dobslaw | Modeling of present-day atmosphere and ocean non-tidal dealiasing errors for future gravity mission simulations[END_REF]. These errors are manifested as north-south striping patterns in the spatial domain (e.g., gridded TWS products). The application of smoothing techniques with the primary aim of removing the stripes can lead to spatial leakages. The spatial averaging introduced by the smoothing kernels such as the Gaussian Kernel in Jekeli (1981) or non-Gaussian Kernels in [START_REF] Kusche | Approximate decorrelation and nonisotropic smoothing of time variable GRACEtype gravity field models[END_REF], results in spatial interference of mass anomalies. These leakage errors do not allow for perfect separation of gravity anomalies, e.g., between land and oceans, and limit the detection of small-scale hydrological signals. The accuracy of GRACE TWS estimation is very important for hydrological applications especially at the basin scale, e.g., to interpret redistribution of water storage or to indicate drought and flood patterns (e.g., [START_REF] Yeh | Remote sensing of ground water storage changes in Illinois using the Gravity Recovery and Climate Experiment (GRACE)[END_REF][START_REF] Longuevergne | GRACE Hydrological estimates for small basins: Evaluating processing approaches on the High Plains Aquifer, USA[END_REF][START_REF] Awange | Exploring hydrometeorological drought patterns over the Greater Horn of Africa (19792014) using remote sensing and reanalysis products[END_REF]. Therefore, better postprocessing of GRACE data must be applied to improve consistencies between various types of products that are usually used for studying the water cycle (e.g., [START_REF] Eicker | Does GRACE see the terrestrial water cycle 'intensifying[END_REF].

Different filtering methods have been proposed to reduce north-south striping errors, such as the isotropic Gaussian filter (Jekeli, 1981) and anisotropic filters (e.g., [START_REF] Swenson | Post-processing removal of correlated errors in GRACE data[END_REF][START_REF] Kusche | Approximate decorrelation and nonisotropic smoothing of time variable GRACEtype gravity field models[END_REF]Klees et al., 2008). A comprehensive review on filtering techniques has been done e.g., by [START_REF] Frappart | Monitoring Water Mass Redistributions on Land and Polar Ice Sheets using the GRACE Gravimetry from Space Mission[END_REF]. The isotropic Gaussian filter Jekeli ( 1981) is a degree-dependent filter in the spectral domain and bell-shaped filter in the spatial domain.

Anisotropic filters, on the other hand, are introduced to deal with the correlated errors between the coefficients of L2 products (e.g., different marginal shapes in the north-south and the eastwest directions). In general, filtering techniques that spatially smooth the L2 signal contents (e.g., [START_REF] Wahr | Accuracy of GRACE mass estimates[END_REF][START_REF] Kusche | Decorrelated GRACE time-variable gravity solutions by GFZ, and their validation using a hydrological model[END_REF] down-weight L2's higher degree and order potential coefficients. Although these filters reduce noises, their main problem is that they also attenuate the signals. In addition, the application of filtering moves gravity anomalies from one region to another region. Generally speaking, after applying a smoothing kernel some parts of the signals inside an area of interest leak out from it or alternatively signals from outside leak into the area of interest (e.g., [START_REF] Chen | Attenuation effect on seasonal basin-scale water storage changes from GRACE time-variable gravity[END_REF][START_REF] Baur | GRACE-derived ice-mass variations over Greenland by accounting for leakage effects[END_REF]. These issues become more critical for basin-scale studies, especially where the sizes of the basins are small in comparison to the spatial resolution of GRACE (e.g., [START_REF] Yeh | Remote sensing of ground water storage changes in Illinois using the Gravity Recovery and Climate Experiment (GRACE)[END_REF][START_REF] Longuevergne | GRACE Hydrological estimates for small basins: Evaluating processing approaches on the High Plains Aquifer, USA[END_REF].

Several methods have been put forward to mitigate spatial leakage effects in TWS estimations from L2 products. These methods can largely be categorised into the following three groups (i) those that numerically estimate the leakages (leakage in and out) using the averaging kernels (e.g., [START_REF] Seo | Simulated estimation of hydrological loads from GRACE[END_REF][START_REF] Baur | GRACE-derived ice-mass variations over Greenland by accounting for leakage effects[END_REF][START_REF] Longuevergne | GRACE Hydrological estimates for small basins: Evaluating processing approaches on the High Plains Aquifer, USA[END_REF], (ii) those that are based on scaling factors derived from synthetic data (e.g., [START_REF] Landerer | Accuracy of scaled GRACE terrestrial water storage estimates[END_REF][START_REF] Long | Deriving scaling factors using a global hydrological model to restore GRACE total water storage changes for China's Yangtze River Basin[END_REF], and (iii) those that use inversion for simultaneous signal separation and leakage reduction (e.g., [START_REF] Wouters | Improved accuracy of GRACE gravity solutions through empirical orthogonal function filtering of spherical harmonics[END_REF][START_REF] Frappart | An independent Component Analysis approach for filtering continental hydrology in the GRACE gravity data[END_REF]Forootan et al., 2014b;[START_REF] Frappart | Monitoring Water Mass Redistributions on Land and Polar Ice Sheets using the GRACE Gravimetry from Space Mission[END_REF]. From the first group, [START_REF] Swenson | Methods for inferring regional surface-mass anomalies from Gravity Recovery and Climate Experiment (GRACE) measurements of time-variable gravity[END_REF] developed an isotropic kernel using a Lagrange multiplier filter to best balance signal and leakage errors over a basin of interest. A non-isotropic Gaussian filter proposed by [START_REF] Han | Non-isotropic filtering of GRACE temporal gravity for geophysical signal enhancement[END_REF] to improve spatial resolution during the filtering process also belongs to this group. In another effort, Harig and Simons (2015) used Slepian-function analysis to decrease leakage effects in Antarctica by maximizing signal energy concentration within the area of interest. The second category uses synthetic data, e.g., from land surface models (LSMs) or hydrological fluxes to derive scaling factors that can be multiplied by GRACE filtered products to recover the lost signals. In this approach, efforts are focused on the application of the same filtering techniques to the synthetic data (that is close enough to the signal contents of GRACE products). Basin-averaged or gridded scale factors are usually estimated as the solution of a least squares adjustment that compares data before and after application of the filter. [START_REF] Landerer | Accuracy of scaled GRACE terrestrial water storage estimates[END_REF] estimated gridded scaling factors for GRACE TWS anomalies to restore the signals lost after applying a regular smoothing filter (a Gaussian smoothing kernel). A similar study that uses a different spatial scale (basin averages) has been performed by [START_REF] Long | Deriving scaling factors using a global hydrological model to restore GRACE total water storage changes for China's Yangtze River Basin[END_REF] who estimated scale factors using a global hydrological model over the Yangtze River Basin in China. A possible drawback of this approach is its dependency on the reliability of the hydrological model used to estimate the desired scale factors. The inversion techniques in (iii) also require a prior information about mass changes within different storage compartments. The dependency of final signal separation results on these information has not been reported yet.

To address the above problems arising from the application of filtering methods, the present study proposes a new filtering method, Kernel Fourier Integration (KeFIn), which is designed to reduce both types of above-mentioned errors using a two-step algorithm. In the first step, the advantages of image processing techniques such as motion filters (e.g., Hichri et al., 2012;[START_REF] Zhang | An effective filtering for GRACE time-variable gravity: Fan filter[END_REF] are exploited to reduce the measurement noise and aliasing of unmodelled high-frequency mass variations. This attempt is designed to keep as much of the higher frequency information as possible. It should be mentioned here that, although the proposed KeFIn filter has less effect on high-frequency signals compared to the existing methods, some signal inferences still exist mainly due to the truncation of degree and order in L2 products. In the second step of the KeFIn filter, the leakage problem is mitigated using an anisotropic kernel to isolate the signals in the basin of interest. The main idea of this step is to combine the Fourier transform and basin kernel functions to increase the strength of the attenuated signals. It will be shown in the following that the KeFIn filter is suited to deal with basins of various shapes and sizes.

The primary objectives of this study is developing a filter for (i) dealing with colored/correlated noise of high-frequency mass variations (i.e., stripes); and (ii) reducing basin scale spatial leakage errors for hydrological applications. These objectives are addressed by introducing novel methodologies discussed in Section 3.1.1 and 3.1.2, respectively. The performance of the introduced filtering method (KeFIn) in terms of leakage reduction is compared with commonly used methods that deal with leakage problem from the basin averaging kernel and the model-based scaling factor groups. For this purpose, both real and synthetic data sets are employed. The purpose of using synthetic data is to provide a more accurate evaluation of the newly proposed method in comparison to existing methods (e.g., [START_REF] Seo | Simulated estimation of hydrological loads from GRACE[END_REF][START_REF] Chen | 2005 drought event in the Amazon River basin as measured by GRACE and estimated by climate models[END_REF].

Therefore, we generate synthetic data in 43 globally distributed basins and use them to examine the performance of the proposed KeFIn and other commonly used filters. These filters are further assessed using water flux observations in the context of the water balance equation (see Equation 1 in Section 2.3), as well as by comparisons with in-situ measurements.

Data

GRACE

Monthly GRACE L2 products along with their full error information are computed at the Technical University of Graz known as the ITSG-Grace2014 gravity field models [START_REF] Mayer-Gürr | ITSG-Grace2014: a new GRACE gravity field release computed in Graz[END_REF]. We use these products and their full covariance errors up to degree and order 60 covering the period 2002-2013 (https://www.tugraz.at/institute/ifg/downloads/gravity-fieldmodels/itsg-grace2014). Degree 1 coefficients are replaced with those estimated by [START_REF] Swenson | Estimating geocentervariations from a combination of GRACE and ocean model output[END_REF] to account for the movement of the Earth's centre of mass. Degree 2 and order 0 (C20) coefficients are replaced by those from Satellite Laser Ranging solutions owing to unquantified large uncertainties in this term (e.g., [START_REF] Chen | Attenuation effect on seasonal basin-scale water storage changes from GRACE time-variable gravity[END_REF]. We also account for the post glacial rebound by incorporating the corrections provided by [START_REF] Geruo | Computations of the viscoelastic response of a 3-D compressible Earth to surface loading: an application to Glacial Isostatic Adjustment in Antarctica and Canada[END_REF]. The L2 gravity fields are then converted to 1 • ×1 • TWS fields following the approach of [START_REF] Wahr | Time variability of the Earth's gravity field' Hydrological and oceanic effects and their possible detection using GRACE[END_REF].

To evaluate the filtering techniques, no smoothing filter is applied at this stage on GRACE L2 products.

Synthetic data

In order to assess the efficiency of different filtering methods considered in this study, they are applied on synthetic data whose advantage is the possibility to unambiguously estimate leakage errors since the applied post-processing techniques must replicate the synthetic input data. For this purpose, the world's 43 major river basins with diverse sizes and shapes located at different places around the Earth are chosen (see Figure 1). A large number of significantly different basins helps us to properly investigate the efficiency and reliability of the newly proposed KeFIn filter. For synthetic TWS data, a summation of monthly (1 • ×1 • ) soil moisture, snow, and the canopy water storage from the Global Land Data Assimilation System (GLDAS) NOAH [START_REF] Rodell | The global land data assimilation system[END_REF](Rodell et al., ) over 2003(Rodell et al., -2013 is used (http://giovanni.sci.gsfc.nasa.gov/) is used (http://giovanni.sci.gsfc.nasa.gov/). Following [START_REF] Wang | Spherical harmonics scaling[END_REF], the TWS fields are converted to potential spherical harmonic coefficients up to degree and order 120. Only those coefficients that are up to degree and order 60 are used to generate similar spectral content as the real GRACE L2 products. These data are perturbed by north-south striping errors using the full covariance matrix of ITSG-Grace2014 products. Using the Cholesky decomposition method, the monthly covariance matrices are split into their upper triangular and their conjugate transpose matrices. By multiplying each of the upper triangular matrices with a column of the unit random matrix, the GRACE-type realizations of monthly errors are generated (see, e.g., [START_REF] Forootan | Separation of global time-variable gravity signals into maximally independent components[END_REF][START_REF] Kusche | Mapping probabilities of extreme continental water storage changes from space gravimetry[END_REF]. GLDAS TWS outputs are also used to compute model-derived scale factors using forward modelling following [START_REF] Long | Deriving scaling factors using a global hydrological model to restore GRACE total water storage changes for China's Yangtze River Basin[END_REF]. These hydrological datasets have also been used to estimate gridded gain factors following [START_REF] Landerer | Accuracy of scaled GRACE terrestrial water storage estimates[END_REF]. Results of these filters will be compared to the KeFIn filtering approach (see Section 4.1).

Auxiliary data sets

Recently developed Mass Concentration blocks (mascons) data (http://grace.jpl.nasa.gov) provided by Jet Propulsion Laboratory (JPL) are used to analyze their correlation to our estimation from L2 products as shown in the Appendix. The monthly JPL RL05M Mascon solution is post-processed liquid water equivalent thickness data using a Coastline Resolution Improvement (CRI) filter to separate the land and ocean portions of mass [START_REF] Wiese | GRACE monthly global water mass grids NETCDF RELEASE 5[END_REF][START_REF] Watkins | Improved methods for observing Earth's time variable mass distribution with GRACE using spherical cap mascons[END_REF]. We apply land-grid-scaling coefficients provided with the data to water equivalent thicknesses in 1 • ×1 • spatial resolution. These filtered data are compared with the results of filters applied in this study.

In addition, the temporal derivative of filtered GRACE data, known as total (hydrological) water fluxes (TWF) is compared with measured precipitation (P ), Evapotranspiration (ET ), and surface water discharge (or runoff, R) through the water balance equation below:

dS/dt = T W F = P -ET -R, (1)
where the dS/dt represents TWF derived from the ITSG-Grace2014 products following the procedure in [START_REF] Eicker | Does GRACE see the terrestrial water cycle 'intensifying[END_REF]. The assessment in Equation 1 requires additional hydrological water flux measurements, which are not easily accessible globally. Eight river basins are selected to perform this assessment, i.e., the Amazon (South America), Mekong (Southeast Asia), Arkansas-White (North America), Ohio (North America), Lachlan (Australia), Namoi (Australia), Lower Mississippi (North America), and Macquarie-Bogan (Australia) basins. We et al., 2010;[START_REF] Schmidt | The Question of Communist Land Degradation: New Evidence from Local Erosion and Basin-Wide Sediment Yield in Southwest China and Southeast Tibet[END_REF].

Each dataset is associated with a level of uncertainty and varies for different basins due to the diverse climatological condition. A number of studies has investigated the validity of above observations over various basins, e.g., [START_REF] Cai | Evaluation of TRMM precipitation data over the Inland River Basins of Northwest China[END_REF], [START_REF] Yan | Evaluation of TRMM Precipitation Product for Meteorological Drought Monitoring in Hai Basin[END_REF], [START_REF] Awange | Exploring hydrometeorological drought patterns over the Greater Horn of Africa (19792014) using remote sensing and reanalysis products[END_REF] for TRMM, as well as [START_REF] Velpuri | A comprehensive evaluation of two MODIS evapotranspiration products over the conterminous United States: Using point and gridded FLUXNET and water balance ET[END_REF], [START_REF] Ramoelo | Validation of Global Evapotranspiration Product (MOD16) using Flux Tower Data in the African Savanna, South Africa[END_REF][START_REF] Miralles | The WACMOS-ET project Part 2: Evaluation of global terrestrial evaporation data sets[END_REF] for MODIS products. Precipitation errors highly depend on temporal and spatial resolution [START_REF] Chen | Assessing objective techniques for gauge-based analyses of global daily precipitation[END_REF]. Uncertainty in measuring precipitation over lands are smaller compared to oceans since satellite data are merged with in-situ stations that are distributed over the continents.

The major source of uncertainty in MOD16 is the misclassification of landcover types from the MODIS land cover products, scaling from flux tower to landscape, and other algorithm limitations [START_REF] Ramoelo | Validation of Global Evapotranspiration Product (MOD16) using Flux Tower Data in the African Savanna, South Africa[END_REF]. Evaluation of MODIS data in previous studies (e.g., [START_REF] Zhang | A continuous satellite-derived global record of land surface evapotranspiration from 1983-2006[END_REF][START_REF] Mu | Improvements to a MODIS global terrestrial evapotranspiration algorithm[END_REF] have shown a good agreement between the data and eddy flux tower observations. The consideration of associated errors to the observation for imbalance problem in water budget closure (using Equation 1) is beyond the scope of this study, and the post-processing is restricted to filtering out the highly noisy measurements.

In-situ Measurements

Groundwater in-situ measurements are used to assess filters' results. To this end, we provide bore stations datasets over the Arkansas-White, Ohio, and Lower Mississippi basins within the Mississippi Basin from USGS and Lachlan, Namoi, and Macquarie-Bogan basins within the Murray-Darling Basin from New South Wales (NSW) Government. The distribution of groundwater in-situ stations is presented in Figure 2. Monthly well measurements are acquired and time series of groundwater storage anomalies are generated. Generally, a specific yield is required to convert well-water levels to variations in groundwater storage (GWS) in terms of equivalent water heights [START_REF] Rodell | Estimating groundwater storage changes in the Mississippi River basin (USA) using GRACE[END_REF]Zaitchik et al., 2008). Following [START_REF] Strassberg | Comparison of seasonal terrestrial water storage variations from GRACE with groundwater level measurements from the High Plains Aquifer (USA)[END_REF], we use an average (0.15) of specific yields ranging from 0.1 to 0.3 (suggested by [START_REF] Gutentag | Geohydrology of the High Plains aquifer in parts of Colorado, Kansas, Nebraska[END_REF] over the Arkansas-White, Ohio, and Lower Mississippi basins, and 0.13 specific yield from the range between 0.115 and 0.2 (suggested by the Australian Bureau of Meteorology (BOM) and [START_REF] Seoane | Regional GRACE-based estimates of water mass variations over Australia: validation and interpretation[END_REF] for the Lachlan, Namoi, and Macquarie-Bogan basins.

Furthermore, we use in-situ soil moisture (SM) measurements obtained from the moisturemonitoring network (http://www.oznet.org.au/), as well as International Soil Moisture Network (https://ismn.geo.tuwien.ac.at/). These data provide long-term records of measured volumetric soil moisture at various soil depths for distributed stations (cf. Figure 2). For each station and each depth, soil moisture anomalies over the study period are calculated. Following [START_REF] Strassberg | Evaluation of Groundwater Storage Monitoring with the GRACE Satellite: Case Study High Plains Aquifer, Central USA[END_REF], data for stations with shallow measurements are upscaled using soil moisture data from deeper stations. We then calculate average soil moisture storage anomalies from all stations within a 1 • ×1 • cell. The same averaging process is done for groundwater measurements. Afterwards, area-weighted anomaly of groundwater and soil moisture are used to achieve GWS+SM. We use these GWS+SM, following [START_REF] Strassberg | Evaluation of Groundwater Storage Monitoring with the GRACE Satellite: Case Study High Plains Aquifer, Central USA[END_REF] and [START_REF] Longuevergne | GRACE Hydrological estimates for small basins: Evaluating processing approaches on the High Plains Aquifer, USA[END_REF], to evaluate the performance of different filters considered in this study. This method does not account for snow water equivalent, canopy, and surface water storages due to their small contribution in TWS over the Mississippi (less than 5%, e.g., [START_REF] Strassberg | Comparison of seasonal terrestrial water storage variations from GRACE with groundwater level measurements from the High Plains Aquifer (USA)[END_REF] and Murray-Darling (less than 6%, e.g., BOM and Burrell et al., 2015) basins. In addition to GWS+SM, we also compare the results with only GWS by computing their correlation coefficients (see details in Section 4.2).

Methods

In this section, first, details of the proposed KeFIn technique are discussed. Afterwards, the other implemented filters including four filters based on the basin averaging approach and two filters that use scale factors' are presented. These techniques are chosen due to their popularity in hydrological studies.

3.1. Kernel Fourier Integration (KeFIn) Filter

The KeFIn Method -First Step

The KeFIn approach follows a straight forward image processing technique, which has been widely applied to geophysical images to enhance their visual interpretation and geological understanding [START_REF] Zhang | Application of image enhancement techniques to potential field data[END_REF]. The application of image enhancement methods is also beneficial for users that are less familiar with processing and filtering the standard GRACE L2

products. The KeFIn includes two processing steps: (1) designing a 2D destriping filter in the spectral domain, and (2), defining an efficient averaging kernel to estimate basin average TWS and at the same time decreasing the leakage-in and -out in the grid domain. A 2-D filter in the spectral domain (Hichri et al., 2012;[START_REF] Zhang | An effective filtering for GRACE time-variable gravity: Fan filter[END_REF] is defined as:

G(u, v) = F (u, v) • H(u, v), (2)
where G(u, v) stands for a Fourier transform of the noisy TWS fields with u and v being spatial frequencies, F denotes a Fourier transform of the ideal (unperturbed) signal (here the 'signal part' or the 'true' TWS values), H is a Fourier transform of a 2-D smoothing kernel to suppress the 'noise' part of the observations, and the dot represents the matrix multiplication. Ideally, F can be estimated by applying an inverse filtering if G and H are known.

In general, however, the information on H does not exist, and its determination usually requires some trial-and-error procedures. Besides, noise in data sets can be amplified leading to the destruction of previous attempts made in reconstructing the TWSs. One solution for restoring F is to use the Wiener Filter (W i ) as F = W i • G, which allows to use an averaging kernel as H to estimate F . Here, a motion filter is used as an averaging kernel (H) to mitigate the south-north stripping problem with different smoothing lengths, which provides us a convolutive filter with different averaging. More detail on creating the kernel with various smoothing lengths can be found e.g., in Bhagat and Gour (2013) (see Equation 5). The impact of smoothing length on the final TWS estimations is presented in Section 4.1.

Thereafter, F can be estimated using H and the Wiener Filter process as:

F (u, v) = |H(u, v)| 2 • G(u, v) |H(u, v)| 2 • H(u, v) + K , ( 3
)
where K is a signal to noise ratio [START_REF] Le Roux | Consistent Wiener filtering: generalized time-frequency masking respecting spectrogram consistency[END_REF]. A suitable estimate for K can be derived as:

K = S G /S F , (4)
where S G is estimated from the power spectral density of the noisy observed signal (G), and S F is derived from the power spectral density of the ideal (unperturbed) signal (F ). The main difference between the new filter and an ordinary Gaussian filter at this stage is the inclusion of the parameter K, which makes Wiener filter more robust and better suited to reduce highfrequency spatial patterns that likely correspond to high magnitude striping patterns. Besides, it introduces a reasonable trade off that minimizes errors of the smoothing process. In order to calculate K in Equation 4, S G is derived from G. For S F , where no information of ideal signal F exists, one can estimate the power spectral density of TWS estimated from a hydrological model and use the mean/median of the estimated powers of S F (see details in [START_REF] Pitas | Digital Image Processing Algorithms[END_REF].

Alternatively one can derive S F by trial-and-error from a range of values (here [0 10]) to control the smoothness of the output, e.g., when the signal is very strong relative to the noise, selecting K ≈ 0 yields less smoothed signals. Different values of K and their impacts on the smoothness of TWS estimations are discussed in Section 4.1. Here, we also use average model TWS estimates from GLDAS NOAH during the study period to compare with the value of K obtained through trial-and-error. The proposed scheme retains most of the high-frequency (spatial) changes that are usually over-smoothed by an ordinary smoothing process [START_REF] Sonka | Image Processing, Analysis, and Machine Vision[END_REF].

The KeFIn Method -Second Step

In the second step of the KeFIn filter, we try to mitigate the problem that arose from the previous stage, i.e., leakage effects caused by spatial smoothing. In what follows, first, spatial averaging and the leakage problem are discussed, then a kernel is defined to reduce the leakage-in and leakage-out errors at the same time. Spatial averaging (Equation 5) is usually applied for improving surface mass anomalies within a specific area [START_REF] Swenson | Methods for inferring regional surface-mass anomalies from Gravity Recovery and Climate Experiment (GRACE) measurements of time-variable gravity[END_REF][START_REF] Longuevergne | GRACE Hydrological estimates for small basins: Evaluating processing approaches on the High Plains Aquifer, USA[END_REF][START_REF] Vishwakarma | Minimizing the effects of filtering on catchment scale GRACE solutions[END_REF],

F R = 1 R F h dΩ, (5)
where,

R = h dΩ, (6)
and F R is the change in vertically integrated water storage averaged over the region of interest,

shown by R, with the integrals done on a sphere. In both equations, h is a basin kernel with values 1 inside the basin and 0 outside of it as,

h(X) =      1 if X ∈ R 0 if X ∈ Ω -R. (7)
X refers to the positions on the surface of the Earth and Ω refers to the entire Earth's surface.

Let us assume that F is derived after applying a filter (that contains smoothing) in step 1. The smoothing moves signals both inside and outside of the basin. In the following, we start by separating the signal F inside and outside the basin and investigate the effects of smoothing leading to F .

The whole water storage changes can be written as a summation of water storage signals inside and outside the basin following [START_REF] Vishwakarma | Minimizing the effects of filtering on catchment scale GRACE solutions[END_REF] represented by the terms F h and F (1 -h), respectively, in Equation 8 as,

F = F h + F (1 -h), = F R + F 1-R . (8)
This is equal to Equation 9 after applying the smoothing procedure from the first step, i.e.,

F = Fl + Fl * , ( 9
)
where Fl is the smoothed signals inside the basin (with leakage out effects) and Fl * refers to the smoothed signals outside the basin (with leakage in effects). By multiplying both sides of Equation 9 by h (Equation 10) and (1 -h) (Equation 11), we achieve the filtered water storage over the region R and outside of it (1 -R).

FR = FlR + E leakage in , (10)
F1-R = Fl * 1-R + E leakage out . (11)
Considering that FlR and Fl1-R are the attenuated signals of F R and F 1-R , [START_REF] Longuevergne | GRACE Hydrological estimates for small basins: Evaluating processing approaches on the High Plains Aquifer, USA[END_REF] showed that they are related using a scaling factor s. For signals inside the basin (the same approach can be used for signals outside the basin), it can be shown that,

F R = s FlR , (12)
s = h dΩ h h dΩ , (13)
with h derived by smoothing h. Equation 10, thus, can be rewritten as,

F R = s ( FR -E leakage in ). ( 14
)
To be able to estimate F R , one needs to calculate the leakage error (E leakage in ) first. To this end, we developed a kernel to account for both leakage in and leakage out errors. The proposed method looks for stronger anomalies outside the basin (for leakage in) and inside the basin (for leakage out). The definition starts by creating a kernel expressed in terms of spherical harmonics as:

υ c lm υ s lm = θ φ ψ(θ, φ) Plm (cos(θ)) cos(mφ) sin(mφ) sin(θ). (15)
In Equation 15, Plm are the normalized associated Legendre functions, υ c lm , υ s lm represent the spherical harmonic coefficients and the summation covers the entire surface of the Earth. The definition of the mask filter ψ is very important and different literatures have found various methods to implement this. For example, [START_REF] Seo | Simulated estimation of hydrological loads from GRACE[END_REF] use a Gaussian filter to smooth mentioned kernel inside a basin (for B 1 and B 2 in their study). [START_REF] Swenson | Estimated accuracies of regional water storage variations inferred from the Gravity Recovery and Climate Experiment (GRACE)[END_REF] applied Lagrange multiplier rather than a Gaussian filter. Here, we use a different definition and instead of simply having a value 1 inside a basin, the method tries to maximize signals concentrated in different regions while decreases their effects on the surrounding signals. For the leakage in effect, ψ contains values outside the basin with special focus on strong anomalies while for the leakage out effect, it considers values inside the basin again with a concentration on strong anomalies. Accordingly, the mask filter ψ is defined through the following procedure.

Note that in the following, we consider F (the smoothed signal from step 1) as a 2D matrix and apply an image processing procedure (as follow) to extract strong signals.

A: The calculated F in the first part of the filtering process is used to create F as a measure of spatial variability of GRACE TWS.

F = ( ( F -min( F )) (max( F ) -min( F ))
). ( 16)

Then, the 2D intensity matrix (I),

I =      1 if F > S b 0 if F < S b , (17)
can be used to identify strong anomalies using the normalised F (given by F ). The B: A high pass filter, e.g., Laplacian filter (Gonzalez andWoods, 1992, 2002) using Equation 18, is applied to intensify strong anomalies (found in [A]) and reduce their effects on surrounding anomalies.

L = 1 sinθ ∂ ∂θ (sinθ ∂I ∂θ ) + 1 sin 2 θ ∂ 2 I ∂φ 2 . ( 18
)
C: Convolving the filtered matrix L with a Gaussian filter (W in Equation 19), which can be applied with different averaging radii. Smoothing is applied because converting the basin kernel from spatial to spectral domain introduces short-wavelength errors due to the Gibbs effect and introduces artificial fluctuations around the high contrast edges [START_REF] Zeng | Partitioned image filtering for reduction of the Gibbs phenomenon[END_REF].

L = W (θ, φ, θ , φ ) L(θ , φ ) dΩ , (19)
In Section 4.1, the impact of the smoothness on the final averaging values is assessed. It should be mentioned here that this step is not restricted to the application of a Gaussian filter, and one can use anisotropic filter such as the DDK smoothing filters proposed by [START_REF] Kusche | Decorrelated GRACE time-variable gravity solutions by GFZ, and their validation using a hydrological model[END_REF]. Nevertheless, in the following we only discuss the application of Gaussian smoothing for the sake of simplicity. The mask filter ψ is then calculated by ψ = 1 + L, which can be used in Equation 15to estimate υ c lm and υ s lm . Figure 3 illustrates a schematic performance of the three steps above.

The final form of the basin kernel (υ) is built as,

υ(θ, φ) = 1 4π ∞ l=0 l m=0 {υ c lm cos(mφ) + υ s lm sin(mφ)}. ( 20
)
The created kernel is multiplied by the smoothed field from the first step to estimate F N using,

F N = F • υ, (21)
where the operator • performs a pixel-wise multiplication. Once F N is computed, it is used rather than F to estimate leakage in and leakage out (Equations 22 and 23). To estimate the leakage in, we only consider F N outside the basin and apply smoothing to capture its effect inside. A similar process can be done to compute the effect of leakage out by only considering anomalies inside the basin. The smoothing in these procedures can be done by applying either the same smoothing procedure as the first step of the proposed filter or using a Gaussian filter, e.g.,

E leakage in = h(θ, φ) 4π W (θ, φ, θ , φ ) (1 -h(θ , φ )) F N (θ , φ ) dΩ , (22)
E leakage out = 1 -h(θ, φ) 4π W (θ, φ, θ , φ ) h(θ , φ ) F N (θ , φ ) dΩ . (23)
The estimated E leakage in is used in Equation 14 to obtain the averaged water storage over the region of interest. The example of the KeFIn filter performance in the second step is presented in Figure 4. Synthetic signals are produced in the spatial domain (Figure 4a) and are smoothed using an ordinary Gaussian filter (Figure 4b). The application of the KeFIn with two different sets of parameters are shown in Figures 4c and4d. The effects of the filter are clearly visible from the reduction of signals interferences caused by leakage. Implementing the filter with various Gaussian filter sizes (r) and different S b (as in Equation 17) yields different results.

Detailed results that indicate the filter's sensitivity to different parameters are presented in Section 4.1. Figure 5 provides a flowchart that summarizes the filter process using the KeFIn algorithm.

Basin Averaging Kernel Methods

Averaging using basin functions or basin kernels is a common approach for estimating basin scale TWS (see e.g., [START_REF] Swenson | Methods for inferring regional surface-mass anomalies from Gravity Recovery and Climate Experiment (GRACE) measurements of time-variable gravity[END_REF]. The kernel h (cf. Equation 7) can be expanded in terms of spherical harmonic coefficients and subsequently combined with L2 potential coefficients to obtain basin averaged GRACE TWS estimates (see e.g., [START_REF] Swenson | Estimated accuracies of regional water storage variations inferred from the Gravity Recovery and Climate Experiment (GRACE)[END_REF], and more details in Section 3.1). Different kernel averaging methods will likely result in different signal attenuation and displaced mass anomalies based on the shape and size of the basins [START_REF] Werth | Evaluation of GRACE filter tools from a hydrological perspective[END_REF]. [START_REF] Swenson | Methods for inferring regional surface-mass anomalies from Gravity Recovery and Climate Experiment (GRACE) measurements of time-variable gravity[END_REF] introduced the spatial averaging kernel for regional studies that try to minimize leakage errors coming from outside into the area of interest by isolating the signals inside the area (see also [START_REF] Swenson | Estimated accuracies of regional water storage variations inferred from the Gravity Recovery and Climate Experiment (GRACE)[END_REF]. Their approach reduces short wavelength effects using a smooth averaging kernel with less power on short wavelengths using Lagrange multiplier rather than applying a Gaussian filter. For the Lagrange Multiplier method, we apply a smoothing radius of 300 km. Furthermore, we use a time-dynamic filter proposed by [START_REF] Seo | Simulated estimation of hydrological loads from GRACE[END_REF]. Here we use filter number three (from four types of their filters), which can be directly applied to GRACE L2 products. This is a dynamic filter that scales spherical harmonic coefficients using the ratio of signal variance and signal plus noise variance that employs a least squares optimum approach. The method is based on the Lagrange Multiplier Method [START_REF] Swenson | Estimated accuracies of regional water storage variations inferred from the Gravity Recovery and Climate Experiment (GRACE)[END_REF] while assuming that the root-mean-square (RMS) of the signal over the target basin is known from GLDAS model (for more details, see [START_REF] Seo | Simulated estimation of hydrological loads from GRACE[END_REF][START_REF] Seo | Terrestrial water mass load changes from Gravity Recovery and Climate Experiment (GRACE)[END_REF]. Here we use GLDAS NOAH for this purpose.

In a different approach, Han and Simons (2008) tried to maximize the ratio of the energy of the function within the target region (h) by constraining regional contributions to global spherical harmonics spectra based on [START_REF] Simons | Localization of the gravity field and the signature of glacial rebound[END_REF]. They argued that the resulted localized coefficients increase the signal-to-noise ratio. This method is also applied in the present study with the spectrum band-limited to spherical harmonic degree and order of 25.We also use a data-driven approach recently introduced by [START_REF] Vishwakarma | Minimizing the effects of filtering on catchment scale GRACE solutions[END_REF], where leakage in and out are separately solved using a catchment mask and a filter kernel. A Gaussian filter of half width radius of 350 km (following [START_REF] Vishwakarma | Minimizing the effects of filtering on catchment scale GRACE solutions[END_REF] is used to suppress the noise before implementing this approach in the present study. The data-driven filter is sensitive to basin sizes in a way that noise increases as the catchment size decrease (see [START_REF] Vishwakarma | Minimizing the effects of filtering on catchment scale GRACE solutions[END_REF], for more details).

Scaling Factor Methods

Landerer and Swenson (2012) suggested the use of a scaling (gain) factor, which can be multiplied with filtered GRACE TWS estimates. In this study, monthly simulations of the GLDAS NOAH are used as synthetic input TWS (a summation of snow water equivalent, canopy water storage, soil layers, and surface water) to estimate scaling factors following [START_REF] Landerer | Accuracy of scaled GRACE terrestrial water storage estimates[END_REF] as in Equation 24, where the goal is to find the scaling factor α by minimising the quadratic sum of difference M between original ( S T ) and filtered ( S F ) GLDAS TWS fields, i.e.,

M = ( S T -α S F ) 2 . ( 24
)
Following [START_REF] Landerer | Accuracy of scaled GRACE terrestrial water storage estimates[END_REF] and [START_REF] Long | Deriving scaling factors using a global hydrological model to restore GRACE total water storage changes for China's Yangtze River Basin[END_REF], synthetic TWS data is converted to spherical harmonics and truncated at degree and order 60. We then apply the destriping procedure after [START_REF] Swenson | Post-processing removal of correlated errors in GRACE data[END_REF] and a 300 km Gaussian filter to smooth high-degree and order noises. The model-derived TWS estimates before ( S T ) and after ( S F ) filtering are used to calculate scaling factors. In this study, two methods of scaling factors at grid points and basin scale are computed and used for comparison with the newly developed KeFIn and other filtering techniques. All filters used in this study are presented in Table 1.

Application Example of the Proposed KeFIn Filter

First, the performance of the KeFIn filter with respect to both leakage-in and leakageout errors is assessed, for which two tests are performed that correspond to each type of error (leakage-in and leakage-out). Setup (i), the signal is only introduced inside a basin and GRACElike TWS noise is added as described in Section 2.2. A 300 km half width radius Gaussian filter (Jekeli, 1981) is then applied to smooth the introduced signals and noises, which causes signal leakage outside the basin. Setup (ii), TWS signals are introduced only outside a basin to assess the leakage-in effects. The KeFIn filter is applied to post process both scenarios as shown in Figure 6. In Figure 6a bottom, the blue line represents the introduced synthetic TWS while the green lines show the signal after the application of a Gaussian filter. In Figure 6a, the results correspond to a cross section at 3 • S that passes the Amazon basin, South America, and in Figure 6b, they correspond to a cross section at 41 • N crossing the Huang He Basin, China.

The results clearly indicate that the Gaussian filter attenuates the original signal and causes leakage-out and leakage-in effects shown in Figures 6a and6b, respectively. The smoothed signals of the KeFIn filter are shown by the red lines, which in both cases better follow the initial TWS (blue lines). It is worth mentioning that if there was no striping noise added to the signal, the red curve (KeFIn) would have closely reproduced the true signal (blue curve). (RMSE) time series of the filters performances using the synthetic data over the basins is the KeFIn filter with other methods (F 1 to F 6 in Table 1). It is clearly visible in Figure 7 that the KeFIn filter works properly in both basins. RMSE values over the Colorado Basin (Figure 7c) suggest that the application of the KeFIn filter (i) successfully decreases leakage error, and

TWS (cm) TWS (cm)

Results

In Section 4.1, various filtering techniques (cf. Table 1) are tested on the synthetic TWS data while in Section 4.2, the results from filtering the real GRACE data are assessed against direct observations of water fluxes through the water balance equation (Equation 1), as well as in-situ groundwater measurements.

Filter Results Based on Synthetic Data

There are two effective factors in each step of the proposed KeFIn filter, which potentially change the final filtering outcomes. The main aim here is to find out which choice yields an optimum performance of the filter in terms of leakage error reduction. Figure 8a contains the results of applying the first step of KeFIn while considering different sizes for the motion filter (controlling the smoothing of north-south stripping error) and K to mitigate the signal attenuation. Each scenario (using Equations 3 and 4) is applied separately to the basins and the average errors for all basins and are represented in Figure 8a. From our investigations, using K from GLDAS provides the best results with ∼14.76% higher leakage error reduction with different filter lengths. Considering K as a constant can lead to a promising result with the value of 1 with 58 mm average error. On the other hand, motion filters with bigger windows better decrease errors, where the optimum value in this study is derived from the 75 degree motion filter size. As mentioned, the first part of the filter deals with colored/correlated noise of high-frequency mass variations (i.e., stripes). In order to investigate the performance of this step of the filter, we compare its results with the widely used destriping algorithm by [START_REF] Swenson | Post-processing removal of correlated errors in GRACE data[END_REF] and DDK smoothing filter following [START_REF] Kusche | Approximate decorrelation and nonisotropic smoothing of time variable GRACEtype gravity field models[END_REF] and [START_REF] Kusche | Decorrelated GRACE time-variable gravity solutions by GFZ, and their validation using a hydrological model[END_REF].

We apply these filters over all basins and illustrate the average results in Table 2. Note that we apply the KeFIn method with best cases of K and motion filter for the comparison presented in Table 2. Based on these results, the first step of the KeFIn filter performs comparable to other filters in terms of RMSE reduction. The level of RMSE reduction, as well as correlation improvements for the KeFIn filter are larger in most of the cases, particularly compared to Gaussian with 250 km radii and DDK3. In addition, we used the same experiment this time for the second part of the filter (cf.

Equations 6 and 8) while applying diverse values of S b and selecting various smoothing radii (half-width radius, r) for the Gaussian filter. Using the best cases of K and motion filter length, we analyze the effects of different S b and r on errors (Figure 8b). In general, results indicate that a higher S b needs lower r to derive better results. Nevertheless, applying the second part of the KeFIn filter with S b = 0.5 and r = 300 km performed better in most of the cases. For comparison, all the filters of F 1 , F 2 , F 3 , F 4 , F 5 , and F 6 (cf. Table 1) as well as the KeFIn filter are then applied on the GRACE-like synthetic TWS fields. A summary of these results is presented in Table 3. For every basin, we estimate F R (averaged signals inside the basin) and F 1-R (averaged signals outside the basin) using each filter and compare the results to initial unperturbed TWS values inside and outside the basins by calculating the RMSE and correlation coefficients. Note that for a better assessment, seasonal variations are removed from time series. The average results for the study period, i.e. 2002-2013, and for all the 43 basins (cf. Figure 1) is given in Table 3. Note that detailed RMSE values for each individual basin can be found in the Appendix. From Table 3, it can be seen that higher correlations, both inside and outside the basin, can be found by applying the KeFIn filter. Estimated measures indicate that the KeFIn filter is more successful in recovering the spatial distribution of the synthetic TWS estimates. Overall, the KeFIn filter performs better both inside and outside the basins with an average of 73.6% TWS recovery from the perturbed synthetic data (cf. Table 3).

⠀戀⤀

Our results further indicate that the KeFIn filter works well over smaller river basins such as the Colorado, Ohio, Lachlan, and the Namoi basins, showing maximum ∼81% TWS recovery from noisy data. We also found that in 35 out of the 43 basins, the proposed filter provides the lowest RMSE (cf. Appendix). Nevertheless, in the other 8 cases, the KeFIn approach still demonstrates a promising performance in terms of RMSE reduction. Overall, Table 3 suggests that the proposed filter performs better in more than 80% of the basins. We further assess the performance of the filters, using independent data sets such as water fluxes. Therefore, TWS changes are evaluated through the water balance equation (cf. Equation 1) using TRMM 3B43-v7 precipitation, AVHRR data to account for evaporation products, and in-situ discharge data over the Amazon, Mekong, Arkansas-White (basins 1 and 31 in Figure 1, respectively), Ohio, Lachlan, Namoi, Lower Mississippi, and Macquarie-Bogan basins (cf. Figure 2). We further assess the results of the different filters against groundwater measurements as mentioned in Section 2.4. TWS estimates after implementing each filter and a summation of groundwater storage (GWS) and soil moisture contents (GWS+SM) are compared in the following basins: Arkansas-White, Ohio, Lachlan, Namoi, Lower Mississippi, and Macquarie-Bogan (cf. Figure 2), where access to in-situ data is provided. For each basin and each filtering method, basin averaged values are compared with GWS+SM. For this purpose, absolute differences between the filtered results and in-situ measurements are illustrated in Figure 11.

(F R -FR ) / C R (F 1-R -F1-R ) / C 1-

䬀攀䘀䤀渀

Similar to the previous section, the minimum errors are found after using the KeFIn filter for these basins. It can be seen from the distribution of error points that the KeFIn results obtain errors with less magnitudes and variances. This indicates the smaller deviations of these results compared to in-situ measurements. Among the other filters, in general, smaller errors are found for F 2 and F 6 . F 2 and F 5 depict less errors over the Ohio Basin and Lachlan Basin. In summary, the KeFIn filter and F 2 better decrease errors over these basins, respectively 38% and 22% (on average) better than the other filters. These show the higher capability of the two filters for reducing errors within smaller basins. For a better comparison, the average errors in Figure 11 for all the basins are shown in Figure 12. Figure 12 illustrates that the proposed KeFIn filter in all the cases has the minimum error (24.13 mm on average). Similar to the two basins discussed earlier in this section, using F 2 , F 3 and to a lesser degree F 3 lead to a higher agreement with observations compared to the other methods (except the KeFIn filter). The results of these filters are much closer to those of the proposed filter in Arkansas-White and Macquarie-Bogan Basins. F 4 seems to have an approximately constant effect on different basins (37.58 mm on average) except for the Ohio Basin. The summary of comparisons between different filtered TWS and in-situ groundwater time series measurements are presented in Table 5. This is performed to show each filter's performance independent against direct observations without incorporating model estimates.

Higher correlations are reported between the KeFIn filter results and in-situ data, which indicates 19.31%, 6.67%, 10.57%, 8.41%, 18.52%, and 6.33% improvements in comparison to F 1 , F 2 , F 3 , F 4 , F 5 , and F 6 , respectively. F 3 , F 6 , and F 4 results are also in good agreement with the [START_REF] Vishwakarma | Minimizing the effects of filtering on catchment scale GRACE solutions[END_REF]. This approach is, however, found to be sensitive to the basin size in a way that noise increases when the catchment size decreases.

Over smaller basins (e.g., Lachlan and Namoi basins), F 2 works significantly better than F 3 and F 4 . This confirms the findings of Han and Simons (2008) that this filter is designed to address leakage errors over basins with a small area (cf. Table 4). In summary, our results

indicate that the KeFIn filter and F 2 are likely better suited to deal with the leakage in small river basins. F 3 , designed by [START_REF] Seo | Simulated estimation of hydrological loads from GRACE[END_REF] and tested over the Mississippi Basin,

shows reliable results over this basin with fewer similar performance in other basins. This likely indicates that filters must be extensively tested over different basins that are of different shapes and sizes with different magnitude and distribution of TWS signals.

We find that the proposed KeFIn filter reduces the leakage errors over ∼82% of the basins with an area less than 1 million km 2 , thus, we conclude it is suitable for leakage error reduction over basins with various sizes and shapes. Comparison with water flux observations indicates that in addition to the KeFIn filter, the recently developed F 6 and F 4 that use a hydrological model to recover GRACE smoothed signals (on a gridded basis), better approximate the derivatives of TWS changes than the other filters. Over the larger basins (e.g., Amazon and Mekong basins), the results of the F 1 and F 5 filters are found to be better than those in the smaller basins. Overall, more consistent leakage reduction within different basins is achieved by the KeFIn filter, F 2 , F 6 , and F 4 considering the results of Figure 10 and Figure 12, as well as Table 3.

Conclusion

In this study, a new GRACE post-processing technique, the so-called KeFIn filtering method, is proposed and its performance in reducing GRACE TWS errors in higher spatial frequencies as well as leakage (in/out) errors is investigated. The KeFIn filtering method successfully mitigates the existing problems with other leakage filtering methods, e.g., the high sensitivity of them to prior models in the scale factor approaches. To demonstrate the benefit of using the KeFln filtering method, two different test scenarios are considered over the 43 river basins of different shapes and sizes. First, all the filtering methods are compared using generated synthetic data with properties similar to real GRACE TWS data within the 43 globally distributed river basins. In addition, we assessed the performance of the filters against water storage changes from water fluxes observation, as well as a summation of observed groundwater storage and soil moisture content over the selected basins. The results show that the KeFIn filter successfully (i) mitigates the amplitude damping caused by smoothing, and (ii) increases flexibility towards a variety of basins (shapes and sizes of basins as well as the magnitude of TWS). It is worth mentioning here that we do not claim that the KeFIn method is able to reduce all possible artificial features appearing in the two steps of the post processing algorithm.

Therefore, further investigations will be done to optimize parameters that are used to define the shape of the KeFIn filter.

Acknowledgement

We would like to thank Dr. Tim R. McVicar, the associate editor of Remote Sensing of Environment, and the three anonymous reviewers for their useful comments, which we used to improve this study. The authors acknowledge the efforts of providers of GLDAS, TRMM, and MODIS products, as well as GRACE ITSG-Grace2014 gravity field models. We would also like to thank USGS, GRDC, ANA, HYBAM and NSW for sharing ground-based data sets. M. Khaki is grateful for the research grant of Curtin International Postgraduate Research Scholarships (CIPRS)/ORD Scholarship provided by Curtin University (Australia). This work is a TIGeR publication.

0 1 Figure 1 :

 11 Figure 1: Shapes, sizes and locations of the world's 43 major river basins (red borders and green areas) used in this study.

Figure 2 :

 2 Figure 2: Distribution of groundwater (red crosses) and soil moisture (cyan circles) stations over the six selected river basins of Arkansas-White, Lower Mississippi, Ohio, Macquarie-Bogan, Namoi, and Lachlan basins.

  threshold S b in Equation 17 is chosen to be a value within [0 1]. Often the median of F can be a good choice for S b . A smaller S b yields a smoother intensity matrix that controls the mass anomalies being considered in the averaging, and which is less weighted. Different values of S b are tested in this study and their results are reported in Section 4.1.

Figure 3 :

 3 Figure 3: A schematic view of the steps for preparing ψ in [A]-[C] described above. (a) shows the initial unperturbed signal, (b) represents the smoothed signal from the first step of the filter (applied with the motion length of 60), (c) is I in step [A] using S b = 0.5, and (d) depicts the kernel ψ created by r = 300 km.

Figure 4 :

 4 Figure 4: Performance of the second step of the KeFIn filter based on synthetic data. (a) Initial TWS anomalies, (b) smoothed TWS using a Gaussian filter with the half-width radius of 500 km. (c) and (d) represent the performance of the KeFIn filter with different factors of S b and r (half width radius in kilometre). In this figure, we show how the KeFIn filter tries to reproduce the signals in (a) based on the smoothed signal (b), which result in (c) and (d).

Figure 5 :

 5 Figure 5: Flowchart of the proposed KeFIn filtering process.

Figure 6 :

 6 Figure 6: Assessing the performance of two filtering techniques on synthetic GRACE-like TWS examples with realistic noise. (a) TWS is introduced in the Amazon River Basin, South America, and (b) TWS is introduced outside of the Huang He River Basin, China. The line plots indicate the TWS after application of Gaussian filter with 300 km radii (green) and the KeFIn filter (red), estimated using the motion length of 60, S b = 0.5, and r = 300 km. Note that the line plot of kernel (black) is also shown in these figures, which are shifted for better visual demonstration. The initial synthetic TWS is represented by the blue lines. Units are in cm.

(

  ii) improved results in relation to other filters. We find approximately 34% RMSE reduction compared to the unperturbed signals by implementing the KeFIn filter. By comparing RMSE values in the Congo basin, again, smaller errors are found for those associated with the KeFIn filter compared to the other six filters applied in this study. This indicates that the KeFIn filter successfully decreased leakage effects based on the GRACE-like artificial data, especially over smaller basins.

Figure 7 :

 7 Figure 7: The KeFIn filter operation over the Colorado (a) and Congo (b) basins using synthetic GRACElike TWS signals and noise. In column (1), the unperturbed water storages are shown; in column (2), the corresponding perturbed water storages are shown, and the results of the KeFIn filtered TWS estimates are presented in column (3). Panel (c) shows the average RMSE results within both basins for the filters listed in

Figure 8 :

 8 Figure 8: Average error (mm) derived after applying the KeFIn filter with different values of K and the motion filter length (a) for the first step of the filter as well as different scenarios that contain S b and r for the second step of the filter (b). (a) indicates that the filter length of larger than 30 km and K between 0 to 2 yield smaller errors, while (b) indicates S b of 0.5 and r = 300 km yield the smallest errors.

Figure 9 :

 9 Figure 9: Comparison between the derivative of filtered TWS (red) and TWF from observations (blue) within the Namoi Basin. Each sub-figure corresponds to one filter and also contains error bars that is computed as the absolute value of difference between GRACE derivatives and the observed TWF.

Figure 10 :

 10 Figure 10: The temporal average of errors defined as derivative of filtered TWS minus observed TWF. Each error bar is estimated after applying the F1 to F6 and KeFIn filters over 8 selected river basins (units are mm).

Figure 11 :

 11 Figure 11: Errors estimated at each epoch after applying the assessed filters F1-F6 and KeFIn on the Ohio (basin number 35, with blue circles) and Lachlan (basin number 41, with red triangles) basins. These values are calculated as differences between in-situ measurements (GWS+SM) and filtered TWS before (E1) and after (E2) removing seasonal effects. The average absolute error is indicated in each sub-figure and for each basin.

Figure 12 :

 12 Figure 12: Average differences between GRACE TWS and observed groundwater plus soil moisture content within 6 river basins. GRACE data are processed using 7 filtering techniques (F1 to F6 and KeFIn filters) over 6 selected river basins (units are mm).

Table 1 :

 1 A summary of the implemented GRACE leakage filtering methods, which are used in this study for comparison with the proposed KeFIn filter.

	Study		Method	Case Study	Evaluation Method	Abbreviation *
	Swenson	and	Lagrange multiplier method	Mississippi River Basin	Using synthetic GRACE data	F1
	Wahr (2002)
	Han and Simons	Localization of Global Geopotential Fields	Java/Sunda trench	Using seismic model based	F2
	(2008)				data
	Seo and Wilson	B1, B2, B3, and B4	Amazon, Mississippi, Lena,	Using synthetic GRACE data	F3
	(2005)			Huang He and Oranje Basins
	Landerer	and	Gridded gain factor	46 globally distributed basins	GLDAS data	F4
	Swenson (2012)
	Landerer	and	Single gain factor	46 globally distributed basins	GLDAS data	F5
	Swenson (2012)
	Vishwakarma et	Data-driven approach	32 globally distributed basins	Closed-loop environment us-	F6
	al. (2016)				ing monthly GLDAS fields
	The present study	Kernel Fourier Integration (KeFIn)	43 globally distributed basins	Using synthetic data and soil	KeF In
					moisture + groundwater data

*

In the last column, the abbreviations are assign to the filters we use in the present study.

Table 1 .

 1

Table 2 :

 2 Average statistics derived after applying different filtering methods over the world's 43 major river basins using synthetic data (after removing seasonal effects) in comparison with the unperturbed synthetic data (F0). Note that the first step of the KeFIn filter is used in this table.

		Gauss (250 km) Gauss (350 km) Gauss (500 km) DDK1 DDK2 DDK3 KeFIn
	RMSE (mm)	78.54	54.13	60.91	57.87	53.19	62.67	52.73
	Correlation	0.73	0.81	0.78	0.83	0.80	0.76	0.81

Table 3 :

 3 Average statistics derived after applying different filtering methods over the world's 43 major river basins using synthetic data in comparison with the unperturbed synthetic data (F0). Averaged signals inside and outside of the basins are calculated using CR = F0 h dΩ and C1-R = F0 (1 -h) dΩ, respectively.

	Method	Inside the Basin	Outside the Basin	TWS improvement (%)
		Correlation RMSE (mm)	Correlation RMSE (mm)

Table 4 :

 4 Correlations between the TWFs as precipitation minus evaporation minus runoff, and the derivatives of TWS changes from each applied filter. The correlation coefficients have been computed at the 95% confidence level.

	Basin	F 1	F 2	F 3	F 4	F 5	F 6 KeF In
	Amazon	0.92 0.93 0.94 0.92 0.91 0.95	0.95
	Mekong	0.85 0.92 0.88 0.88 0.89 0.91	0.93
	Arkansas-White	0.78 0.82 0.81 0.75 0.73 0.75	0.88
	Ohio	0.76 0.82 0.74 0.82 0.81 0.78	0.85
	Lachlan	0.80 0.86 0.82 0.73 0.75 0.84	0.89
	Namoi	0.72 0.87 0.78 0.80 0.82 0.81	0.91
	Lower Mississippi	0.77 0.78 0.79 0.81 0.80 0.78	0.84
	Macquarie-Bogan	0.79 0.85 0.81 0.78 0.74 0.69	0.92
	4.2.2. Comparisons with Groundwater and Soil Moisture

Table 5 :

 5 Correlations between the filtered results and in-situ measured groundwater time series. We find here that those filters based on the averaging kernel, especially F 2(Han and Simons, 2008) and F 3[START_REF] Seo | Simulated estimation of hydrological loads from GRACE[END_REF], deal better with leakage errors over smaller basins compared to those based on scaling factor (F 4 and F 5 ;[START_REF] Landerer | Accuracy of scaled GRACE terrestrial water storage estimates[END_REF]. Nevertheless, in general, F 2 , F 6[START_REF] Vishwakarma | Minimizing the effects of filtering on catchment scale GRACE solutions[END_REF], F 3 , and F 4 perform better than F 1 and F 5 in most of the cases. The grid-based F 4 is found to better reduce leakage errors in comparison to the single gain factor F 5 . Between basin average kernel methods, in general, F 6 and F 2 perform better compared to F 1 . The results confirm that F 6 reduces leakage errors better than other basin average techniques when it is applied over larger basins as mentioned in

	Basin	F 1	F 2	F 3	F 4	F 5	F 6 KeF In
	Arkansas-White	0.78 0.76 0.73 0.69 0.63 0.75	0.81
	Ohio	0.76 0.82 0.73 0.63 0.69 0.84	0.85
	Lachlan	0.69 0.71 0.75 0.67 0.68 0.78	0.83
	Namoi	0.59 0.74 0.64 0.75 0.58 0.66	0.80
	Lower Mississippi	0.54 0.78 0.73 0.66 0.67 0.72	0.78
	Macquarie-Bogan	0.77 0.85 0.81 0.76 0.73 0.82	0.88
	5. Discussion
	Evaluation of the proposed KeFIn filter against common techniques (cf. Table 1) using

different datasets suggests that this filter successfully removes striping and reduces leakage errors over basins of different shapes and sizes. Other filters show a different level of improvements within the world's major 43 basins

(Figures 9, 11, and 12)

.

Table A1 :

 A1 Summary of RMSE (mm) estimated using the unperturbed basin averaged synthetic TWS and the perturbed TWS after using different filtering methods over the 43 river basins. Note that the basins are sorted according to their area.

	Basin	Area (million km 2 )	F1	F2	F3	F4	F5	F6	KeF In
	1 (Amazon)	6.97	31.25 31.83	31.19	30.98 30.88 31.06	30.69
	2 (Ob)	4.40	26.64 24.25	28.77	27.64 28.55 22.93	23.79
	3 (Yenisey)	4.09	29.76 23.17	26.10	21.94 21.11 19.44	17.63
	4 (Lena)	3.99	31.94 36.15 27.56 32.99 33.70 30.71	29.57
	5 (Congo)	3.81	25.24 23.52	24.96	24.19 25.60 21.59	20.47
	6 (Mackenzie)	2.88	24.60 29.23	31.42	26.63 28.56 23.86	22.18
	7 (Parana)	2.64	37.97 31.83	38.18	30.27 32.71 26.97	26.68
	8 (Nile)	2.48	34.17 34.01	33.86	34.15 33.45 37.36	32.79
	9 (Mississippi)	2.35	42.93 38.52	37.51	38.22 43.37 39.83	37.20
	10 (Niger)	2.11	34.56 33.01	29.34	34.93 33.38 27.84	27.78
	11 (Amur)	1.85	52.03 49.35 46.19 50.33 50.24 47.73	48.52
	12 (Yangtze)	1.81	39.90 36.93	38.56	40.20 41.81 36.68	35.75
	13 (Yukon)	1.58	37.91 40.27	38.63	38.10 39.67 37.14	36.69
	14 (Nelson)	1.43	31.50 24.22	30.43	22.99 26.21 24.12	21.41
	15 (Volga)	1.38	30.00 34.84	31.33	32.01 33.71 28.22	28.93
	16 (St. Lawrence)	1.27	39.53 39.97	33.82	36.02 37.68 34.14	32.65
	17 (Lake Eyre)	1.12	24.10 26.49	24.94	17.51 29.60 19.45	16.45
	18 (Zambezi)	1.12	29.54 28.10	31.67	34.98 33.10 29.05	27.67
	19 (Murray Darling)	1.01	46.66 41.51	38.89	40.94 43.42 38.72	37.84
	20 (Danube)	0.93	36.67 35.97	37.39	39.31 41.72 31.40	29.20
	21 (Ganges, Brahmaputra)	0.92	28.92 17.77 33.88	26.03 25.19 29.75	28.25
	22 (Indus)	0.91	41.31 33.39	36.57 32.04 34.04 35.50	33.69
	23 (Orange)	0.90	18.71 14.96	21.82	13.01 16.08 11.67	7.94
	24 (North West Coast)	0.80	16.85 17.97	18.81	22.10 17.58 19.39	16.39
	25 (Huang He)	0.78	33.86 28.70	30.77 23.09 27.11 24.98	23.30
	26 (Sumatra)	0.76	32.46 27.08	28.43	34.61 34.38 28.03	26.19
	27 (Euphrates and Tigris)	0.74	35.91 22.20	19.75	22.00 24.36 19.08	17.53
	28 (Orinoco)	0.73	42.57 35.90	34.65	38.94 35.99 32.69	32.42
	29 (Tocantins)	0.71	25.32 15.02 16.05	18.99 20.29 22.72	20.65
	30 (Ayeyarwady)	0.69	34.75 38.21	36.09	34.17 35.64 35.74	33.97
	31 (Mekong)	0.68	34.46 35.78	32.27	33.06 36.82 38.03	31.78
	32 (Kalahari Stampriet)	0.67	36.01 34.75	37.51 32.32 39.25 35.16	34.29
	33 (Dnieper)	0.65	24.25 23.18	29.80	25.14 25.59 24.09	22.84
	34 (Colorado)	0.63	23.12 19.32	23.05	22.95 23.87 21.79	18.05
	35 (Ohio)	0.52	23.31 22.54	25.06	21.82 26.96 23.38	20.46
	36 (Sirdaryo)	0.51	32.56 27.98	30.35	25.47 25.63 29.40	24.74
	37 (Central East Coast)	0.49	40.21 42.51	37.21	38.64 39.23 41.20	36.31
	38 (Western Mediterranean)	0.45	31.64 28.42	28.59	32.44 36.54 37.91	27.06
	39 (Namoi)	0.43	21.80 12.80	14.33	17.09 25.31 19.53	12.43
	40 (Kamchatka)	0.40	43.06 33.69	33.62	40.80 37.73 38.89	32.90
	41 (Lachlan)	0.08	34.05 32.42	28.11	25.79 32.05 28.41	24.46
	42 (Yalu)	0.03	13.07 12.33	11.40	12.76 18.40 14.73	8.82
	43 (Lower Mississippi)	0.01	23.94 18.79	20.24	25.13 28.42 21.13	18.49

Han, S.-C., Simons, F.J., 2008. Spatiospectral localization of global geopotential fields from the Gravity Recovery and Climate Experiment (GRACE) reveals the coseismic gravity change owing to the 2004 Sumatra-Andaman earthquake, J. Geophys. Res., 113, B01405, http://dx.doi.org/10.1029/2007JB004927. Harig, C., Simons, F. J., 2015. Accelerated West Antarctic ice mass loss continues to outpace East Antarctic gains. Earth and Planetary Science Letters, 415, 134141. http://dx.doi.org/10.1016/j.epsl.2015.01.029. Henck, A.H., Montgomery, D.R., Huntington, K.W., Liang, C., 2010. Monsoon control of effective discharge, Yunnan and Tibet, Geology, v. 38, p. 975-978. Hichri, S., Benzarti, F., Amiri, H., 2012. Robust Noise Filtering in Image Sequences. International Journal of Computer Applications (0975 8887), V 50 No.18. Huffman, G.J., Adler, R.F., Bolvin, D.T., Gu, G., Nelkin, E.J., Bowman, K. P., Hong, Y., Stocker, E.F., Wolff, D.B., 2007. The TRMM multi-satellite precipitation analysis: quasiglobal, multi-year, combined-sensor precipitation estimates at fine scale. Journal. Hydrometeorol., 8 (1), pp. 3855. Jekeli, C., 1981. Alternative methods to smooth the Earth's gravity field. Technical report 327, Department of Geodesy and Science and Surveying, Ohio State University. Khaki, M., Hoteit, I., Kuhn, M., Awange, J., Forootan, E., van Dijk, A.I.J.M., Schumacher, M., Pattiaratchi, C., 2017a. Assessing sequential data assimilation techniques for integrating GRACE data into a hydrological model, Advances in Water Resources, Volume 107, Pages 301-316, ISSN 0309-1708, http://dx.doi.org/10.1016/j.advwatres.2017.07.001. Khaki, M., Schumacher, M., J., Forootan, Kuhn, M., Awange, E., van Dijk, A.I.J.M., 2017b. Accounting for Spatial Correlation Errors in the Assimilation of GRACE into Hydrological Models through localization, Advances in Water Resources, Available online 1 August 2017, ISSN 0309-1708, https://doi.org/10.1016/j.advwatres.2017.07.024. Klees, R., Revtova, E. A., Gunter, B. C., Ditmar, P., Oudman, E., Winsemius, H. C., 2008.

The design of an optimal filter for monthly GRACE gravity models. Geophysical Journal International, 175, 2, 417-432. http://dx.doi.org/10.1111International, 175, 2, 417-432. http://dx.doi.org/10. /j.1365International, 175, 2, 417-432. http://dx.doi.org/10. -246X.2008.03922.x. .03922.x.

To this end, we calculate TWF (from Equation 1) over each basin (see Section 2.3). Figure 9, for example, shows the results of this comparison within the Namoi Basin. The figure also contains error bars for every filter representing the differences between the observed TWF and those derived by estimating the temporal derivative of filtered TWS change. It can be seen that the results of the KeFIn filter are much closer to the observed TWF with the smallest average error of 11.13 mm and overall 13% higher correlation in comparison with the other filters.

Average error estimates within different basins corresponding to each filter are illustrated in Figure 10. Errors after applying the KeFIn filter are found to be the smallest in all the assessed basins. We find F 2 , F 4 , and to a lesser degree F 6 to be efficient in most of the cases, especially over the Ohio Basin. More details on results can be found in Table 4, in which correlations between the TWFs (estimated as precipitation minus evaporation minus runoff) and the derivatives of TWS changes that are filtered by all implemented filtering methods are represented. Maximum correlations are calculated for the proposed filter with 0.89 average correlation. A higher correlation is achieved from all the filters over the Amazon and Mekong basins, which can be due to their stronger signals compared to other basins. Results from F 2 , F 3 , and F 6 are found to have larger correlations to TWFs than those from F 1 and F 5 .

Appendix

The following table shows the basin averaged RMSE values calculated by each filtering technique. The results in the table are temporally averaged (between 2002 and 2013), and indicate that the KeFIn filtering method works better compared to other filters in 35 out of the 43 basins, especially over smaller basins.