remote sensing Determination of Weak Terrestrial Water Storage Changes from GRACE in the Interior of the Tibetan Plateau

: Time series of the Gravity Recovery and Climate Experiment (GRACE) satellite mission have been successfully used to reveal changes in terrestrial water storage (TWS) in many parts of the world. This has been hindered in the interior of the Tibetan Plateau since the derived TWS changes there are very sensitive to the selections of different available GRACE solutions, and ﬁlters to remove north-south-oriented (N-S) stripe features in the observations. This has resulted in controversial distributions of the TWS changes in previous studies. In this paper, we produce aggregated hydrology signals (AHS) of TWS changes from 2003 to 2009 in the Tibetan Plateau and test a large set of GRACE solution-ﬁlter combinations and mascon models to identify the best combination or mascon model whose ﬁltered results match our AHS. We ﬁnd that the application of a destriping ﬁlter is indispensable to remove correlated errors shown as N-S stripes. Three best-performing destriping ﬁlters are identiﬁed and, combined with two best-performing solutions, they represent the most reliable solution-ﬁlter combinations for determination of weak terrestrial water storage changes in the interior of the Tibetan Plateau from GRACE. In turn, more than 100 other tested solution-ﬁlter combinations and mascon solutions lead to very different distributions of the TWS changes inside and outside the plateau that partly disagree largely with the AHS. This is mainly attributed to less effective suppression of N-S stripe noises. Our results also show that the most effective destriping is performed within a maximum degree and order of 60 for GRACE spherical harmonic solutions. The results inside the plateau show one single anomaly in the TWS trend when additional smoothing with a 340-km-radius Gaussian ﬁlter is applied. We suggest using our identiﬁed best solution-ﬁlter combinations for the determination of TWS changes in the Tibetan Plateau and adjacent areas during the whole GRACE operation time span from 2002 to 2017 as well as the succeeding GRACE-FO mission. for supplying


Introduction
The Tibetan Plateau (TP) (25 •~4 0 • N and 75 •~1 05 • E) is both the highest and the largest plateau in the world. Since it is regarded as 'the water tower of Asia' and is very sensitive to global change, it is of major importance to understand the terrestrial water storage (TWS) changes in the plateau, not only for the utilization of water resources in the plateau and surrounding Asian countries but also for global change studies.
Terrestrial in situ hydrology observations to determine TWS changes have rarely been made due to the harsh weather conditions and the challenging geography in the  [19] CSR SH RL05 08/2002-12/2016 Cham P4M6 G250 3(a) 2, W/N Jing et al. [13] 3 SH models 01/2003-12/2016 -LOESS 5(a-c) 1, M/N Liu et al. [20] GRGS denote commonly used destriping methods by Duan et al. [22], Swenson and Wahr [3], and Chambers [23] respectively, PxMy takes polynomial degree to be x and start order for destriping to be y; G300 represents a Gaussian filter [4] with an averaging radius of 300 km, and LOESS stands for local regression smoothing [13]; numbers in the last column-the number of the anomalies of TWS changes inside the plateau, E/N-the single anomaly or the larger one of anomalies appears in the eastern Plateau by north, and W/N and M/N indicate that the single anomaly appears in the west and middle by north, respectively. Maximum d/o is 60 for the GRACE SH solutions but the Maximum d/o used by Liu et al. [20] is 80, while the size for GRACE mascon models CSR_M, JPL_M and GSFC_M are 1 • , 3 • and 1 • , respectively.
Erkan et al. [14], for example, used the CSR SH solutions [24] together with the destriping method by Duan et al. [22] to find two signal anomalies in the plateau, a smaller one in the west and a larger one in the east. Jiao et al. [16] instead used the level-3 data products derived from JPL SH solutions in combination with the destriping method by Swenson and Wahr [3] and found only one anomaly in the eastern plateau. If no destriping filters are applied, the number of the derived TWS anomalies usually increases [12,18,20]. We also noticed, solely by visual inspection, large discrepancies between the CSR_M, JPL_M and GSFC_M mascon model results in Jing et al. [13] and Loomis et al. [21].
We can summarize at least two reasons for the different previous results and the challenge to extract TWS changes in the TP: (1) Measurement noises, orbit errors and different strategies for producing GRACE solutions condition the fact that solutions from different agencies vary conspicuously; (2) Compared to other places, the hydrological conditions and their changes in the TP and adjacent areas are so special that the magnitudes of the TWS changes inside the plateau are small, and so can be easily falsified by the superimposed interference signals from nearby strong TWS signals-especially when using destriping and spatial smoothing filters for SH solutions. Hence, an adequate postprocessing of GRACE data must be tailored for clear determination of TWS changes in the TP.
The controversial TWS changes summarized in Table 1 may confuse readers and make them wonder which TWS change results are reliable in the TP. Based on currently existing GRACE solutions (including SH solutions and mascon models) and various destriping filters (e.g., destriping methods by Swenson and Wahr [3] and Duan et al. [22]), our aim is to identify the best GRACE solution-filter combinations and to obtain most reliable and reasonable TWS changes for the TP. This would close discussions on this controversial subject and provide guidelines for future studies on TWS changes in the TP and adjacent areas.
The criterion for choosing a certain solution-filter combination should be given by independent data serving as some kind of 'ground truth'. Considering the harsh weather conditions, challenging geography and thus rare in situ measurement of TWS changes in the TP, an aggregated TWS hydrology signal of the few available GRACE-independent observations of contributing signals is suggested herein. GRACE-derived TWS changes in the TP consist of changes in various water storage components, including soil moisture, glacier mass, snow water storage, lake water, permafrost and groundwater [2]. These contributions can be obtained from previous studies based on multiple altimetric missions, optical satellite imagery, hydrological models (e.g., GLDAS [25]) and changes in permafrost degradation calculated from an Active-Layer Depth (ALD) model [14,26]. The results from previous works are obtained through rigorous methods and techniques, and they usually have relatively small errors [27][28][29][30] or have been proven a good applicability in the TP [31,32]. Given the small errors, and that any additional water storage components in the TP have not been identified yet, the aggregated hydrology signals (AHS) of these different results should represent the currently best available 'ground truth' criterion for GRACE-derived TWS changes with relatively low spatial resolution (300-400 km). Moreover, many previous studies that also used a combination of GRACE and satellite altimetry and hydrological models for water storage change investigations got very good results [2,12,16,[33][34][35].
We test about 130 combinations of GRACE solutions and destriping filters and some mascon models until the derived TWS changes agree with the AHS of the same period. We note that all these solutions, destriping methods and models are known and well established, thus we do not aim to introduce a new filtering approach or model.
In the next section, the used GRACE solutions and destriping filters are briefly introduced. In Section 3, the calculation of the GRACE-derived TWS changes and the AHS are shown, followed by the presentation and discussion of the results in Section 4. Finally, we summarize our findings in Section 5.

GRACE Solutions
We apply fifteen different newly released monthly GRACE solutions, i.e., twelve SH solutions and three mascon models. The time span covers the period from April 2002 to June 2017. Twelve newly updated SH solutions are ITSG [36], CSR [24], COST-G [37], JPL [38], GFZ [39], Tongji [40], WHU [41], HUST [42], GRGS [43], IGG [44], AIUB [45] and XISM [46]. The different SH solutions are provided with Stokes coefficients of different maximum degrees such as 60, 90 and 96. We made tests to ensure the proper maximum degree since the Stokes coefficients become noisier with higher degrees, which needs an effective destriping filter to depress the N-S stripe noises. The glacial isostatic adjustment (GIA) effects are removed from the monthly Stokes coefficients according to the GIA model ICE-6G_C(VM5a) [47].
The three global equal area mascon models are CSR_M [6], JPL_M [48] and GSFC_M [49]. CSR_M is derived from GRACE information only [6], while JPL_M contains next to GRACE information and a-priori information of mass variations as constraints, which are from land, ocean, ice, inland sea, large inland lake, earthquake and glacial isostatic adjustment [48]. GSFC_M uses the signal auto-covariance matrix as a regularization constraint, and is constructed by spatio-temporal constraint equations [49,50]. To be applicable for comparison in this study, we decompose the mascon grid EWT data into SH Stokes coefficients from degree 1 to 90.
Considering the absence of degree-1 coefficients and inaccuracy of C20 and C30 coefficients for SH solutions, we complement GRACE SH solutions with newly released C10, C11 and S11 coefficients [51], and replace the native GRACE C20 and C30 coefficients with the latest SLR-derived version [52]. We note that the C30 coefficients were replaced when the time series passed March 2012, where the SLR-derived C30 starts in the data set according to Loomis et al. [52].

Destriping Filters
The monthly Stokes coefficients of each GRACE SH solution are destriped with three frequently used destriping methods proposed by Duan et al. [22], Swenson and Wahr [3] and Chambers [23], respectively. For each of the three methods, a polynomial of selected degree x is used in a window of width w, fitted to all the coefficients of the same parity in degrees ≥y for a certain order (≥y), and the coefficients used in fitting are subtracted by the fitting values. In this way, the correlated errors which manifest as N-S stripes can be reduced in the coefficients within the window. Thus, a series of the destriping filters named 'Duan PxMy', 'S&W PxMy' and 'Cham PxMy' in this study can be devised and used for x = [ 2,3,4] and y = 8. These values were found to perform well in earlier studies [53]. The filters 'Duan P4M10', 'Duan P4M15', 'S&W P3M10' and 'S&W P3M15' are additionally considered in order to observe the effects of different y values.
Destriping filters 'Duan PxMy' and 'S&W PxMy' use moving filtering windows. According to Duan et al. [22], the width w for 'Duan PxMy' filtering is calculated using Equation (1), dependent on the degree l and order m of the coefficients to be filtered and the parameters (γ, p, A, K) taking empirical values (0.1, 3, 30, 15), However, for 'S&W PxMy' filtering, w is calculated using Equation (2), dependent on the order m of the coefficients, and the parameters (A, K) taking the values (30,10), For 'Cham PxMy' filtering, the windows include all the coefficients of the same parity in degrees from y to 90 (or 96) for a certain order (also ≥y).
Piretzidis et al. [54] developed a method of identifying and selectively destriping SH coefficients with correlated errors by utilizing the capabilities of machine learning algorithms (MLAs), which can be used to calculate monthly changes in TWS [54], and this destriping method is marked by 'MLAs' here.
The monthly Stokes coefficients can also be destriped and smoothed simultaneously for each GRACE gravity SH solution with the non-isotropic filtering approach proposed by Kusche et al. [5]. We choose the filter DDK8 for four gravity SH solutions (ITSG/CSR/JPL/GFZ) and the much stronger filters DDK2 for the preferred ITSG model to show their performances, respectively. According to Kusche et al. [5], the filtering characteristics of the DDK2 filter are approximately equivalent to a Gaussian filter with an averaging radius of 340 km.

GRACE-Derived TWS Change
According to Wahr et al. [55], the monthly TWS changes (in EWT) can be derived at an any grid site with co-latitude θ and longitude φ on 1.0 • × 1.0 • grids covering our study area with 2l + 1 1 + k l × (∆c lm cos mφ + ∆s lm sin mφ) P lm (cos θ), (3) in which ∆c lm and ∆s lm are the changes in Stokes coefficients from GRACE gravity SH solutions or decomposed from three mascon models, P lm (cos θ) is the lth degree and mth order normalized associated Legendre polynomial, k l is the degree l elastic load Love number of the Preliminary Reference Earth Model for potential perturbation [56], ρ ave is the average density of the earth, ρ w is the water density and a is the radius of the Earth. The derived TWS changes ∆T TWS (θ, φ) are expressed as EWT converted to the unit of millimeter. In order to analyze the features of the monthly TWS changes in the TP and adjacent areas, the trend rate, annual change at each grid point, together with a semi-annual and a 161 days aliasing variation (from the S2 semidiurnal solar tide), are simultaneously estimated using a least-square linear regression.

Aggregated Hydrology Signal
The AHS (denoted by trend rates of EWT changes) from 2003 to 2009 (and 2016) represent the total TWS changes from soil moisture (SM), glacier ice (ICE), lake water (LAKE) and permafrost (PM) in the TP and adjacent areas. Groundwater storage (GWS) contributions have not yet been assessed independent of GRACE data in this region and thus are not included in the AHS. Signal changes caused by snow and ICE on glaciers were coupled following Gardner et al. [28], which we marked as ICE, while internal testing yielded that snow out of glaciers has a negligible impact. Hence, snow is not added in the AHS. Further, we assume that the AHS represent the basic features of realistic TWS changes inside the plateau except for the eastern plateau where the increase in GWS signals is too strong to be neglected as a component of TWS changes [2]. Hence, the eastern plateau is omitted when determining the best GRACE solution-filter combination.
For comparisons of the GRACE-derived TWS changes with the AHS, the contribution of each hydrology component is also derived with Equation (3), summed till degree 60 (and 90, for evaluating the reliability of the simulated TWS changes derived from different maximum d/o). However, when deriving the contributions, mass change for the hydrology component K (= SM, ICE, LAKE or PM) can be represented as ∆W K , and thus, the changes in Stokes coefficients ∆c K lm , ∆s K lm are given by decomposition of corresponding grid data ∆W K designed referring to xiang et al. [2] in our surveyed area (Ω) using with dΩ = sin θ dθ dφ . The integrations of Equation (4) can be easily solved by summing the contributions of all the grids for SM, ICE, LAKE and PM, respectively. The hydrology components can be given by trend rates of EWT changes from 2003 to 2009. The grid SM data are from three different land surface models: Noah, Variable Infiltration Capacity (VIC) and Catchment Land Surface Model (CLSM) of GLDAS-2.1 which is forced with a combination of model and observation data from 2000 to present [25]. The LAKE data for 219 lakes are combined from Zhang et al. [27] and Farinotti et al. [29], and the ICE data are taken from Gardner et al. [28]-all are based on the analysis of ICESat-1 data. Since grid PM data are not available during our surveyed period, we follow Xiang et al. [2] and use trends of active-layer depth in the plateau from 1980 to 2001 by Oelke and Zhang [26] and Erkan et al. [14]. For comparison purposes over a longer period, the hydrology components are also given by trend rates of EWT changes from 2002 to 2016, using the grid SM data from Noah, VIC and CLSM of GLDAS, PM data based on Xiang et al. [2], LAKE data for 61 lakes derived from Li et al. [57] and LEGOS HYDROWEB [30] based on multi-source remote sensing, and the ICE data from Brun et al. [58] using time series of digital elevation models derived from ASTER optical satellite stereo-imagery.  Tables 2 and A1), which are intended to show their pattern and magnitude agreements of observed and aggregated models. PCC-values, RMSE-values and IOA-values are calculated using Equations (A1)-(A3) respectively. Table 2. PCC-, RMSE-and IOA-values between the GRACE-derived TWS trends of selected GRACE solution centers (as in Figures 1a-l and 2a-h) and the aggregated hydrology signals in our chosen comparison area inside the plateau (Figure 3f). Additional results for 7 more solution centers can be found in Table A1. The best-performing filters are separated by two dashed lines. In order to verify the effectiveness of the selected best solution-filter combination(s) over a longer period, we also obtained PCC-, RMSE-and NSE-values between GRACE-  In an area of the plateau interior delimited by red lines in Figure 3f, the AHS can be used as a ground truth for selection of the favorable GRACE solution(s) and destriping filter(s) for the TP interior. This marked target area is distant from the surrounding (glaciated) mountains where strong hydrology signals may impact the GRACE-derived TWS changes inside the plateau through signal leakage. It also excludes the eastern part of the plateau where possible obvious GWS changes [2] could not be modeled. In Tables 2 and B1, [25], (b) ICESat-1 based glacier measurements [28], (c) ICESat-1 lake measurements and LEGOS Hydroweb (e.g., [27,29,30]), (d) one ALD model [14,26], (e) ICE-6G_C(VM5a) model [47] and (f) aggregated hydrology signals representing the sum of four hydrological components (a-d).

Filter Dependence of the GRACE-Derived TWS Trends
In Figure 1a-k, the results of the GRACE-derived TWS trends from the best performing ITSG solution show obvious N-S variability depending on the selected destriping method, and E-W variability depending on the selected filter parameter (polynomial degree x = {2,3,4}). The results can be compared with those with original N-S stripes for the not-destriped case (Figure 1l). Particularly, they are compared with the AHS in the selected region inside the plateau (Figure 3f) to check their pattern and magnitude agreement. No obvious N-S stripes and better pattern and magnitude agreement (e.g., PCCvalues ≥0.66, RMSE-values≤5.95 and IOA-values ≥0.76) would imply a reasonable selection of the destriping filter(s).
In an area of the plateau interior delimited by red lines in Figure 3f, the AHS can be used as a ground truth for selection of the favorable GRACE solution(s) and destriping filter(s) for the TP interior. This marked target area is distant from the surrounding (glaciated) mountains where strong hydrology signals may impact the GRACE-derived TWS changes inside the plateau through signal leakage. It also excludes the eastern part of the plateau where possible obvious GWS changes [2] could not be modeled.
In Tables 2 and A1,   We compare the destriping method of Chambers [23] with a non-moving window but with variable widths depending on the start order y (Figure 1g-i) with the non-destriped case (Figure 1l). For filters with x = 2, 3 and 4, visible N-S stripes remain inside and outside the plateau. As we also find much smaller PCC-and IOA-values, and much larger RMSEvalues than those mentioned above, we think that these destriping filters are less effective. In Figure 1k, the results of the DDK8 filter also show some N-S stripes especially outside the plateau, and have a smaller PCC-value of 0.34, a smaller IOA-value of 0.53, and a larger RMSE-value of 9.03. This weaker filtering was already mentioned by Kusche et al. [5]. In Figure 1j, although there are no obvious N-S stripes inside the plateau, PCC-and IOA-value are small, the RMSE-value is so large that they also imply a less effective destriping filter.
We have further tested four destriping filters 'Duan P4M10', 'Duan P4M15', 'S&W P3M10' and 'S&W P3M15' having higher start orders than our previous tests. These results ( Figure A1) are very similar in the pattern and the amplitude to those for the filters 'Duan P4M8' and 'S&W P3M8'. However, the PCC-and IOA-values are smaller and RMSEvalues are larger than those for 'Duan P4M8' and 'S&W P3M8', respectively, and become rapidly smaller (for PCC-and IOA-values) and larger (for RMSE-values) as the start orders become larger from 8 to 15, thus we suggest that the start order y = 8 is adequate for the tested destriping filters, which was also suggested by Swenson and Wahr [3] when they introduced the filter. Figure 2 shows TWS trends from the five GRACE SH solutions ITSG, CSR, COST-G, JPL and GFZ filtered with one of the best performing filters, 'S&W P3M8', and from the three mascon models CSR_M, JPL_M and GSFC_M. The results can be compared with the AHS (Figure 3f) in the selected region inside the plateau.

Solution Dependence of the GRACE-Derived TWS Trends
In Figure 2a Table 2). The three mascon model results (Figure 2f-h) show much fewer similarities and agreement to the AHS than those derived from the SH solutions. This is already indicated by their small PCC-values between −0.04 and 0.13, small IOA-values between 0.37 and 0.42, as well as by their large RMSE-values between 8.75 and 11.48. The single anomaly found in the middle of the plateau obviously disagrees with the anomalies of the AHS. This could be due to less effective suppression of the N-S stripes in the processing of the mascon solutions, which will be discussed in Sections 4.3 and 4.4.
Comparing with other SH solutions listed in Table A1 and combining them with filter 'S&W P3M8', ITSG and CSR solutions also show the highest degree of correlation and agreement but smallest bias with the AHS. When considering the RMSE-values only, Tongji and WHU solutions are also good, but values of the other two indicators are slightly worse than those of SH solutions in Table 1. Other SH solutions in Table A1 perform much worse than the ones discussed here. When using the other two best performing destriping filters (i.e., 'Duan P4M8' and 'S&W P2M8'), ITSG and CSR solutions still perform best.

Determination of GRACE-Derived TWS Changes inside the Plateau
Here, the GRACE-derived TWS changes inside the plateau, before and after smoothing with a 340-km-radius Gaussian filter, are determined. Furthermore, we identify negative signals (~25.2 mm/yr) in the Tien Shan due to increased ice melting, decreasing soil moisture, and negative signals from Northwest India (~55.0 mm/yr) to Bengal Basin (~34.0 mm/yr) due to excessive use of groundwater [2], superimposing the effect of the increasing ICE melting along the Himalayas as well as soil moisture reduction there.
The anomalies inside the plateau of the TWS trends from the three mascon models (Figure 2f-h) obviously disagree with the AHS (Figure 3f). Moreover, the anomaly patterns of the three mascon models and their magnitudes disagree with each other inside and outside the plateau. For example, the negative anomaly due to the over-exploitation of groundwater in Northwest India has a magnitude of 82 mm/yr for JPL_M, which is more than 80 percent larger than the 45 mm/yr determined for GSFC_M. The one for CSR_M is 66 mm/yr, also much smaller than that for JPL_M.

Annual Changes
The annual amplitudes ( Figure A2a-d) and annual phases ( Figure A3a-d) from our best performing GRACE SH models ITSG and CSR together with the destriping filters 'Duan P4M8', 'S&W P2M8' and 'S&W P3M8' (except for CSR together with S&W P2M8) yield very consistent results. The annual amplitude and annual phase for CSR together with 'Duan P4M8' are not shown here as they are nearly the same as those for CSR together with 'S&W P3M8'. Amplitudes of~130 mm,~60 mm and~100 mm appear in the west, middle, southeast of the plateau, respectively, with almost the same values for the best solution-filter combinations. Positive phases in the middle of the plateau accompanied with negative phases also coincide very well.
We also show annual changes from other destriping filters ('MLAs', 'DDK8' and 'Cham P4M8' with ITSG SH solution) and three mascon models in Figure A2e-j and in Figure A3e-j, which vary extraordinarily. The results are very different from those derived from our best performing solution-filter combinations, suggesting that they are less applicable for TWS change investigations in the TP.

TWS Changes with Smoothing
We now apply additional Gaussian filter smoothing to our best performing GRACE solution-filter combinations. Very similar results are obtained for trend rates (Figure 4a-d), annual amplitudes ( Figure A4a-d), and annual phases ( Figure A5a-d), respectively. A single positive anomaly in the trend rates is found inside the plateau with the center in the east around the Three-River-Source Region and Qaidam Basin, dragging a small tail showing slight increase in the West Kunlun region, which agrees with previous works (e.g., [2,16,53]). The pattern of the positive anomaly is also, to a large degree, supported by the Gaussian filter smoothed AHS (Figure 4l Figure 4g where two anomalies appear with one in the east and the other in the west. Obviously, those are different from one single positive anomaly with the center in the east as for the determined smoothed trends from the best combinations (Figure 4a-d). Such single anomalies in the center of the plateau, as found with less effective destriping filters and the three mascon models, present signal patterns with shapes and amplitudes similar to Figure 4k where destriping was omitted and the result is only smoothed by a 340-km-radius Gaussian filter. Hence, these results may result from not effectively suppressed stripe noise. These results as well as the occurrence of two anomalies (Figure 4g) are also not supported by the Gaussian filter smoothed AHS (Figure 4l). the east around the Three-River-Source Region and Qaidam Basin, dragging a small tail showing slight increase in the West Kunlun region, which agrees with previous works (e.g., [2,16,53]). The pattern of the positive anomaly is also, to a large degree, supported by the Gaussian filter smoothed AHS (Figure 4l), although the positions and amplitudes of the trend anomalies (Figure 4a-d) are not completely consistent with those in Figure 4l in some special regions, such as in the Three-River-Source Region and Qaidam Basin and Northwest India, due to absence of groundwater information in these regions.  Figure 1 but for some selected GRACE solutions and destriping filters smoothed by a 340-km-radius Gaussian filter (a-j), compared with the signals from the ITSG solution without destriping (k) and aggregated hydrology signals (l), but smoothed by the same Gaussian filter. The Gaussian filter is not used for ITSG+DDK2 (f), because the DDK2 filter has the smoothing characteristics of a 340-km-radius Gaussian filter according to Kusche et al. [5]. The purple lines in (a-d)  Figure 1 but for some selected GRACE solutions and destriping filters smoothed by a 340-km-radius Gaussian filter (a-j), compared with the signals from the ITSG solution without destriping (k) and aggregated hydrology signals (l), but smoothed by the same Gaussian filter. The Gaussian filter is not used for ITSG+DDK2 (f), because the DDK2 filter has the smoothing characteristics of a 340-km-radius Gaussian filter according to Kusche et al. [5]. The purple lines in (a-d) delimit the range of the the Three-River-Source Region (below) and Qaidam Basin (above), and the dark brown rectangle represents the West Kunlun region.

Assessment of the Effectiveness of Destriping Filters and GRACE Solutions
We investigate the effectiveness of destriping for the tested destriping filters and three mascon models by a combined analysis of the filtered and unfiltered (Figure 1l) TWS trends.
First, we focus on the unfiltered results (Figure 1l), in which the anomalies inside and outside the plateau are elongated in the north-south direction, and thus are extremely impacted by the N-S striping noise. According to Wahr et al. [55], a Gaussian filter deals with the stochastic noises in the higher-degree Stokes coefficients, making the filtering equivalent to a spatial smoothing of the TWS trends. The N-S stripes should be one kind of system-noise errors, which should be suppressed with a dedicated destriping filter such as 'Duan P4M8', 'S&W P2M8' and 'S&W P3M8'. Therefore, although the N-S stripes are not apparently visible in the smoothed result (Figure 4k), we cannot expect that the effects of N-S stripes are fully suppressed when only a Gaussian filter is used. We thus consider the TWS trends shown in Figure 4k, smoothed but not destriped, as a benchmark for indicating the large effects of N-S stripes.
Since N-S stripes in the TWS trends from the five best solution-filter combinations (i.e., 'ITSG+Duan P4M8', 'ITSG+S&W P2M8', 'ITSG+S&W P3M8', 'CSR+Duan P4M8' and 'CSR+S&W P3M8') before smoothing (Figures 1c-e and 2a,b) are already almost imperceptible, the corresponding Gaussian smoothed TWS trends (Figure 4a-d) certainly show different patterns from our chosen benchmark result (Figure 4k). Note that the TWS trends from best solution-filter combinations before smoothing also agree well with the AHS (Figure 3f), showing the larger PCC-and IOA-values, and smaller RMSE-values. Figures 1c-e, 2a,b and 4a-d thus present reasonable GRACE-derived TWS trends before and after smoothing, respectively, in which the N-S stripes are effectively suppressed. Due to extreme similarity with those from 'CSR+S&W P3M8', TWS trends from 'CSR+Duan P4M8' are thus not shown in these Figures.
For the combination series 'ITSG+Cham PxM8', obvious stripes in the TWS trends before smoothing are seen in Figure 1g-i with very small PCC-and IOA-values, and very large RMSE-values. The pattern of the TWS trends after smoothing (Figure 4g) is not supported by the smoothed AHS (Figure 4l).
From the DDK series of filtered data products, we select two representative filters, DDK8 and DDK2 with weak and strong filtering, respectively. As seen in Figure 1k, the TWS trends of DDK8 are less effectively destriped, with anomalies stretched in the N-S direction visually and small PCC-and IOA-values as well as large RMSE-values. For the DDK2 filter, equivalent to a 340-km-radius Gaussian filter [5], the similarity of the derived TWS trends (Figure 4f) to the benchmark results (Figure 4k) implies that the original N-S stripes seem to be smoothed just by a Gaussian filter, but not effectively suppressed. As for the NLAs filter, although there are no obvious N-S stripes in the TWS trends before smoothing (Figure 1j Figure A6 are smaller than those for the time span 2003 to 2009. This is due to limited availability of monitoring data of glacier melting and lake water changes as well as the lack of groundwater storage changes data for the longer period until 2016. For the three mascon models, the anomalies inside the plateau of the TWS trends before smoothing (Figure 2f-h) cannot match the three positive anomalies and one negative anomaly (Figure 3f), which shows bad pattern agreement of the results with the AHS. In addition, as mentioned above, the anomaly patterns of the three mascon models and their magnitudes do not match with each other inside and outside the plateau. Furthermore, the smoothed TWS trends (Figure 4h-j) show very similar patterns with those of the benchmark results (Figure 4k), which strongly implies that large effects due to N-S stripes are included. This, in turn, means that, although not visible, N-S stripes have already not been effectively suppressed in the TWS trends before smoothing (Figure 2f-h). This is probably due to some smoothing functions in the inversion processing of mascons. We thus argue that the tested mascon models are less reliable in the TP area.

Additional Tests and Discussions
Although the contributions of each hydrology component in the AHS have been well determined, as a criterion for selecting, it is necessary to investigate the influence of its uncertainties on identifying the best solution-filter combination or mascon model. Therefore, we estimated the error sigmas for each contribution (Figure A7a-d) referring to Xiang et al. [2] and calculated the error sigmas (σ) for the AHS ( Figure A7e) using Equation (A4) and a square root function. Further, four extreme cases of AHS ( Figure A7f-i) considering ±1σ and ±2σ as their uncertainties were simulated. Based on those, we conducted the same experiments as outlined in Section 4.1 but using the four extreme cases of AHS.
We can find in Figure A7e that the σ for AHS in the plateau interior (delimited by red lines, also in Figure 3f) are relatively small, avoiding large signal leakage from strong hydrology signals around. Therefore, when considering the four extreme cases 'AHS ± 1σ' and 'AHS ± 2σ', their signal patterns in Figure A7f-i are very similar to those of the AHS in Figure 3f.
The PCC-, RMSE-and IOA-values between the GRACE-derived TWS trends and the four extreme cases of AHS are listed in Tables A2-A5. In Tables A2 and A3, the PCC-, RMSE-and IOA-values are for the cases of 'AHS+1σ' and 'AHS-1σ', respectively. Similar to Table 2 (Tables A2 and A3). The PCC-value describes consistent proportional increases or decreases around the respective means of two models but only few differences between their magnitudes [59], thus reflecting mainly the similarity of signal distribution. Changes of ±1σ in AHS do not change the signal pattern much ( Figure A7f,g), thus the PCC-value is not altered much, but they change the differences between the observed signals and the AHS. These differences affect RMSE-and IOA-values, however, and thus they change somewhat significantly. The 'AHS+1σ' case (Table A2) yields smaller RMSE-values and better IOA agreement than the AHS, while the 'AHS-1σ' case (Table A3) results in larger RMSE-values and worse IOA agreement than the AHS. That means the 'AHS+1σ' would be closer to the TWS changes actually happening in the interior of the TP than the AHS, and much closer than the 'AHS-1σ' case. When adding or subtracting 2σ from AHS, the 'AHS-2σ' case (Table A5) performs even worse than the 'AHS-1σ' case, while in the 'AHS+2σ' case (Table A4), the indicator values provided both better and worse agreement than the AHS for different solutions and indicators used. These differing results of the cases may be due to the absence of groundwater increase in the TP [2] contributing to the AHS. This hypothesis should be addressed in a future study.
Concerning the validity of the SH solution and mascon models in the TP, Figure 2 shows striking disagreements between SH solutions and mascon models. TWS changes from the three mascon models do not reflect changes of soil moisture, glacier, lake water and permafrost in the interior of the TP. Mascon models appear to be applicable only to regions with simple or easily-distinguishable and large TWS changes, but our investigations show that they are not applicable to the TP where complex components contributions (e.g., soil moisture, glaciers, lakes, permafrost and groundwater) distributed in various places constitute the TWS changes, especially when parts are difficult to observe independent of GRACE, such as groundwater and permafrost [2]. The JPL_M mascon model takes advantage of near-global geophysical models to prevent longitudinal 'stripe patterns' in the solutions which are caused by GRACE's poor observability of east-west gradients [48]. While the GLDAS NOAH model is the only hydrological constraint used in the TP [48], GLDAS NOAH neither defines the ice-covered region nor depicts the sites and distribution of the lakes and other components contributions [25], which introduces regional errors. JPL_M uses 3 • equal-area spherical cap mascons, each usually covering an area with scattered small lakes, glaciers, permafrost, un-known groundwater and/or other components, which yield the large total water changes. As a mascon consists of a flat disk mascon unit mass, it is, in essence, a form of isotropic smoothing filter which removes longitudinal stripe patterns quite effectively visually, but in reality, it smoothes both stripe errors and the actual signal simultaneously. This explains the similarity of signal patterns between smoothed TWS trends (Figure 4i) and those of the benchmark results (Figure 4k).
CSR_M is developed with 1 • resolution hexagonal tiles mascons, based on regularized spherical harmonic GRACE solutions derived from GRACE information only [6]. The developers emphasize that the regularized spherical harmonic solutions benefit from the L-curve method with Tikhonov regularization, significantly reducing striping errors [6,60]. GSFC_M provides global mascon solution with 1 • resolution, resulting from an approach that combines an iterative solution strategy with specially designed regularization matrices for the mascon inversion, which is also reported to be highly effective by the developers [49]. Although CSR_M and GSFC_M are successfully used in many regions [6], the destriping effectiveness of their inversion processing is not supported in the TP given our results above.
As we know, some GRACE SH solutions provide Stokes coefficients d/o up to 96 (e.g., CSR) or even 120 (i.e., ITSG). In order to investigate the effects of higher d/o Stokes coefficients (e.g., from 61 to 90) on the GRACE-derived TWS changes, we conducted the same experiments as mentioned above in Sections 3 and 4 but using Stokes coefficients maximum d/o up to 90, and show corresponding results in Figures A8-A10.
It is found that nearly all results show obvious N-S stripes in the TWS trends from GRACE SH solutions regardless of different destriping filters applied ( Figure A8, with exception of Figure A8c Although there are no stripes found in the TWS trends from the three mascon solutions with higher d/o ( Figure A9f-h), all PCC-, RMSE-and IOA-values between the GRACE TWS signals and the AHS ( Figure A10f) suggest a strong influence of N-S stripes remaining.
Although this paper is devoted to the determination of the GRACE-derived TWS changes inside the plateau, Figures 1, 2 and 4 show that the TWS changes outside the plateau also strongly vary depending on the selected destriping filters as well as GRACE solutions. It appears that our best GRACE solution-filter combinations for the area inside the plateau could still be applicable for the areas outside the plateau. Hence, determined TWS changes there, as shown in Figures 1c-e, 2a,b and 4a-d, would reasonably reflect the glacier melting in the surrounding mountains of the Tien Shan, Himalayas and southeast of the plateau, as well as the over-exploitation of groundwater from Northwest India to the Bengal Basin. Note that the determined TWS changes outside the plateau also have similar characteristics of patterns and magnitudes in all best GRACE solution-filter combinations.
In view of the fact that GRACE Follow-On (GRACE-FO) data are of the same type and format as GRACE data, and are affected by the same sources of error, our findings above are considered to be applicable to GRACE-FO missions too. Several studies show, meanwhile, the excellent continuation of the GRACE data with GRACE-FO data [61,62]. GRACE-FO data are not included in this study due to limited availability of AHS components which hinders the extension of our investigation time span to the whole period of GRACE and GRACE-FO together.
Our results of the GRACE-derived TWS changes could be influenced by signal leakages due to the truncation of higher harmonic degree terms and necessary filtering (e.g., using a Gaussian filter in this study). However, it is beyond the scope of this study to remove the leakage effects to restore the realistic signals of TWS changes. Concerning the solid earth dynamics in the TP and adjacent areas, we corrected the effect of GIA on the GRACE-derived TWS changes with the GIA model ICE-6G_C(VM5a) carefully. The broad-scale tectonic uplift in the TP would be isostatically compensated by an increasing mass deficiency at depth and has little net effects on our results [63], while the erosion rates are too small and conservative [2]. Both effects are thus not considered in this study.

Conclusions
We identified the best GRACE solution-filter combination for identification of terrestrial water storage changes in the Tibetan Plateau. 'Duan P4M8', 'S&W P2M8' and 'S&W P3M8' are found to be best destriping filters as they remove the spatially correlated errors for all twelve GRACE SH solutions tested in here. These three destriping filters, together with the two GRACE SH solutions of ITSG and CSR, can be considered the best solution-filter combinations, since the inversed TWS changes from 2003 to 2009 inside the Tibetan Plateau agree well with the aggregated hydrology signals, with an exception of the CSR solution together with the 'S&W P2M8' filter. This is also confirmed when extending the end of the time span from 2009 to 2016.
The TWS changes, including trends and annual signals with and without smoothing by a 340-km-radius Gaussian filter, are well determined with our best GRACE solution-filter combinations, which yield very similar features in the signal pattern and very similar magnitudes inside and outside the plateau. We find three positive anomalies and two negative TWS trend anomalies inside the plateau which can also be identified in an aggregated hydrology signal. Inspection of smoothed results shows one single anomaly inside the plateau, elongated from east to west with its center in the east of the plateau and a small tail showing slight increase in the West Kunlun region, as found by Jiao et al. [16] and Xiang et al. [2,53]. We suggest that the best solution-filter combinations can also be applied to the areas outside the plateau and other observation time spans, and thus can help to reasonably determine TWS changes inside and outside the plateau during the GRACE lifetime from 2002 to 2017, as well as its Follow-On (GRACE-FO) mission.
Very different TWS changes inside and outside the plateau from those found with our best solution-filter combinations are derived with the Chambers, DDK8/DDK2 and MLAs filters, as well as the three mascon models. We think this is very likely due to the fact that the N-S stripes are not effectively suppressed, and thus these filters and products are recommended to be used only very carefully when investigating TWS changes in the Tibetan Plateau and adjacent areas.
Finally, noting that current destriping filters cannot depress the N-S stripes effectively for Stokes coefficients of d/o higher than 60, Stokes coefficients of only d/o 1-60 of GRACE SH solutions are suggested to be used for investigation of TWS changes.  The mascon model can be obtained from http: //www2.csr.utexas.edu/grace/RL06_mascons.html, the mascon model JPL_M from https://grace. jpl.nasa.gov/data/get-data/jpl_global_mascons/, and the mascon model GSFC_M from https: //earth.gsfc.nasa.gov/geo/data/grace-mascons. Degree-1 and C20 data are available via GRACE TN-13 (GRACE Technical Note 13, https://podaac-tools.jpl.nasa.gov/drive/files/allData/grace/docs/ TN-13_GEOC_CSR_RL06.txt) and GRACE TN-14 (NASA GSFC SLR C20 and C30 solutions, https: //podaac-tools.jpl.nasa.gov/drive/files/allData/grace/docs/TN-14_C30_C20_GSFC_SLR.txt), respectively. The ICE-6G_C model is available at http://www.atmosp.physics.utoronto.ca/~peltier/ data.php. and where O i and A i are the trends of TWS and AHS on grids, n is the number of grids in target area, and O is the mean of all O i . The PCC-value always has values between −1 and 1. A value of ±1 indicates a perfect degree of association between the two variables, a + (plus) sign indicates a positive relationship and a − (minus) sign indicates a negative relationship. When the correlation coefficient value approaches 0, the relationship between the two variables will be weaker. The RMSE-value summarizes the mean difference in the units of the observed and aggregated signals, which is a measure of how much the differences between the models are Remote Sens. 2022, 14, 544 20 of 33 spread out. The IOA-value ranges from 0 to 1. The closer the IOA-value is to 1, the better the model-data agreement is.

Appendix A.2. Error Estimation of the Aggregated Hydrology Signal
The AHS in this paper are derived from soil moisture (SM), glacier ice (ICE), lake water (LAKE) and permafrost (PM). Therefore, the error variances (σ 2 AHS ) for the AHS can be calculated by Referring to Xiang et al. [2], error sigmas at grids for SM (σ SM ) are given by the standard error of the mean (SEM) of Noah, VIC and CLSM models, error sigmas for ICE (σ ICE ) are calculated with uncertainties of glacier regions provided by Gardner et al. [28], while those for LAKE (σ LAKE ) and PM (σ PM ) are given by absolute values of 10 and 100 percent of the signals, respectively.     Figure A2 but additionally smoothed by a 340-km-radius Gaussian filter. The Gaussian filter is not used for ITSG+DDK2 (f). Figure A4. Same as Figure A2 but additionally smoothed by a 340-km-radius Gaussian filter. The Gaussian filter is not used for ITSG+DDK2 (f).