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

Abstract

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 filters 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-filter combinations and mascon models to identify the best combination or mascon model whose filtered results match our AHS. We find that the application of a destriping filter is indispensable to remove correlated errors shown as N-S stripes. Three best-performing destriping filters are identified and, combined with two best-performing solutions, they represent the most reliable solution-filter 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-filter 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 filter is applied. We suggest using our identified best solution-filter 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.

1. Introduction

The Tibetan Plateau (TP) (25°~40°N and 75°~105°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 plateau. Hence, time-variable gravity measurements of the Gravity Recovery and Climate Experiment (GRACE) twin-satellite mission since 2002 were warmly welcomed as they were seen as bright prospect in monitoring the regional TWS changes there (e.g., [1,2]).

Two types of time series of GRACE data products are provided which can be used for the estimation of TWS changes, (1) spherical harmonic (SH) gravity solutions consisting of monthly Stokes coefficients provided by different institutions, and (2) mass concentration (mascon) models, also provided by different institutions, which directly give the TWS changes in equivalent water thickness (EWT) on mascon grids.

When the TWS changes are derived based on a GRACE SH gravity solution, the precise determination of TWS changes is hampered on one hand by measurement noises, which increase with the harmonic degree, and on the other hand due to the prominent north-south-oriented (N-S) stripe effects from the orbital configuration of the twin satellites. It is thus indispensable to apply adequate post-processing with a destriping filter (e.g., the filter proposed by Swenson and Wahr [3]) and a smoothing filter (e.g., Gaussian filter proposed by Jekeli [4]), successively. Kusche et al. [5] proposed a non-isotropic filtering approach which can conduct destriping and smooth filtering simultaneously. In a mascon model, the stripe errors and measurement errors are effectively suppressed according to the developers and thus the users are advised to avoid any additional empirical destriping or smoothing (e.g., [6,7]).

While GRACE data products have been successfully used to derive TWS changes in many parts of the world (e.g., [8,9,10,11]), clear identification of TWS changes in the TP from GRACE is challenging. The magnitudes of potential TWS anomalies are small and in the order of about 15 mm/yr [2]. These are surrounded by much larger TWS anomalies due mainly to glacier melting in the mountains and groundwater extraction from Northwest India to the Bengal Basin in the order of 25-50 mm/yr [2]. Previous studies providing TWS changes based on GRACE products from different agencies and using various data-processing strategies thus show significant discrepancies over the TP (e.g., [12,13]). A dozen publications show TWS changes in the plateau interior with very different features, which includes the number, position, and magnitude of TWS anomalies (see

Table 1. Overview of terrestrial water storage trend rates inside the Tibetan Plateau from previous works.

Reference	GRACE Solution/Model	Time Span	Destriping Filter	Smoothing Filter	Shown in Figures	Anomalies
Matsuo and Heki [1]	CSR SH RL04	05/2003-04/2009	Cham P5M11	G400	3(b)	1, E/N
Erkan et al. [14]	CSR SH RL04	01/2003-12/2009	Duan P2M?	G300	8	2, E/N
Yi and Sun [15]	CSR SH RL05	01/2003-12/2012	Cham P4M6	G300	6(a)	1, W/N
Jiao et al. [16]	JPL Level-3	01/2003-12/2012	S&W P?M?	G300	1	1, E/N
Xiang et al. [2]	3 SH models	01/2003-12/2009	S&W P2M8	G340	A.3(a)	1, E/N
Pan et al. [17]	CSR SH RL05	04/2002-12/2014	Cham P4M6	G350	3	1, W/N
Guo et al. [18]	GRGS SH RL03	01/2003-12/2012	–	G300	3(d)	3, E/N
Feng et al. [12]	CSR SH RL05	04/2002-12/2014	–	G100	12(c)	5, W/N, E/N, M/N, M/S, E/S
Zou er al. [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 SH RL03	01/2003-12/2014	–	G250	6(d)	3, M/N, M/S, M/M
Jing et al. [13]	JPL_M	01/2003-12/2016	–	–	5(d)	1, M/N
Jing et al. [13]	CSR_M	01/2003-12/2016	–	–	5(e)	3, M/N
Loomis et al. [21]	GSFC_M	01/2003-07/2016	–	–	2(a)	1, M/N

Note: CSR = Center of Space Research; JPL = Jet Propulsion Laboratory; GFZ = German Research Centre for Geosciences; GRGS = Groupe de Recherche de Géodésie Spatiale; GSFC = Goddard Space Flight Center; RL0x = solution release number; 3 SH models = CSR, GFZ and JPL SH RL05 solutions; ‘Duan’, ‘S&W’ and ‘Cham’ 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.

). This may be partially attributed to different time periods or different averaging radius values of the Gaussian filters. However, their effects are far less than those from the selection of the GRACE data type and solution(s), and especially the chosen destriping filter(s). As Table 1 shows, different GRACE solutions using various data-processing filters, or even the same solution using different data-processing, lead to very controversial patterns of the GRACE-derived TWS changes in the interior of the TP: The number of mass increase anomalies varies between one and five, and they are located in the east or west by north or south or in the center of the plateau.

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.

2. GRACE Solutions and Destriping Filters

2.1. 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].

2.2. 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 ($\gamma ,p,A,K$) taking empirical values (0.1, 3, 30, 15),

$$w=max\left\{A{e}^{-\frac{{\left[\left(1-\gamma \right)m^p+\gamma l^p\right]}^{1/p}}{K}}+1, 5\right\}$$

(1)

However, for ‘S&W PxMy’ filtering, w is calculated using Equation (2), dependent on the order m of the coefficients, and the parameters $\left(A,K\right)$ taking the values (30, 10),

$$w=max\left\{A{e}^{-\frac{m}{K}}+1, 5\right\}$$

(2)

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.

3. Calculation of GRACE-Derived TWS Change and Aggregated Hydrology Signal

3.1. 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 $\theta$ and longitude $\varphi$ on $1.0°\times 1.0°$ grids covering our study area with

$$\Delta T^{TWS}\left(\theta ,\varphi \right)=\frac{a{\rho }_{ave}}{3{\rho }_w}{\sum }_{l=0}^{\infty }{\sum }_{m=0}^l\frac{2l+1}{1+k_l}\times \left(\Delta c_{lm}\cos m\varphi +\Delta s_{lm}\sin m\varphi \right){\tilde{P}}_{lm}\left(\cos \theta \right),$$

(3)

in which $\Delta c_{lm}$ and $\Delta s_{lm}$ are the changes in Stokes coefficients from GRACE gravity SH solutions or decomposed from three mascon models, ${\tilde{P}}_{lm}\left(\cos \theta \right)$ 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], ${\rho }_{ave}$ is the average density of the earth, ${\rho }_w$ is the water density and a is the radius of the Earth. The derived TWS changes $\Delta T^{TWS}\left(\theta ,\varphi \right)$ 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.

3.2. 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) in Xiang et al. [2], 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 $\Delta W^K$, and thus, the changes in spherical harmonic coefficients $\Delta c_{lm}^K$, $\Delta s_{lm}^K$ are given by decomposition of corresponding grid data $\Delta W^K$ designed referring to xiang et al. [2] in our surveyed area ($\Omega$) using

$$\{\begin{array}{c}\Delta c_{lm}^K=\frac{1}{4\pi }{\int }_{\Omega }\Delta W^K\left({\theta }^{\prime },{\varphi }^{\prime }\right){\tilde{P}}_{lm}\left(\cos {\theta }^{\prime }\right)\cos m{\varphi }^{\prime }d{\Omega }^{\prime }\\ \Delta s_{lm}^K=\frac{1}{4\pi }{\int }_{\Omega }\Delta W^K\left({\theta }^{\prime },{\varphi }^{\prime }\right){\tilde{P}}_{lm}\left(\cos {\theta }^{\prime }\right)\sin m{\varphi }^{\prime }d{\Omega }^{\prime }\end{array}$$

(4)

with $d{\Omega }^{\prime }=\sin {\theta }^{\prime }d{\theta }^{\prime }d{\varphi }^{\prime }$. 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.

4. Results and Discussion

4.1. Selection of the Solution-Filter Combination

4.1.1. Correlation Check

We calculated trends in the TWS changes from 2003 to 2009 in the TP and adjacent areas, based on the 130 solution-filter combinations outlined above and three mascon solutions. The best solution-filter combination(s) is determined by applying Pearson’s correlation coefficients (PCC-values), Root Mean Square Errors (RMSE-values) and Index Of Agreement values (IOA-values) between the GRACE-derived trends of TWS changes and the AHS during the same period in a selected region of the plateau interior (as shown in

Table 2. PCC-, RMSE- and IOA-values between the GRACE-derived TWS trends of selected GRACE solution centers (as in Figure 1a–l and Figure 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.

Destriping Filters	ITSG			CSR			COST-G			JPL			GFZ
	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA
Duan P2M8	0.39	9.79	0.60	0.35	10.16	0.57	0.37	9.91	0.58	0.29	11.52	0.54	0.41	10.85	0.59
Duan P3M8	0.56	7.96	0.69	0.54	8.08	0.67	0.55	7.95	0.68	0.44	9.65	0.61	0.57	9.07	0.66
Duan P4M8	0.66	5.81	0.79	0.67	5.95	0.77	0.65	6.10	0.76	0.63	6.22	0.75	0.63	7.36	0.73
S&W P2M8	0.66	5.89	0.76	0.62	6.38	0.71	0.64	6.33	0.71	0.56	6.90	0.68	0.61	7.52	0.70
S&W P3M8	0.68	5.72	0.79	0.67	5.94	0.76	0.66	6.05	0.76	0.63	6.16	0.75	0.64	7.21	0.74
S&W P4M8	0.25	7.79	0.52	0.21	8.07	0.49	0.22	8.01	0.49	0.20	7.90	0.47	0.28	8.70	0.53
Cham P2M8	0.02	14.57	0.39	−0.17	19.23	0.30	−0.02	14.92	0.37	0.11	15.67	0.41	−0.11	25.32	0.29
Cham P3M8	0.03	14.57	0.41	−0.15	19.58	0.31	−0.01	14.60	0.41	0.07	15.11	0.43	−0.10	19.40	0.36
Cham P4M8	0.03	14.63	0.41	−0.03	15.88	0.38	0.04	13.93	0.44	0.15	14.64	0.46	0.01	16.58	0.42
MLAs	0.23	8.14	0.50	0.45	7.91	0.56	-	-	-	0.37	8.45	0.54	0.24	10.89	0.50
DDK8	0.34	9.03	0.53	−0.05	12.17	0.39	-	-	-	0.28	11.31	0.50	0.07	12.70	0.44
not destriped	0.04	13.41	0.43	−0.13	17.88	0.32	0.08	11.62	0.43	0.30	16.48	0.48	−0.14	24.21	0.33
CSR_M						JPL_M					GSFC_M
PCC			RMSE		IOA	PCC		RMSE		IOA	PCC		RMSE		IOA
0.13			8.75		0.42	−0.04		11.48		0.37	−0.02		11.37		0.42

and

Table A1. Same as Table 2 but for GRACE-derived TWS trends derived from the other 7 different GRACE SH solutions.

Destriping Filters	Tongji			WHU			HUST			GRGS
	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA
Duan P2M8	0.27	9.09	0.53	0.28	8.92	0.55	0.26	9.85	0.53	0.38	9.76	0.58
Duan P3M8	0.47	7.20	0.65	0.49	6.87	0.68	0.44	8.05	0.62	0.51	8.57	0.64
Duan P4M8	0.58	5.69	0.71	0.55	5.77	0.70	0.62	5.83	0.75	0.57	6.79	0.70
S&W P2M8	0.56	5.86	0.66	0.54	5.89	0.67	0.51	6.49	0.65	0.61	6.78	0.68
S&W P3M8	0.57	5.77	0.71	0.54	5.90	0.70	0.53	6.38	0.69	0.63	6.48	0.72
S&W P4M8	0.09	8.15	0.43	0.05	8.06	0.39	0.17	8.06	0.47	0.15	8.30	0.45
Cham P2M8	0.27	11.35	0.52	0.12	12.50	0.41	−0.23	18.95	0.28	0.28	10.39	0.50
Cham P3M8	0.27	10.69	0.52	0.12	12.71	0.42	−0.21	18.12	0.30	0.27	10.75	0.50
Cham P4M8	0.33	9.67	0.55	0.12	13.46	0.42	−0.16	16.00	0.31	0.28	10.77	0.51
notdestriped	0.30	11.83	0.54	0.15	11.58	0.45	−0.20	16.97	0.28	0.21	10.68	0.47
Destriping Filters	IGG				AIUB				XISM
	PCC	RMSE	IOA		PCC	RMSE	IOA		PCC	RMSE	IOA
Duan P2M8	0.28	10.42	0.54		0.22	9.51	0.51		0.24	9.70	0.51
Duan P3M8	0.39	8.82	0.60		0.43	7.50	0.62		0.36	8.30	0.57
Duan P4M8	0.50	6.65	0.69		0.45	6.44	0.62		0.46	6.96	0.67
S&W P2M8	0.41	7.26	0.61		0.48	6.36	0.62		0.39	7.18	0.60
S&W P3M8	0.55	6.34	0.70		0.50	6.32	0.66		0.51	6.42	0.69
S&W P4M8	0.22	7.80	0.50		0.06	8.40	0.41		0.27	7.39	0.52
Cham P2M8	−0.04	16.33	0.35		0.07	21.65	0.33		−0.07	16.51	0.35
Cham P3M8	0.12	13.71	0.45		0.05	13.77	0.39		0.05	13.31	0.42
Cham P4M8	0.36	9.85	0.58		0.03	13.15	0.41		0.36	9.83	0.59
notdestriped	0.20	17.24	0.43		0.17	14.42	0.45		0.00	16.78	0.36

), 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.

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-derived trends of TWS changes based on five selected solution-filter combinations and AHS from 2002 to 2016 (Figure A6). Similarly, for exploring the effects of higher d/o Stokes coefficients (i.e., from 61 to 90) on the inversed TWS changes, we repeated the calculations mentioned above, but with Stokes coefficients truncated at maximum d/o 90 (Figure A8, Figure A9 and Figure A10).

Figure 1 and Figure 2 and Table 2 and Table A1 show that the results can vary drastically. The ITSG and CSR solution together with ‘Duan P4M8’, ‘S&W P2M8’ and ‘S&W P3M8’ result as best combinations for their analogous outstanding performance with PPC-values larger than 0.66, RMSE-values smaller than 5.95 and IOA-values larger than 0.76, except for CSR together with S&W P2M8. We show only ITSG results since results for other SH solutions are quite similar to those of the ITSG solution (Figure 1), and for similar reasons, we only show results derived using the ‘S&W P3M8’ filter (Figure 2). The variability of the results in the two figures will be discussed in Section 4.1.2 and Section 4.1.3.

The AHS from 2003 to 2009 in the TP and adjacent areas are shown in Figure 3f. Large negative signals are found in the Tien Shan (with magnitudes of ~35.0 mm/yr), the western/central/eastern Himalayas (~26.0 mm/yr, ~28.0 mm/yr, ~35.0 mm/yr) and the mountains southeast of the plateau (~28.0 mm/yr) mainly due to considerable ice melting. However, smaller signals consisting of three positive anomalies and one negative anomaly are found in the plateau interior. The western positive anomaly (~14.0 mm/yr) is a result of increased ice and soil moisture accumulation in the West Kunlun Mountains and slightly increasing lake water according to Figure 3a–c. The northeastern positive anomaly (~9.0 mm/yr) results from a combination of increased soil moisture and lake water mainly in Qinghai Lake. The positive anomaly (~11.5 mm/yr) to the south of the central TP is a comprehensive effect of lake water increase (~15.0 mm/yr) in the center of the plateau, the decrease in soil moisture to the southwest of the plateau, as well as a relatively small positive signal caused by spherical harmonic truncation from ICE. The prominent negative anomaly (~13.0 mm/yr) is mainly due to glacial melting in east Kunlun and the inner plateau [28] and permafrost degradation in the plateau, as well as the effects from soil moisture decreasing there.

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 Table 2 and Table A1, the PCC-, RMSE- and IOA-values are given for all tested solution-filter combinations and mascon models. The larger PCC- and IOA-values (in general PCC-values > 0.60, IOA-values > 0.70), and the smaller RMSE-values (in general RMSE-values < 6.40, except for GFZ) are obtained for combinations of destriping filters ‘Duan P4M8’, ‘S&W P2M8’ and ‘S&W P3M8’ with five GRACE SH solutions listed in Table 2. The agreement of GRACE-derived TWS trend rates with the AHS in the plateau interior is supported by visual comparison in Figure 1c–e and Figure 2a–e with Figure 3f. All other destriping filters generally result in much smaller PCC- and IOA-values, and much larger RMSE-values, implying that they are probably less effective, which we will verify in Section 4.3. All three mascon models also lead to much smaller PCC- and IOA-values, and much larger RMSE-values (in Figure 2f–h), showing that they disagree extremely with the AHS. Possible causes for this are discussed in Section 4.4.

4.1.2. 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., PCC-values $\ge$ 0.66, RMSE-values $\le$ 5.95 and IOA-values $\ge$ 0.76) would imply a reasonable selection of the destriping filter(s).

Figure 1a–f show results applying the destriping methods by Duan et al. [22] and Swenson and Wahr [3], respectively. Both apply moving windows of variable width. For the filters with x = 2, 3 and 4, the N-S stripes are effectively depressed in the GRACE-derived trends both inside and outside the plateau. However, better pattern and magnitude agreement of the results with the AHS (Figure 3f) is found only for ‘Duan P4M8’, ‘S&W P2M8’ and ‘S&W P3M8’, with larger PCC-values of 0.66, 0.66 and 0.68, and IOA-values of 0.79, 0.76 and 0.79, as well as smaller RMSE-values of 5.81, 5.89 and 5.72, respectively.

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 RMSE-values 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 RMSE-values 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.

4.1.3. Solution Dependence of the GRACE-Derived TWS Trends

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.

In Figure 2a–e, the patterns of the TWS change signals from all five GRACE SH solutions compare visually well with the AHS (Figure 3f). Their PCC-values are between 0.63 and 0.68, while IOA-values are between 0.74 and 0.79, and RMSE-values between 5.72 and 7.21. The best solutions are found to be ITSG and CSR. Both show the highest correlation and agreement but smallest bias with the AHS (also shown in 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 Section 4.3 and Section 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.

4.2. 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.

4.2.1. TWS Changes without Smoothing

Trend Rates

Based on the above pattern and magnitude agreements check, the visual comparison of the derived TWS changes with the AHS, and the visual check for N-S stripes, we determine the reliable TWS changes. Generally, the combinations of all twelve GRACE SH solutions with the three destriping filters ‘Duan P4M8’, ‘S&W P2M8’ and ‘S&W P3M8’ results in larger PCC- and IOA-values, and smaller RMSE-values than with other destriping filters, as seen in Table 2 and Table A1. PCC-values for these combinations are usually larger than 0.50, IOA-values larger than 0.60 and RMSE-values smaller than 7.0. Particularly, ITSG and CSR with the three destriping filters ‘Duan P4M8’, ‘S&W P2M8’ and ‘S&W P3M8’ (except for CSR together with ‘S&W P2M8’) can be used to determine TWS trends. Using these two solutions usually yields the largest PCC- and IOA-values as well as smallest RMSE-values among the twelve GRACE SH solutions independent of the chosen destriping filters. PCC-values for these best combinations are larger than 0.66, IOA-values larger than 0.76 and RMSE-values smaller than 5.95. Both the patterns and the amplitudes of TWS trend signals as shown in Figure 1c–e and Figure 2a,b are consistent with each other. For example, there are three positive and two negative anomalies inside the plateau in all the subfigures. The combination of ITSG SH solution and ‘S&W P3M8’ filter yields a magnitude of ~11.0 mm/yr for the western positive anomaly, ~14.0 mm/yr for the southern positive anomaly and ~17.5 mm/yr for the northeastern positive anomaly. The negative anomaly among the three positive anomalies is ~11.0 mm/yr. Another negative anomaly of ~43.0 mm/yr is mainly due to both increased ice melting and soil moisture reduction in the mountains southeast of the plateau.

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.

4.2.2. 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), 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.

For comparison we show trend rates (Figure 4e–j), annual amplitudes (Figure A4e–j), and annual phases (Figure A5e–j) from the ITSG SH solutions with filter ‘MLAs’, ‘DDK2’ and ‘Cham P4M8’, and the three mascon models, with the same additional Gaussian smoothing. They have different patterns and magnitudes than those determined with the best GRACE solution-filter combinations inside and outside the plateau (e.g., Figure 4a–d). For example, in the plateau interior, the single-anomalies of the smoothed trends (Figure 4e–j) have their centers in the center of the plateau, with an exception of 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).

4.3. 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 (Figure 1c–e and Figure 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. Figure 1c–e, Figure 2a,b and Figure 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), the PCC- and IOA-value are quite small, and RMSE-values are quite large. The positive anomaly inside the plateau after smoothing (Figure 4e) does not show its tail and has its center somewhat in the center of the plateau, which resembles the benchmark result (Figure 4k). Of course, this is not supported by the smoothed AHS (Figure 4l).

Larger PCC- and IOA-values, and smaller RMSE-values are achieved if GRACE SH solutions are combined with the best-performing filters ‘Duan P4M8’, ‘S&W P2M8’ and ‘S&W P3M8’ (Figure 1c–e and Figure 2a–e, Table 2 and Table A1), whereas largest values yield for ITSG and CSR (except ‘S&W P2M8’) which also match the AHS very well before and after smoothing with a Gaussian filter. When the time span is extended from ending in December 2009 to ending in December 2016, a similar performance is found for the above filters and GRACE solutions. However, the PCC-, RMSE- and IOA-values shown in 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.

4.4. 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 ($\sigma$) for the AHS (Figure A7e) using Equation (A4) and a square root function. Further, four extreme cases of AHS (Figure A7f–i) considering $\pm 1\sigma$ and $\pm 2\sigma$ 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 $\sigma$ 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 $\pm 1\sigma$’ and ‘AHS $\pm 2\sigma$’, 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

Table A2. Same as Table 2 but for the ‘AHS+$1\sigma$’ case instead of aggregated hydrology signals (AHS).

Destriping Filters	ITSG			CSR			COST-G			JPL			GFZ
	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA
Duan P2M8	0.45	8.43	0.66	0.43	8.61	0.64	0.43	8.40	0.65	0.37	9.91	0.60	0.46	9.29	0.64
Duan P3M8	0.58	6.90	0.74	0.58	6.72	0.74	0.58	6.64	0.74	0.50	8.26	0.67	0.58	7.71	0.72
Duan P4M8	0.64	5.48	0.79	0.66	5.13	0.80	0.62	5.41	0.78	0.64	5.26	0.80	0.60	6.46	0.76
S&W P2M8	0.66	5.04	0.80	0.65	5.12	0.78	0.65	5.16	0.78	0.62	5.44	0.77	0.60	6.24	0.75
S&W P3M8	0.64	5.41	0.79	0.65	5.13	0.80	0.63	5.36	0.78	0.64	5.25	0.79	0.59	6.36	0.76
S&W P4M8	0.15	7.70	0.46	0.13	7.64	0.46	0.13	7.67	0.45	0.14	7.45	0.44	0.18	8.13	0.49
Cham P2M8	0.04	13.83	0.37	−0.18	18.65	0.25	0.00	13.98	0.35	0.19	14.31	0.44	−0.11	24.74	0.25
Cham P3M8	0.05	13.78	0.39	−0.17	19.03	0.26	0.01	13.65	0.38	0.11	13.89	0.43	−0.11	18.65	0.31
Cham P4M8	0.05	13.85	0.39	−0.02	15.03	0.35	0.07	12.91	0.43	0.21	13.31	0.48	0.00	15.66	0.39
MLAs	0.21	7.71	0.52	0.48	6.42	0.64	-	-	-	0.37	7.14	0.59	0.18	10.09	0.48
DDK8	0.40	7.49	0.62	−0.05	11.06	0.37	-	-	-	0.35	9.42	0.57	0.08	11.22	0.45
not destriped	0.01	12.73	0.40	−0.14	17.27	0.27	0.09	10.32	0.44	0.39	14.96	0.52	−0.17	23.63	0.28
CSR_M						JPL_M					GSFC_M
PCC			RMSE		IOA	PCC		RMSE		IOA	PCC		RMSE		IOA
0.15			7.58		0.46	0.01		9.88		0.40	−0.01		9.85		0.43

,

Table A3. Same as Table A2 but for the ‘AHS-$1\sigma$’ case.

Destriping Filters	ITSG			CSR			COST-G			JPL			GFZ
	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA
Duan P2M8	0.32	11.71	0.53	0.27	12.20	0.50	0.30	11.93	0.51	0.21	13.56	0.48	0.36	12.87	0.53
Duan P3M8	0.52	9.79	0.62	0.49	10.10	0.58	0.52	9.95	0.59	0.38	11.60	0.54	0.55	11.02	0.60
Duan P4M8	0.66	7.34	0.72	0.66	7.81	0.68	0.65	7.85	0.68	0.60	8.14	0.66	0.64	9.12	0.67
S&W P2M8	0.64	7.78	0.67	0.58	8.47	0.61	0.61	8.37	0.61	0.50	9.06	0.57	0.60	9.52	0.62
S&W P3M8	0.69	7.26	0.73	0.67	7.79	0.68	0.67	7.81	0.68	0.61	8.05	0.66	0.66	8.95	0.68
S&W P4M8	0.32	8.88	0.52	0.27	9.41	0.47	0.30	9.27	0.48	0.25	9.26	0.45	0.36	10.09	0.52
Cham P2M8	0.00	15.80	0.40	−0.15	20.20	0.35	−0.05	16.32	0.38	0.03	17.40	0.39	−0.10	26.21	0.33
Cham P3M8	0.01	15.86	0.42	−0.14	20.52	0.36	−0.03	16.01	0.42	0.02	16.74	0.42	−0.08	20.52	0.39
Cham P4M8	0.02	15.91	0.43	−0.05	17.18	0.40	0.02	15.42	0.44	0.08	16.38	0.45	0.01	17.92	0.44
MLAs	0.24	9.48	0.46	0.42	10.02	0.47	-	-	-	0.35	10.40	0.47	0.28	12.33	0.50
DDK8	0.27	11.11	0.44	−0.05	13.79	0.39	-	-	-	0.20	13.54	0.44	0.06	14.61	0.42
not destriped	0.07	14.63	0.45	−0.12	18.91	0.35	0.08	13.42	0.41	0.21	18.32	0.45	−0.10	25.10	0.38
CSR_M						JPL_M					GSFC_M
PCC			RMSE		IOA	PCC		RMSE		IOA	PCC		RMSE		IOA
0.11			10.58		0.38	−0.09		13.51		0.34	−0.03		13.34		0.40

,

Table A4. Same as Table A2 but for the ‘AHS+$2\sigma$’ case.

Destriping Filters	ITSG			CSR			COST-G			JPL			GFZ
	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA
Duan P2M8	0.49	7.93	0.69	0.49	7.85	0.68	0.48	7.72	0.68	0.44	8.95	0.64	0.49	8.46	0.67
Duan P3M8	0.58	6.95	0.74	0.59	6.44	0.76	0.58	6.44	0.75	0.53	7.75	0.70	0.56	7.29	0.73
Duan P4M8	0.58	6.55	0.72	0.61	5.81	0.75	0.57	6.15	0.72	0.64	5.77	0.76	0.54	6.76	0.73
S&W P2M8	0.64	5.72	0.76	0.65	5.32	0.77	0.63	5.45	0.75	0.65	5.31	0.78	0.57	6.16	0.75
S&W P3M8	0.57	6.51	0.71	0.60	5.83	0.74	0.56	6.12	0.72	0.62	5.81	0.75	0.52	6.74	0.72
S&W P4M8	0.05	8.62	0.37	0.04	8.24	0.38	0.03	8.37	0.37	0.07	8.07	0.36	0.08	8.55	0.41
Cham P2M8	0.06	13.67	0.35	−0.18	18.50	0.19	0.03	13.59	0.33	0.27	13.44	0.47	−0.11	24.49	0.22
Cham P3M8	0.07	13.56	0.37	−0.17	18.91	0.19	0.03	13.27	0.33	0.16	13.19	0.42	−0.13	18.33	0.26
Cham P4M8	0.06	13.63	0.36	0.00	14.69	0.31	0.09	12.48	0.39	0.27	12.50	0.50	−0.01	15.22	0.34
MLAs	0.18	8.30	0.49	0.48	6.04	0.67	-	-	-	0.36	6.87	0.60	0.11	10.07	0.45
DDK8	0.45	6.87	0.68	−0.04	10.64	0.32	-	-	-	0.41	8.14	0.64	0.08	10.34	0.44
not destriped	−0.01	12.68	0.36	−0.15	17.12	0.23	0.09	9.73	0.42	0.47	13.87	0.56	−0.21	23.39	0.23
CSR_M						JPL_M					GSFC_M
PCC			RMSE		IOA	PCC		RMSE		IOA	PCC		RMSE		IOA
0.16			7.41		0.45	0.06		8.94		0.41	0.00		9.02		0.41

and

Table A5. Same as Table A2 but for the ‘AHS-2$\sigma$’ case.

Destriping Filters	ITSG			CSR			COST-G			JPL			GFZ
	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA	PCC	RMSE	IOA
Duan P2M8	0.26	13.96	0.48	0.20	14.53	0.45	0.23	14.25	0.46	0.13	15.87	0.44	0.31	15.17	0.48
Duan P3M8	0.47	12.03	0.55	0.43	12.46	0.51	0.47	12.29	0.52	0.32	13.87	0.49	0.51	13.32	0.54
Duan P4M8	0.65	9.52	0.64	0.64	10.16	0.59	0.64	10.13	0.59	0.55	10.50	0.56	0.64	11.34	0.61
S&W P2M8	0.60	10.14	0.57	0.53	10.92	0.51	0.57	10.80	0.52	0.43	11.54	0.47	0.57	11.89	0.55
S&W P3M8	0.68	9.44	0.65	0.65	10.13	0.59	0.67	10.09	0.60	0.57	10.40	0.57	0.66	11.16	0.61
S&W P4M8	0.38	10.65	0.49	0.32	11.33	0.43	0.36	11.15	0.45	0.29	11.21	0.41	0.42	12.01	0.50
Cham P2M8	-0.02	17.43	0.41	−0.14	21.52	0.38	−0.07	18.08	0.38	−0.04	19.40	0.37	−0.09	27.37	0.36
Cham P3M8	-0.01	17.52	0.42	−0.12	21.80	0.39	−0.04	17.78	0.42	−0.02	18.68	0.41	−0.06	21.97	0.42
Cham P4M8	0.00	17.57	0.43	−0.06	18.83	0.41	0.00	17.27	0.44	0.03	18.40	0.43	0.02	19.60	0.45
MLAs	0.24	11.39	0.42	0.38	12.44	0.40	-	-	-	0.33	12.71	0.41	0.32	14.21	0.48
DDK8	0.21	13.49	0.38	−0.06	15.78	0.38	-	-	-	0.14	15.99	0.39	0.05	16.79	0.39
not destriped	0.08	16.27	0.45	−0.11	20.30	0.38	0.07	15.54	0.38	0.13	20.41	0.43	−0.06	26.28	0.41
CSR_M						JPL_M					GSFC_M
PCC			RMSE		IOA	PCC		RMSE		IOA	PCC		RMSE		IOA
0.10			12.81		0.33	−0.12		15.81		0.31	−0.03		15.60		0.37

. In Table A2 and Table A3, the PCC-, RMSE- and IOA-values are for the cases of ‘AHS+$1\sigma$’ and ‘AHS-$1\sigma$’, respectively. Similar to Table 2, the larger PCC- and IOA-values and smaller RMSE-values are still obtained for combinations of destriping filters ‘Duan P4M8’, ‘S&W P2M8’ and ‘S&W P3M8’ separated by two dashed lines in the Table A2 and Table A3, while their combinations with ITSG and CSR SH models show the highest correlation and agreement. All other destriping filters and mascon solutions generally result in much smaller PCC- and IOA-values, and much larger RMSE-values. This confirms the selected solutions and destriping filters in our study and verifies the effectiveness of our designed AHS criterion.

The $\pm 1\sigma$ variation of AHS leads to only small changes in PCC-values, but alters both RMSE- and IOA-values remarkably clearly (Table A2 and Table 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 $\pm 1\sigma$ 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\sigma$’ case (Table A2) yields smaller RMSE-values and better IOA agreement than the AHS, while the ‘AHS-$1\sigma$’ case (Table A3) results in larger RMSE-values and worse IOA agreement than the AHS. That means the ‘AHS+$1\sigma$’ 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\sigma$’ case. When adding or subtracting $2\sigma$ from AHS, the ‘AHS-$2\sigma$’ case (Table A5) performs even worse than the ‘AHS-$1\sigma$’ case, while in the ‘AHS+$2\sigma$’ 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, unknown 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 Section 3 and Section 4 but using Stokes coefficients maximum d/o up to 90, and show corresponding results in Figure A8, Figure A9 and Figure 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,f) or different SH solutions (Figure A9a–e). All the results show poor pattern and magnitude agreements with the AHS (Figure A10f). Thus, rather low PCC- and IOA-values and rather high RMSE-values are expected, which implies heavy N-S stripe noises in the TWS change trends. In other words, current destriping filters cannot depress the N-S stripes effectively for Stokes coefficients of higher d/o (i.e., 60–90). The destriping filters Cham PxM8 (x = 2,3,4) showed the worst effectiveness in depressing N-S stripes, because the related TWS trends (Figure A8g–i) resemble the pattern of unfiltered results (Figure A8l), showing nearly no valid signal at all.

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, Figure 1, Figure 2 and Figure 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 Figure 1c–e, Figure 2a,b and Figure 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.

5. 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.