An Iterative ICA-Based Reconstruction Method to Produce Consistent Time-Variable Total Water Storage Fields Using GRACE and Swarm Satellite Data

Observing global terrestrial water storage changes (TWSCs) from (inter-)seasonal to (multi-)decade time-scales is very important to understand the Earth as a system under natural and anthropogenic climate change. The primary goal of the Gravity Recovery And Climate Experiment (GRACE) satellite mission (2002–2017) and its follow-on mission (GRACE-FO, 2018–onward) is to provide time-variable gravity fields, which can be converted to TWSCs with ∼ 300 km spatial resolution; however, the one year data gap between GRACE and GRACE-FO represents a critical discontinuity, which cannot be replaced by alternative data or model with the same quality. To fill this gap, we applied time-variable gravity fields (2013–onward) from the Swarm Earth explorer mission with low spatial resolution of ∼ 1500 km. A novel iterative reconstruction approach was formulated based on the independent component analysis (ICA) that combines the GRACE and Swarm fields. The reconstructed TWSC fields of 2003–2018 were compared with a commonly applied reconstruction technique and GRACE-FO TWSC fields, whose results indicate a considerable noise reduction and long-term consistency improvement of the iterative ICA reconstruction technique. They were applied to evaluate trends and seasonal mass changes (of 2003–2018) within the world’s 33 largest river basins.


Introduction
The Gravity Recovery And Climate Experiment (GRACE, 2002(GRACE, -2017 satellite mission [? ] was a joint United States (National Aeronautics and Space Administration, NASA) and German (Deutsche Zentrum für Luft-und Raumfahrt, DLR) space mission, which provided estimates of variations in the gravity field arising from mass movements within the Earth system. GRACE was launched on 17th March 2002, and consisted of two identical spacecraft, one following the other, in the same orbit separated by roughly 220 km. The variations in the distance between the two GRACE satellites were measured by a K/Ka-Band Ranging (KBR) system, which can largely be attributed to global mass changes within the Earth's interior, at its surface, and in the atmosphere.
Over any typical 30-day span (the nominal data accumulation interval of GRACE gravity field solutions) of the non-repeating orbit configurations, a dense ground track grid can be obtained (see e.g., [? ]); therefore, usually (almost) one month of the pre-processed GRACE along-track KBR data were used to compute time-variable gravity field solutions that are known as the GRACE level 2 data (GRACE L2, [? ]). These L2 data were used to provide satellite-based estimates of global and regional terrestrial water storage change (TWSC, a vertical summation of water storage changes) with a few hundred kilometers spatial resolution, e.g., over Asia (e.g., [? ? ? ]), Africa (e.g., [? ? ]), and Australia (e.g., [? ? ? ]). On a global scale, TWSCs are discussed, e.g., in [? ? ? ? ? ]. All these studies came to the same conclusion that GRACE L2 products are suitable for studying large scale TWSCs on decadal, annual, and inter-annual time scales. A review of GRACE applications in climate studies is provided in [? ].
Initially, GRACE was targeted to cover a five-year period, which was exceeded by ten years till October 2017. The GRACE Follow-On (GRACE-FO, [? ]) mission was launched in May 2018, but suffered a failure of the main instrument processing unit between July and October 2018. This has led to approximately one year of data gap. Many studies have been working to bridge the gap between GRACE and GRACE-FO missions using alternative observations. However, there is no single satellite mission or proxy observation which is able to fill this gap with comparable quality (see e.g., [? ? ? ]).
Geodetic satellite laser ranging (SLR) is another geodetic remote sensing observation technique that provides information about temporal variations of the Earth's gravity field at the lowest degrees of the spherical harmonic spectrum (e.g., up to degree and order 6), and covers a longer time series than GRACE mission (see e.g., [? ? ? ? ? ]). The long time span of collecting data by SLR missions provides an opportunity to determine the Earth's gravity field before GRACE and during the gap between GRACE and GRACE-FO. Talpe et al. [? ] combined high-resolution gravity fields derived from GRACE with the low-resolution fields from SLR/DORIS (Doppler Orbitography and Radiopositioning Integrated by Satellite, DORIS) tracking data, to reconstruct global time-variable gravity fields for the period of 1993-2013. In [? ], GRACE observations were combined with SLR data to reconstruct a consistent long-term estimation of the Earth's gravity field.
In addition to GRACE and SLR missions, satellites in low earth orbit (LEO), namely the European Space Agency (ESA)'s Swarm Earth explorer mission, provides Earth gravity field models from alternative sources of gravimetric data. Swarm consists of three identical satellites (launched on 22 November 2013) [? ], where two satellites are orbiting the Earth side-by-side in nearly polar orbits currently (November 2019) at the altitude of ∼ 430 km (Swarm-Alpha and Swarm-Charlie, separated by 15 seconds in time at their equatorial crossing) and one is flying at higher altitude of ∼ 500 km (Swarm-Bravo). Swarm's main mission objective is to measure the Earth's magnetic field; however, since the Swarm satellites are equipped with geodetic-quality GPS receivers, accelerometers, star-tracker assemblies, and laser retro-reflectors, their high-low satellite-to-satellite tracking (hl-SST) observing system provides the opportunity of computing time-variable gravity field solutions as shown by [? ? ], and [? ]. Recently, [? ] provided a long-term time series of monthly gravity field solutions as a combination of SLR data, GPS and K/Ka-band observations of the GRACE mission, and GPS observations of the three Swarm satellites. In their study the monthly gravity field derived from Swarm were used to fill the gap between GRACE and GRACE-FO missions. GRACE L2 data are provided up to degree and order 96 [? ]. Since the fields are noisy, isotropic (e.g., Gaussian in [? ]), or anisotropic and de-striping or decorrelation filters (e.g., [? ? ? ? ], with ∼ 300-900 km half radius based on the application) have often been applied to achieve TWSC fields of few hundred km resolution. Swarm's monthly gravity fields, in contrast, have been computed up to degree and order 70, but contain significant time-variable signal only up to about degree 13 (∼1500 km or longer resolution), and the fields are relatively noisy [? ], where filters with the smoothing radius between 500-833 km were used to estimate TWSC fields with few thousand km resolution [? ]. This fact limits the application of Swarm data for, e.g., continental hydrology, where higher resolution estimations are desirable [? ]; therefore, improving the spatial resolution and reducing the level of noise in Swarm fields are required to enhance their application for hydrological applications. This topic is addressed in this study.
In order to fill the gap between GRACE and GRACE-FO, in a recent study, [? ] developed an empirical approach to predict GRACE-like TWSC fields using the climate inputs as indicators. Their study addresses the prediction of TWSC in large river basins but it is restricted to the regions that contain high correlation between climate variations and water storage changes. The empirical technique that they applied also requires a learning phase that might be relatively sensitive to the episodic variations caused by natural (e.g., floods and droughts) and anthropogenic factors (e.g., water storage over extractions).
In this study, two reconstruction formulations (Approach 1 and Approach 2) are proposed, which can be used to merge GRACE and Swarm time-variable gravity fields. The proposed approach of this study entirely relies on gravity field outputs, which can be derived from different satellite missions. For both formulations, one can apply statistical decomposition techniques such as principal component analysis (PCA, [? ]) and independent component analysis (ICA, [? ? ]) or their extensions (Complex PCA and Complex ICA, as in [? ]). Here, we prefer ICA over PCA, since it performs better to separate TWSC signal from noise (see, e.g., [? ? ]). Besides, ICA's independent modes can better be related to physical phenomena as shown by e.g., [? ? ? ]. Approach 1 follows a common data reconstruction procedure similar to [? ] and [? ], but in the spatial grid domain (TWSC), although Approach 2 is novel and formulated as an iterative ICA-based reconstruction, using the whole period of available GRACE and Swarm TWSC time series, covering 2003-2018.
Validating GRACE TWSC fields was performed by comparisons with geophysically relevant signals from independent data sources. For example, hydrological models were used to assess GRACE data over land area (e.g., [? ? ]); however, many studies, e.g., [? ? ], indicate that the trend and seasonality of these models are uncertain and such comparisons might be performed with care. Other studies, e.g., [? ? ? ] suggested the use of in-situ data such as ocean bottom pressure and global navigation satellite systems (GNSS)-derived vertical loads to assess GRACE signals. These validations are more suited to evaluate low degree spectral (long wavelength spatial) variations. Over (large) lakes, satellite altimetry and in-situ data can be used to assess GRACE data; however, a reliable reduction of hydrological signals surrounding these lakes is a challenging task (see, e.g., [? ? ? ]). In this study, satellite altimetry data over selected lakes are used as an indicator of changes in trend and/or seasonality of TWSCs derived from the reconstructed data. In addition, TWSCs derived from GRACE-FO is used to evaluate the reconstructed TWSC in this study, for some selected months after the year 2018.
In what follows, data sources of this study are described in Section ??, and the two alternative reconstruction techniques are formulated in Section ??. Our goal was to illustrate that the optimal combination of the GRACE and Swarm (three-satellite constellation) temporal gravity field solutions improves the resulting retrieved TWSC fields compared with mass estimation from only one type of satellite mission pleonastic. In Section ??, the main results of this study are presented including a comparison between the two approaches used for reconstruction, an evaluation of the reconstructed TWSCs in comparison with that of the GRACE-FO product, as well as an evaluation of mass changes (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) within the world's 33 major river basins. Finally in Section ??, the main conclusions and an outlook are provided. This paper has two appendices that provide detailed information about the efficiency of the presented reconstruction approach in decreasing the level of noise and spatial leakage in GRACE and Swarm data.

Data
In this section, satellite geodetic data to be used to estimate global TWSC fields and their post-processing are described.

GRACE L2 Fields
The GRACE mission consists of two identical spacecraft (GRACE-A and GRACE-B) following each other in tandem in a single orbital plane. The twins are equipped with a K/Ka-band ranging system (KBR) to provide inter-satellite ranges. In addition, both satellites carry GPS receivers that are used for orbit determination, as well as accelerometers to remove non-gravitational forces prior to gravity modeling. The time-variable GRACE gravity field solutions, known as GRACE level 2 (L2) products, are provided by a number of institutions, using different processing approaches, background models, and assumptions.
For this study, the latest release of monthly GRACE L2 product (RL06, [? ]), covering April 2002 to September 2016, and the previous release (RL05), covering November 2016 to January 2017, expressed as spherical harmonic coefficients up to degree and order 96, were downloaded from the Center for Space Research (CSR, http://www2.csr.utexas.edu/grace/). The biases between RL05 and RL06 data will be eliminated by implementing the reconstruction procedure introduced in Section ??. To evaluate the reconstruction results of this, GRACE-FO data of 2019 were downloaded from the CSR website. It is worth noting here that, due to the aging batteries on the GRACE satellites, no ranging data were collected during certain orbit periods, and hence no gravity fields could be computed. These gaps occur approximately every 5-6 months, mostly after 2011, and last for 4-5 weeks. The months with no GRACE data are June and July 2002, June 2003, January and June 2011, May and October 2012, March, August and September 2013, February, July, and December 2014, June, October, and November 2015 October in 2016.

Swarm L2 Gravity Field Models
The three satellites of ESA's Earth's Magnetic Field and Environment Explorer (Swarm) mission, which were launched in November 2013, were designed specifically to accurately measure the Earth's magnetic and electric fields [? ]. Global positioning satellite (GPS) receiver, one aboard each Swarm satellite, allows them to produce precise orbit solutions accurate to centimeters in conjunction with post processing of the orbital data using gravity field models [? ].
The Swarm level 1b products (processed raw satellite measurements) are released by ESA to represent the best estimate of the actual values of the electric and magnetic field at each point in space and time. The Swarm Satellite Constellation Application and Research Facility (SCARF) [? ], in the form of a consortium of several research institutions, use level 1b data and auxiliary sources to provide Swarm level 2 data [? ]. The background force models for the Swarm processing are defined corresponding to GRACE models (e.g., as in case of GRACE, the same de-aliasing products to reduce short-term variations in the atmosphere and ocean masses).
A number of institutes have produced gravity field models from Swarm L2 data, with different approaches and different kinematic orbits. Four official centers that provide Swarm gravity field are: Astronomical  In this study, the combined monthly Swarm L2 gravity model provided by (ASU, [? ]) are used, which can be downloaded from http://www.asu.cas.cz/~bezdek/vyzkum/geopotencial/ in terms of spherical harmonic coefficients up to degree and order 40, covering the period of December 2013 to December 2018. These gravity field models were derived from the GPS data of three satellites (Swarm-Alpha, Swarm-Charlie, and Swarm-Bravo) through the acceleration approach [? ], and the kinematic orbits for gravity inversion in this model were provided by IFG, Graz (ftp://ftp.tugraz.at/ outgoing/ITSG/tvgogo/orbits/). Though these fields might contain a higher level of noise, compared to the official ESA's Swarm L2 data, they are purely based on one processing principal, which makes the choice consistent with that of GRACE (in Section ??). Changing the Swarm solutions does not considerably change the final reconstruction results (corresponding evaluations are not shown here).

Treatment of the Degree One and Two Coefficients
Degree one coefficients or reference frame motion is not observed by the GRACE or GRACE-FO low-low K/Ka-Band ranging measurements. In this study, we augmented them by those derived from [? ], which will have impact on the seasonality and to a much lesser degree on the trends computed from TWSC over the land.
It is also known that the estimation of degree 2 coefficients of GRACE L2 is unreliable using KBR ranging only, which is primarily due to the measurement strategy (i.e., only along-track satellite-to-satellite ranging is available for gravity inversion). Here, we replaced these coefficients by the ones derived from SLR data following [? ]. The degree one and two coefficients of the Swarm fields were also replaced to be consistent with the treatment we used for GRACE L2 products.

Data Post-Processing
The higher degrees of GRACE-derived spherical harmonic coefficients are strongly affected by correlated noise. Thus, a smoothing operator needed to be applied to suppress the effect of this noise in the GRACE-derived TWSC maps. In this study, an isotropic Gaussian filter [? ] with the radius of 500 km was applied to filter TWSC derived from GRACE L2 products. By this choice, we provided end-users the opportunity to implement their own treatment of error correlations. It is worth mentioning here that we did not exclude the GRACE fields with worse quality than normal, e.g., those affected by the GRACE's known battery problem. These fields can be replaced within the reconstruction procedure. The same post processing, as for GRACE, is implemented to convert GRACE-FO data of 2019 to TWSC fields.
To reduce the high magnitude noise in the Swarm fields, a smoothing filter needs also to be applied. Da Encarnação et al. [? ] suggested the smoothing filters with half-width radius of between 500 to 833 km are effective to reduce the noise of degree 12 to 20 coefficients. They noted that this choice results in a higher correlation (than other options) between Swarm and GRACE estimates. In this study, the same Gaussian filter (that was used for GRACE) but with 1000 km half-wavelength radius was applied to smooth Swarm fields. It is worth mentioning that, since Swarm solutions will be replaced by the iterative reconstruction values, one could also select the radius of 500 km to smooth the Swarm L2 solutions. The choice of smaller radius leads to a higher noise in Swarm data but it does not introduce a considerable impact on the final reconstruction estimates.

Glacial Isostatic Adjustment (GIA) Corrections
Besides water storage variations, time-variable gravity fields also contain signals associated with the glacial isostatic adjustment (GIA), which reflects the Earth's visco-elastic response to the deglaciation of the last ice age sheets. GIA manifests itself as a trend at the decadal or longer time-scales. For hydrological applications, the GIA effect is normally removed as a linear trend based on the output of GIA models, for which the GIA forward model of [? ] is used in this study. It is worth mentioning here that this correction is needed for the hydrological interpretation of mass estimations from time-variable gravity fields and does not have considerable impact on the reconstruction results themselves. GIA corrections are applied after implementing the reconstruction that merged GRACE and Swarm gravity field solutions.

Conversion to the TWSC Fields and Reducing the Leakage
Formulation of Wahr, et al. [? ] is used to convert the time-variable gravity fields to TWSC. Grid spacing was selected to be global 1 • × 1 • similar to that of GRACE Tellus official data provider (https://grace.jpl.nasa.gov/data/get-data/). Our TWSC fields contain the signal over both land and the oceans. To compute basin averages within the world's 33 largest river basin the leakage reduction and averaging approach in [? ] was used that simultaneously minimizes the summation of leakage-in and leakage-out contributions.

Method
Independent component analysis (ICA, [? ], chapter 4) was applied here to relate GRACE and Swarm time-variable gravity fields. ICA is an extension of the principal component analysis (PCA, [? ? ], chapter 3) technique, which is often used for statistical reconstruction. We first describe the formulation of PCA and ICA, then two reconstruction approaches are introduced in Sections ?? and ??, whereas both can be applied using either PCA or ICA methods; however, we prefer ICA (on PCA) reconstruction as is described later in this section.
PCA or its equivalent singular vector decomposition (SVD) technique. It is a second-order decomposition method, because it uses eigenvalue decomposition of the auto-covariance matrix built on multivariate samples to provide orthogonal eigenvectors that optimally explain the total variance [? ]. Let us consider that the TWSC estimates after removing their temporal mean are stored in the data matrix X. The dimension of X is t × s, where t represents the number of time epochs and s stands for the number of grid points (for 1 • × 1 • fields from either of missions, one has 64,800 grid points). If one treats GRACE and Swarm data individually, X contains 154 fields of GRACE (2003GRACE ( -2017 and 61 of Swarm fields (2013)(2014)(2015)(2016)(2017)(2018). PCA or SVD is defined by expanding X as: whereĒ contains the spatially orthogonal vectors (empirical orthogonal functions, EOFs), associated with temporally uncorrelated time series, known as principal components (PCs), and stored in the columns of matrix P t×m . Thus, the columns ofĒ s×m correspond to the unit-norm eigenvectors e 1 , e 2 , ..., e m of the auto-covariance matrix C X = E{X T X} (E(.) being the expectation function), i.e., E TĒ = I. Columns of P and E orĒ are arranged with respect to the magnitude of the corresponding singular-values λ 1 , λ 2 , ..., λ t , in which λ 1 > λ 2 > ... > λ t , and t is the maximum number of singular-values corresponding to the auto-covariance matrix C X , i.e., the number of time epochs.
In the case of SVD in Equation (??),Ē remains unchanged, butP contains normalized PCs in its columns, i.e.,P TP = I, and Λ is diagonal and holds the singular values λ ordered according to their magnitude (P =PΛ). If needed, scaled eigenvectors can be derived as E =ĒΛ. The total variance of X can be derived by summing the square of singular values, i.e., ∑ λ i 2 , i = 1, ..., t; therefore, the portion of the total variance presented by the first m modes can be computed as: The data matrix X can be re-estimated using the first m PCA components as: Independent Component Analysis (ICA) technique. ICA is a higher (than second) order statistical decomposition technique, which incorporates more statistical information (than PCA) of available data samples in its computational procedure to derive empirical components that are statistically independent. The ICA technique that was applied here follows the formulation of [? ?
], in which EOFs and PCs of Equation (??) are rotated using an orthogonal rotation matrix R in an optimization procedure to estimate components that are as statistically independent as possible. The statistical decomposition of the data matrix X using ICA can be formulated as: where PR and ER are rotated components and always only one of them is independent. For the temporal ICA technique, columns of PR, and for the spatial ICA method, columns of ER are statistically as independent as possible (details of the temporal and spatial ICA techniques and the approach to estimate R can be found in [? ], Chapter 4 and Section ?? of this paper). Our experiments indicate that both temporal and spatial ICA techniques provide similar reconstruction results; therefore, in the following we do not distinguish between these two alternatives and we simply call the method ICA. The advantages of applying ICA (instead of PCA) for reconstruction include: (1) ICA can better separate signal from noise, especially when the distribution of data sets are not Gaussian (this is the case for TWSC time series as shown by [? ? ? ]); (2) Theoretically, independent components can better extract the physically separated phenomena such as trend, seasonality (see, e.g., [? ]) and variations of TWSCs due to large-scale teleconnection indices (see e.g.,[? ? ]); therefore, relying on the ICA modes offers more representative basis to combine various types of data; (3) One usually needs a smaller number of ICA modes (than the number of PCA modes) to reconstruct a certain portion of variance in the original data [? ], which makes the computation of reconstruction more stable. The first m ICA components can be used to re-estimate the data matrix X as:

Joint Diagonalization (JD)
Different criteria exist that can be used to measure independence and equivalently derive the ICA transformation (see examples in [? ? ? ]). Here, following [? ], the fourth-order cumulant, which is a tensorial quantity used to estimate R. A cumulant tensor Q(M) (e.g., [? ]) is a four-dimensional tensor with an arbitrary n × n matrix M = (m ij ) (with zeros everywhere except 1 at index (i, j)), or where x is a random vector of n × 1. The associated cumulant matrix Q(M) is defined as: with Σ being the auto-covariance matrix of x. The desired rotation matrix R that drives ICA in Equation (??), also minimizes the off-diagonal entries of the cumulant tensor Q (consisting the matrices Q 1 to Q m ) in Equation (??) can be estimated as: where the operator diag(.) produces a matrix with all off-diagonal elements being set to zero, and . represents the Frobenius norm or the Hilbert-Schmidt norm. The estimated rotation matrixR in Equation (??) is defined as a product of k rotation matricesV m (R = ∏ k m=1V m ), where eachV m diagonalizes a cumulant matrix (Q m ). Cardoso and Souloumiac [? ] reported that since the cumulants matrices in Equation (??) are estimated based on length-limited samples, they only approximately share the same eigenstructure and do not fully commute. As a result, the exact solution of the joint diagonalization in Equation (??) does not perfectly diagonalize the fourth-order cumulant matrix of Equation (??); therefore, the independent components that are defined following this computation are statistically 'as independent as possible' rather than 'absolutely' independent.

Data Reconstruction-Approach 1
The first reconstruction approach (shown as 'Approach 1' in the results section (Section ??)) follows the commonly used statistical decomposition-based reconstruction, which often is implemented using the PCA approach to link two geophysical fields, for example see [? ? ? ? ].
Here, we followed the methodology by [? ], who applied a reconstruction to combine GRACE time-variable fields with low degree and order time-variable gravity fields from SLR. Their implementation considers PCA to relate GRACE and SLR fields in the spectral domain. The spatial patterns of GRACE were then fitted to SLR fields to extend the fields to the past two decades; however, here, we replaced PCA with ICA and applied the reconstruction in the grid domain to combine GRACE and Swarm fields. Our motivation to select ICA was described before, but selecting the grid domain for reconstruction was due to the fact that the computation of statistics (or empirical functions) using 64,800 grid entries was numerically more stable than using a much smaller number of coefficients in the spectral domain. Another reason was that scaling the TWSC fields can be done by the area of these grids or the cosine of their mid-latitude in the spatial domain, which was more straightforward than in the spectral domain with no unique best way to account for the influence of coefficients of different degree and order.
Approach 1 of this study consists of three steps, i.e., Step-1: ICA (Equation (??)) is applied on the GRACE TWSC fields X GRACE ; Step-2: Spatial patterns of GRACE from the previous step are fitted to those of Swarm TWSC (X Swarm ) while considering the errors (described in Section ??); Step-3: The new temporal values are used to reconstruct TWSC fields and replace those of Swarm. Here is the summary of the reconstruction algorithm, which is labeled by 'Approach 1' in the Results Section: X GRACE : n 1 = 154 × 64, 800 X Swarm : n 2 = 61 × 64, 800 Step-1: X GRACE = PRR TĒT = P E T . Step-2: Swarm with m being defined by setting v = 95% in Equation (??), M represents errors of Swarm TWSC fields, and Λ contains the first m eigenvalues of GRACE TWSC fields.
Step-3: Replace the Swarm's TWSC fields by the reconstruction of Equation (??), and usingˆ P from Step 3 andĒ from GRACE data. Our numerical experiments (not shown here) indicate the selected 95% threshold results in reducing noise and the higher cut-off values yield TWSC fields with higher level of noise.

Iterative Reconstruction-Approach 2
For the second reconstruction approach ('Approach 2' in the Results Section, (Section ??)), we took advantage of the useful property of ICA (one can replace ICA with PCA), i.e., its ability to extract dominant variance of spatio-temporal data sets through its independent modes. Here, we assumed that both smoothed GRACE-and Swarm-derived TWSC fields (after removing their temporal mean) were stored in the matrix X, making it 215 × 64, 800. The rows of X were sorted with respect to the time at which the GRACE and Swarm fields were available; therefore, here, we did not distinguish between the fields derived from GRACE and Swarm, though we knew that the spectral resolution of Swarm gravity fields was smaller than that of GRACE and their noise level was higher.
The proposed approach in this section consists of three steps, i.e., Step-1: ICA (Equation (??)) is applied on the data matrix X to derive the first rough estimation of the variability represented by the time-variable gravity fields of these two missions; Step-2: The dominant modes of the previous step are used in Equation (??) to calculate the reconstruction ofX but replace the values of the Swarm fields only; Step-3: The new field is inserted into Step-1 and Step-2 to calculate new ICA modes and re-evaluate the process, which can be repeated until convergence, i.e., until ICA modes do not considerably change from one iteration to the next. This iterative procedure is very advantageous, because it uses all the information in GRACE and Swarm fields at once, i.e., the data matrix X consists of TWSC values from both missions. It also harmonizes their TWSC fields, because after a certain number of iterations the auto-covariance matrix will be stable and that is why the ICA modes do not change considerably after a certain iteration. We will show (in Section ??) that the reconstruction results of the Approach 2 contain less noise than that of Approach 1 and therefore they might be more useful for hydrological applications.
However, one important remaining point is the decision on the number of relevant ICA modes to be retained to reconstruct the signal in Step-3. Here, a simple variance cut-off rule was used that means only those dominant modes that correspond to 95% of the total variance (Equation (??)) in the data matrix were retained. Here is the summary of the reconstruction algorithm, which is marked as 'Approach 2' in the Results Section: X GRACE : n 1 = 154 × 64, 800 X Swarm : n 2 = 61 × 64, 800 Initial: X = sort([X GRACE ; X SWARM ], time) where sort(., time) is an operator that sorts any matrix according to the ascending values of time in year. Rows of X are weighted by the inverse of covariance matrices in terms of TWSC.
Step-1: X = PRR TĒT . Step-2: Step-3: Find the values ofX that correspond to the Swarm fields and replace the initial values of X and start from Step-1. GRACE and/or GRACE-FO data that are missing or those with considerably high level of noise can also be replaced in this step. This iteration stops when the ICA of Step-1 is not found to be significantly different from previous iteration, i.e., the scaled root-mean-squared-errors between modes is less than 10 −6 .
It is worth mentioning that we have observed the general evolution of the ICA reconstruction through different iterations. A summary of our main results is presented in what follows. By starting the initial values of TWSC from original GRACE and Swarm data, the variance of the dominant ICA modes will generally increase, whereas the amplitude of the less dominant modes will decrease. This is coherent with the interpretation that the initial guess overestimates high frequency signals and underestimates large scale signals by disregarding any local information around the low resolution Swarm fields. This means that the damped amplitude of TWSC signals derived from Swarm because of its lower resolution (compared to GRACE) causes high frequency localized anomalies in the long time series. Keeping this in mind, by looking at how the ICA structure is changed in each iteration, we can see that the modification in each iteration is to diminish the differences between GRACE and Swarm estimations. This is visualized here in Figure ?? by showing the third and fourth independent components that contain both long-term and seasonal frequencies. Starting at the end of 2013, by introducing the Swarm fields, considerable noise can be detected in the ICs, which is eliminated by iterations. It is worth mentioning that IC3 and IC4 are shown here because they are not dominated by trends and their changes through iterations can be better visualized. Our assessments also indicate that the magnitude of noise in other ICs have also been reduced by iterations (results are not shown). For example, IC1 and IC2 represent a long-term linear trend and annual cycles, respectively, for which the original estimations were considerably noisy but after iterations the magnitude of the noise is considerably reduced. By such reduction in the magnitude of noise the consistency between TWSC estimations from GRACE and Swarm increases ([? ], see also the numerical results of Section ??). This emphasizes that the iterative reconstruction approach optimizes the combination of gravity field solutions with different level of accuracy; therefore, this method can effectively replace contemporary approaches to estimate optimal relative data weights between models and combine them (numerical results are presented in Appendix 1 and 2).

Estimation of Reconstructions Errors
In order to estimate uncertainties of the reconstruction from Approach 1 and Approach 2, we followed the Bootstrap approach used similar to that in [? ]. First, n realizations of the data matrix X i , i : 1; ...; n were generated by considering the original data matrix X and adding errors from GRACE and Swarm TWSC fields. The error fields were generated randomly but by considering the TWSC covariances. Applying the three steps of Approach 1 and 2 on these realizations yields n sets of reconstruction results, from which the uncertainty of independent modes was estimated. In the numerical implementation of this study, we selected n to be 1000. These error fields were used to estimated error budget of basin average estimations in Section ??.

Results
In this section, first we assessed the reconstruction results from Approach 1 and 2 (Section ??) described in the previous section. We show that the iterative reconstruction yields TWSC fields that contain a lower level of noise, and therefore, in Section ??, these fields were used to assess water mass changes within the world's 33 largest river basins. All GRACE fields were smoothed by the Gaussian filter with 500 km half-width radius. Our numerical evaluations indicate that the application of alternative filters such as the DDK anisotropic filters in [? ] does not alter the general results of reconstruction, but it might slightly change the basin averaged mass estimations of Section ??.

Evaluating the Reconstruction Results of Approach 1 and 2
First, we compared the results of reconstructions for two arbitrary selected epochs, for which both GRACE and Swarm fields are available. Figure ??A1,B1 indicate GRACE fields in two selected dates in the year 2014 and 2015, where Figure ??A2,B2 show Swarm derived mass estimation for similar dates. By comparing A1 and B1 with A2 and B2, respectively, one can detect similar large-scale mass anomalies, for example over Amazon, Greenland, and North America. The positive-negative anomalies of the northern and southern hemispheres are also well reflected in the Swarm TWSC fields (compare the fields on the first column on left with the second column from left); however, even by a simple visual comparison of the Swarm and GRACE fields, it is clear that mass changes from Swarm data are over-smoothed, i.e., overall a damping of ∼ 1.5 to 2 in terms of standard deviation (StD) was estimated to match Swarm TWSC fields over land to that of GRACE. In general, the anomalies seen over land and the oceans cannot be easily interpreted since they are unrealistically extended over the land area, such as those anomalies over the Australian continent in A2 or those detected over South Asia and Africa in B2, which do not match the results of GRACE in A1 and B1, respectively. This is partly due to the fact that the smoothing radius for denoising GRACE data was 500 km, and 1000 km for Swarm; however, even, GRACE TWSC fields that are smoothed by a filter of 1000 km radius indicate more details than those of Swarm's (results are not shown here). This happens because the GRACE hl-SST observations are more sensitive to mass redistribution justifying the fact that GRACE data are computed up to higher degree and order potential coefficients (than Swarm).
The results of reconstruction from Approach 1 and 2 are shown in Figure ??A3,B3 and Figure  ??A4,B4, respectively. These fields can be compared directly to those of GRACE in A1 and B1. In general, we find less magnitude of noise in the TWSC fields of Approach 2 (iterative reconstruction) than those of Approach 1. We show this by computing standard deviation (StD) of differences between the reconstructed values and original GRACE data over land and oceans, where the mean values are reported in Table ??. This was theoretically expected that the iterations in Approach 2 reduces the noise part (any spatially/temporally random behavior) of GRACE and Swarm data. The reason is that the dominant eigenvectors are oriented to retrieve the 'signal' part of data and the dominant part of 'noise' is discarded (by eigenvalue decomposition and the joint diagonalization of the forth-order cumulants). During the iterations, the orientation of ICs are modified so that their evolution becomes less random. As a result, Approach 2 yields spatially and temporally homogenized reconstruction fields as described in the end of Section ??.  Table 1. Mean of the standard deviation of differences between the reconstructed values and the original GRACE TWSC-terrestrial water storage changes data over the 33 major river basins (marked as land) and rest of the globe (marked as ocean).

2014.7082
Ocean ( Figure ??A1,A2, and the reconstruction results of Approach 1 are shown in Figure ??B1,B2. Those of Approach 2 are shown in Figure ??C1,C2. By comparing the extension of anomalies derived from original Swarm data and reconstructions over land, we see that both reconstruction techniques localize the mass anomalies and their patterns are closer to what we expect from GRACE mission. We also notice that the results of Approach 2 better keep the extension of original signal in Swarm. For example, by a visual comparison of Figure ??A1,B1,C1, one can find higher similarity between the dipole anomaly over South America derived from Figure ??A1,C1, compared to that of Figure ??B1, where the extension of anomaly in the northern part of South America is considerably reduced. Another example is the TWSC fields extended over the middle and west parts of Africa, which are detected in Figure ??A1,C1 (amplitude of ∼ 150 mm) but the amplitude considerably decreases in Figure ??B1 (amplitude of ∼ 50 mm). Considering the selected results in 2018 (Figure ??A2,B2,C2), mass anomalies in the Antarctic are found to be considerably underestimated in the reconstruction results of Approach 1 (an average amplitude of ∼ 185 mm from Swarm (Figure ??A2) and Approach 2 (Figure ??C2) versus ∼ 110 mm from Approach 1). To make sure that the reconstruction approach works properly, when the GRACE-FO data are available, we extended the time of investigation to end of 2019. Figure ?? compares Swarm ( Figure  ??A1,A2), the reconstruction of Approach 2 (Figure ??B1,B2), and GRACE-FO (Figure ??C1,C2) that correspond to 2019.2890 (March 2019) and 2019.6246 (July 2019), respectively. Here, we only show the hydrological signals (i.e., oceans ans polar regions are masked). A visual comparison indicates that the reconstructed fields (Figure ??B1,B2) are very close to those of GRACE-FO (in Figure ??C1,C2), compared to those of Swarm (Figure ??A1,A2). An average correlation coefficient of 0.89 are found between the reconstructions and GRACE-FO during 2019. An evaluation of the noise level and leakage error for smoothed GRACE and Swarm data as well as the reconstructed fields is presented in Appendix 2. The results show that the reconstruction considerably improves the accuracy of TWSC estimations, i.e., the level of noise, is reduced by 400% on average and the median of the leakage errors for all the world's 33 largest river basins is reduced by 51%.
To assess the error behavior of our reconstruction results, Figure ?? shows the differences during 2019.3726 (April 2019) and 2019.7082 (August 2019), where the magnitude of differences over the land area is found to be less than 10 mm (which is acceptable for smoothed TWSC fields). This justifies the suitability of these fields in hydrological studies.
In the light of our results in the previous section, we believe that the reconstruction results from Approach 2 have better quality, i.e., contain less noise and the mass anomalies are well localized. Therefore, these fields are used to assess the global mass redistribution of 2003-2018 in terms of trends, seasonal and sub-seasonal cycles, as well as those related to teleconnection evens, see similar studies, e.g., [? ? ? ? ? ]. To illustrate the benefit of the iterations in Approach 2 to reduce the noise in TWSC fields, a comparison between the results of the first iteration with those of 250th iteration in terms of basin-averaged TWSC is shown in Figure A1. An illustration of TWSC time series on two selected regions in the USA (latitude: 40.5 and longitude: −85.5) and South Africa (latitude: −9.5 and longitude: 24.5) is shown in Figure A2. Linear trend, acceleration corresponding to inter-annual or longer climate episodes (with respect to the time 2010.5) along with annual and semi-annual cycles were extracted from basin-averaged TWSC time series within the 33 basins using a linear model (l(t) = a t + b t 2 + c cos(2πt) + d sin(2πt) + e cos(4πt) + f sin(4πt) + (t), where l(t) contains the basin averaged TWSC; t is the time in year minus 2010.5; a, b, c, d, e are coefficients to be estimated using a least squares fit; (t) contains the residuals). The significance of these coefficients, e.g., a for the linear trend and b for non-linear trend (or known as acceleration) was evaluated using a Fisher test as in [? ], and the results are summarized in Table ??. In the following, we report some features of the basin averaged TWSC that are relevant to make sure the reconstructed TWSC fields are reliable to be used for large-scale hydrological assessments. The results are reported according to different continents. North America: Six river basins are shown in Figure ??, for which statistically significant linear trends are found within the Yukon, Mackenzie, Mississippi, and St. Lawrence River Basins. The non-linear trend (simply called acceleration in what follows) component is not found to be statistically significant in any of the six basins in North America. Positive linear trend is found in the St. Lawrence and Mississippi River Basins that are consistent with independent measurements. Here, we considered altimetry measurements (from https://ipad.fas.usda.gov/cropexplorer/global_reservoir/) to test whether surface water storage (as suggested by [? ]) is a dominant contributor to TWSC trends within the St. Lawrence River Basin containing the Great Lakes of Superior, Huron, Michigan, Erie, and Ontario. Our results indicate that the level of these five lakes increased during 2003-2018 as shown in Figure ??. It seems that the increase in surface water causes the positive trends derived from our reconstructed fields. It is worth mentioning here that the use of altimetry data in this study was to provide an evidence of changes in trends and seasonality rather than providing a quantitative assessment of these components, because for the latter, one needs to implement a sophisticated leakage reduction (e.g., [? ]), which will be performed in future.
Linear trends in four other river basins, Yukon, Nelson, Mackenzie, and Colorado are found to be negative, which can be associated to both irrigation and global warming. For example, Lena and Mackenzie are two glacier-dominant basins in the Arctic region influenced by the climate warming over the past four decades, and consequently, caused decreasing in water storage in these regions [? ? ? ]. Scanlon et al. [? ] indicated that Colorado River Basin experienced multiyear intense droughts between 1980 to 2015, and storage depletion in this basin is dominated by reductions in surface reservoir and soil moisture storage, as well as reductions in groundwater storage due to imbalances between water supplies and demands during the climate extremes.  South America: Four basins are assessed in South America, with no significant trends in the three basins of Amazon, Orinoco, and Tocantins (except Parana). The magnitude of seasonal TWSC considerably decreased in all the basins, which can be attributed to the drought events, for which, e.g., [? ? ] indicated that tropical and sub-tropical South America has experienced several extreme hydrological droughts in the last decades, e.g., 2005, 2010, and 2016 droughts, due to higher than usual sea surface temperatures in the Tropical Pacific (including El Niño events) and Atlantic Oscillations; however, the hydrological consequences due to the drought of 2016 are reported to be more severe than previous ones (see e.g., [? ? ? ]).
The increasing linear trend in water storage within the Parana River Basin, despite of no significant rainfall increases, can be explained by increasing the mean annual discharge during the last 2 decades by converting forests to agricultural land (see, e.g., [? ? ]). The small decrease in the magnitude of seasonal TWSC within the Parana River Basin during 2017-2018 can be attributed to the extreme drought of 2016 (see the trend values in Table ??). Africa: Six river basins are assessed in the African continent (see Figure ??). The linear trends in the Niger, Nile, Okavango, and Orange River Basins, and the acceleration in Congo are found to be statistically significant. Figure ?? shows that the level of three major lakes within the Nile River Basin from altimetry measurements. For all three lakes, a positive trend is found during the period of 2003-2018, which is correlated with our finding from the reconstructed TWSC fields within this basin.
The Niger River Basin is the third longest river in Africa and located across several east-west oriented climate zones, which causes a complex composition of geographical factors. These factors, as well as the impact of climate changes [? ? ] can be the reason of intense pressure on water resources in this region [? ]. Many studies such as [? ? ? ] indicated that because of the changes in land cover and the expansion of agriculture in Nigerian Sahel, discharge regime has been changes and as a result the groundwater storage increased in this basin, which is reflected in the reconstructed TWSC estimates.
The positive trends in TWSCs of Okavango can be associated to the heavy rainfalls, seasonal floods, and geographical location of Okavango Delta, which cause the flow of water into this Delta. The time period of the seasonal floods in Okavango Delta are varied in different years; however, it generally occurs between March and May every year just after the local summer rains (between November and March) in the region of the Delta [? ? ]. Although the positive linear trends in TWSC is found in Africa, decreasing in the magnitude of seasonal TWSC in 2017-2018 is caused by the low net-fluxes after the extreme drought of 2015/2016 (see, e.g., [? ]). This should be noted here that a separate interpretation of changes in precipitation or evaporation on TWSC should be performed with care. For example, [? ] indicate how changes in moisture transport impact hydrological cycle changes.
As a result, not necessarily, an increase or decrease of fluxes would change the seasonal amplitude of the TWSC in the same direction.  Europe, West Asia, and East Asia: Among eight major river basins of the Europe and West and East Asia, statistically significant linear TWSC trends can only be observed within the Euphrates (with negative trend) and Amur (with positive trend) River Basins. Voss et al. [? ] and Forootan et al. [? ] reported a significant decline in groundwater, soil moisture, and surface water storage within the Euphrates River Basin between 2003-2013, which was associated to both irrigation and long-term net-precipitation deficit. Increasing trend in TWSC of Amur River Basin can be due to the increases in deforestation of coniferous forests, which contributed to the increases in rivers' flow and flood events in this basin (see, e.g., [? ]). A negative linear trend can be seen in reconstructed TWSCs within the Lena River Basin, which is intensified during 2014-2018. Previous studies, e.g., [? ? ], argue that warmer temperatures due to climate change have impacted the hydrological regime of rivers within this basin via changes in precipitation type (rain replacing snow). They also hypothesize that other factors, such as the melting of permafrost, glaciers, and/or changes in groundwater conditions to have contributed to this trend; therefore, our results confirm their independent investigations. It also provides an evidence of the contribution of water storage decline in the detected change in hydrological regime. Figure ?? indicates that the reconstructed TWSCs within the Ob and Yenisei River Basins do not contain considerable linear trends, and small trend values are detected within Dnieper, Danube, and Caspian-Volga as part of the European Russia. Grigoriev et al. [? ] analysed TWSC in the European part of Russia and found a non-significant reduction in the TWSC between 2002-2015, which was caused by a decline in surface and groundwater water storage. A negative linear trend is detected in TWSC of the Aral and Mekong River Basins. In recent decades, irrigation is highly developed in both of these basins [? ? ]. Analysis of the water flow development in Aral Basin indicated that the water diversion and irrigation schemes in this region have considerably increased evapo-transpiration and consequently, decreased net water flux from the atmosphere to the surface of the Aral Basin [? ]. These factors cause a reduction in the magnitude of seasonal amplitude of TWSCs within both of these basins. Sun et al. [? ] used GRACE data, flood potential index, as well as precipitation and other meteorological data to generate monthly terrestrial water storage anomalies and evaluated the flood potential in the Yangtze River Basin. Their results indicated that the amounts of water storage increased from 2003 to 2014, which confirms the increasing linear trend in the reconstructed TWSC of this study (numerical values of trends and accelerations are reported in Table ??). Figure 12. Basin average of the reconstructed TWSC fields within the seven major river basins of South Asia covering 2003-2018. The red lines show the long-term linear trend (L-T (mm/yr)) and acceleration (Acc (mm/yr 2 )) components fitted to the TWSC estimates.

Australian:
The two major basins of Murray and Lake Eyre are assessed here with no significant linear trend; however, the acceleration component is found to be negative for both basins. Lakshmi et al. [? ] analyzed the variability of individual water balance components in the major river basins of the world between 2002 and 2015. They showed that the maximum precipitation anomaly, as well as runoff and soil moisture anomaly within the Murray River Basin happened during 2011. Their result agree with the small positive trend and negative acceleration component of reconstructed TWSC within the Murray River Basin. Figure ?? shows that the the peak of TWS changes in Murray River Basin was happened in 2012. The associated time lag between maximum precipitation anomaly in 2011 (shows by [? ]), and the maximum TWS changes in 2012 could be influenced by different factors such as climatic conditions, vegetation density, runoff path-length, different recharge characteristics, permeability of sediments, and the abundance of surface water bodies (see e.g., [? ? ? ]). Forootan et al. [? ] argued that the increasing rainfall trends over the last 30 years, such as heavy rainfall in 2010 and 2011 over Murray and Lake Eyre River Basins, were mainly due to La Niña events. A decrease in 2012-2018 TWSC caused a negative acceleration in the reconstructed TWSC estimates of these two basins, which can be a response to the Millennium meteorological drought of 2000-2009 (see, e.g., [? ] and www.water.vic.gov.au/climate-change/millennium-drought-report).

Summary and Conclusions
In this study, two statistical reconstruction formulations (Approach 1 and Approach 2) are provided, which can be used to merge time-variable satellite gravity solutions with distinct accuracy and signal spectral contents. Here, these approaches are implemented to merge the total water storage change (TWSC) products of two mission GRACE (2003GRACE ( -2017 and the three-satellite constellation Swarm (2013-now). This integration bridges the one year gap between GRACE and GRACE Follow-On mission, providing a consistent 1 • × 1 • fields covering 2003-2018.
The formulation of Approach 1 follows a common data reconstruction procedure similar to [? ], and Approach 2 is an iterative reconstruction using the whole period of available GRACE and Swarm data. For both formulations, we applied independent component analysis (ICA, [? ? ]) to reduce the noise in original fields. We tested these techniques numerically in various ways, i.e., in terms of global empirical base functions (shown in Figure ??), in the global grid domain (Figures ?? and ??), and finally in terms of basin averages (Section ??).
All investigations show that both of the reconstruction techniques decrease the noise level in Swarm gravity fields and improve their spatial resolution. The results of Approach 2 provide more homogenized (long-term) TWSC fields, which are recommended for studies addressing global mass redistribution. Otherwise, the considerably high level of noise in Swarm data and their inconsistencies with GRACE(-FO) does not allows their applications to study long term hydrological changes in river basins (see e.g., Figure A1). Comparisons of the reconstruction results (Approach 2) with those of GRACE-FO in 2019 show that the magnitude of differences is smaller than 10 mm; therefore, these TWSC fields are used to monitor hydrological changes of the world's 33 largest river basins covering 2003-2018, where we find: • positive linear temporary trends in the St. Lawrence and Mississippi River Basins, caused by changes in surface water compartment and confirmed by altimetry measurements as an independent comparison; • negative linear temporary trends in TWSC of the Yukon, Nelson, Mackenzie, and Colorado River Basin, which can be associated to both irrigation and global warming; • a considerable decrease in the magnitude of the seasonal cycle within the four investigated basins in South America (Amazon, Orinoco, Tocantins, and Parana) that is attributed to the drought events of 2016-2018; • a mixed positive and negative temporary trends in the basins within Africa, where that of Nile was dominated by the surface water in lakes such as Victoria, but that of Niger was associated to groundwater; and those of Congo, Okavango, Orange, and Zambezi were influenced by an increase in the rainfall; • the impact of anthropogenic impacts on Euphrates, Ganges, Indus, and Brahmaputra, Murray (negative trend) as well as Amur (positive trend).
The values of these inter-annual climate episode driven trends in the TWSCs are reported in Table ??.
The proposed reconstruction methodology can potentially replace the contemporary techniques trying to optimally derive the relative data weights between two different data types. The reconstructed fields of Approach 2, which was used to estimate basin averages within the world's 33 largest river basin, will be provided upon request.

Appendix A. Impact of Iterations in Approach 2 on Independent Components
A comparison between TWSC estimations derived from the iteration 1 and 250 of Approach 2 was performed. The basin averaged results are shown in Figure ?? for three arbitrary selected river basins of Indus, Mackenzie, and Orinoco. It is clear from the graphs that the noise in Swarm data disturb estimation of trend, acceleration, and even interpretation of the seasonal scale changes in TWS; therefore, Swarm TWSC fields without implementing a reasonable post-processing cannot be used along-side GRACE data to assess hydrological changes within river basins. The noise has however almost vanished in the 250th iteration, where the red curves indicate smooth evolution of TWSC in the river basins as interpreted in Figures ?? to ??. Figure A1. A comparison between basin averaged TWSC estimations derived from the iteration 1 and 250 of the reconstruction-Approach 2. The results are shown for three arbitrary selected river basins.
In order to illustrate the impact of the iterative reconstruction on the TWSC estimations, here we select two arbitrary points in the USA (latitude: 40.5 and longitude: −85.5) and South Africa (latitude: −9.5 and longitude: 24.5). The time series of TWSC from the reconstruction were compared with those of GRACE, GRACE-FO, Swarm. The results in Figure ?? show that the reconstruction, not only fills the gap between GRACE and GRACE-FO, it also homogenizes the temporal evolution of TWSC. This is clear by comparing the original TWSC from Swarm (marked by stars) with those of reconstruction (marked by empty circles). Here, the gap periods are highlighted by colored boxes.

Appendix B. Noise and Leakage Estimation of the Reconstructed TWSC Fields
In this Appendix, we first assessed the level of noise in the original TWSC (i.e., smooth GRACE and Swarm TWSC) and reconstructed fields (compare Figure ?? top with bottom). For this estimation, the noise level is considered to be corresponded with 5% of the less dominant eigenvalues as discussed in Approach 1 and 2. The results indicate that the noise pattern of > 10 mm (which is a common noise level expectation from smoothed GRACE TWSC fields) in terms of TWSC vanishes in the reconstruction fields to < 3 mm. In Figure ??, an estimation of the spatial leakage error (computed following [? ]) is shown for three river basins of different size, i.e., Amazon as the largest river basin, and the area of Amur and Aral is relatively smaller. Numerical results indicate that the leakage error during the Swarm era (2013-2018) has been considerably decreased. The error in Amazon is decreased from an average of ∼ 3 mm to less than ∼ 1 mm, while in Amur this is found to be decreased from ∼ 9 mm on average to ∼ 4 mm, and finally in Aral, from ∼ 19 mm to ∼ 7 mm. An overall overview of the enhancement in the estimation of basin-averaged spatial leakage is shown in Figure ??, where the estimations correspond to the original GRACE (covering 2003-2014) and Swarm data (2013)(2014)(2015)(2016)(2017)(2018). We also show the estimated standard deviation of leakage errors for the same periods but from the reconstructed results (Approach 2). Our numerical investigations indicate that the standard deviation (StD) of leakage errors that correspond to the Swarm data (with an average of ∼ 15 mm) is ∼ three times bigger than that of GRACE (with an average of ∼ 5 mm covering 2003-2014). After the reconstruction, very comparable level of leakage errors is derived over the GRACE and Swarm era (i.e., on average ∼ 3.5 mm).

Acknowledgments:
The authors are grateful to the comments provided by the Associate Editor and four anonymous reviewers, which are used to improve the quality of this work. We acknowledge data product providers of GRACE, Swarm and altimetry products used in this study. E. Forootan thanks DAAD for the fund 57445178. N. Mehrnegar is grateful to the Cardiff University's Vice Chancellor fund, which supports her PhD research. C.K. Shum, Yu Zhang, and Chaoyang