A Combined Use of TSVD and Tikhonov Regularization for Mass Flux Solution in Tibetan Plateau

Abstract

Limited by the Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On (GRACE-FO) measurement principle and sensors, the spatial resolution of mass flux solutions is about 2–3° in mid-latitudes at monthly intervals. To retrieve a mass flux solution in the Tibetan Plateau (TP) with better visual spatial resolution, we combined truncated singular value decomposition (TSVD) and Tikhonov regularization to solve for a mascon modeling. The monthly mass flux parameters resolved at 1° are smoothed to about 2° by truncating the eigen-spectrum of the normal equation (i.e., using the TSVD approach), and then Tikhonov regularization is applied to the truncated normal equation. As a result, the terms beyond the native resolution of GRACE/GRACE-FO data are truncated, and the errors in higher degree and order components are dampened by Tikhonov regularization. In terms of root mean squared errors, the improvements are 27.2% and 12.7% for the combined method over TSVD and Tikhonov regularization, respectively. We confirm a decreasing secular trend with −5.6 ± 4.2 Gt/year for the entire TP and provide maps with 1° resolution from April 2002 to April 2019, generated with the combined TSVD and Tikhonov regularization method.

1. Introduction

Since the launch of Gravity Recovery and Climate Experiment Follow-On (GRACE-FO) in May 2018 [1], the number of monthly products of time-variable gravity fields has grown, which provides the opportunity to further investigate large-scale secular and seasonal geophysical signals [2,3,4]. Known as the Third Pole, the Tibetan Plateau holds the largest number of glaciers, acting as an important contributor to major rivers in Asia, such as the Yangtze and the Ganges [5,6]. Due to the geographic complexity, remote sensing analyses, such as employing GRACE data, the digital elevation model (DEM) differencing approach, or converting height changes from laser altimetry (e.g., from Ice, Cloud, and Elevation Satellite, ICESat) to mass variations, are used in estimating the mass variation of the Tibetan Plateau (TP). For example, Brun et al. calculated the glacial mass loss at a trend of −16.3 ± 3.5 Gt/year between 2000 and 2016 in High Mountain Asia with a DEM [7] and Jacob et al. reported a less negative rate of −4 ± 20 Gt/year based on GRACE data [8]. The differences in mass trends that we find stress the need for a better understanding of the mass variations in the TP with high spatial resolution.

Due to the weak sensitivity of its observations in the east-west direction, measurement errors, and the ill-posed inversion problem, unconstrained GRACE-based spherical harmonic solutions show unreasonable north-south stripes [9]. Several post-processing techniques, such as applying Gaussian filters [10,11], decorrelation filters [12], and DDK filters [13,14], have been developed to constrain the errors in high degree and order. However, they tend to cause an attenuation of real geophysical signals [15]. Another option of recovering mass variations from GRACE data while avoiding the use of post-processing techniques is to use mass concentration approaches applied to GRACE inter-satellite range-rate measurements [16], or range-acceleration data [2,17], or spherical harmonic coefficients [8,18]. However, the normal equations that are formed within the mass concentration approach are typically ill-conditioned, and either truncated singular value decomposition (TSVD) or Tikhonov regularization is required to derive a stable solution. Moreover, the native spatial resolution of the mass flux solutions derived from GRACE data is about 2–3° (63,000 km²) at an error level of 2 cm in terms of equivalent water height (EWH) at monthly intervals [19].

It is well known that TSVD and Tikhonov regularization modify the eigen-spectrum of the normal equations in different ways, and combining the two methods adds one additional degree of freedom. More degrees of freedom for turning the normal equations are used in the multi-parameter regularization approach [20]. Although we cannot claim that the combination will automatically enable results superior to carefully tuned TSVD or Tikhonov regularization solutions, we hypothesize that it provides more flexibility when the gravity field is parameterized through mascons. The goal of the study is to derive an improved solution via the combined use of TSVD and Tikhonov regularization with respect to individual use of them: truncating the eigen-spectrum of the normal equation beyond the native spatial resolution of GRACE/GRACE-FO by TSVD, and then dampening the errors in higher degree and order components by Tikhonov regularization. The rest of this paper is organized as follows: Section 2 presents the combined TSVD and Tikhonov method after introducing the mascon modeling; Section 3 gives the details of total mass variations and spatial distribution of mass flux solution at 1° resolution in the TP and the conclusions are drawn in Section 4.

2. Data and Methods

2.1. GRACE Data

Since the launch of the GRACE mission in 2002, the GRACE Science Data System (SDS) institutions, namely, the Center for Space Research (CSR), the Jet Propulsion Laboratory (JPL), and the GeoForschungsZentrum (GFZ), have been updating gravity monthly solutions during release 01 to 06. Since the spatial, temporal, and spectral differences of these SDS products were shown to be within a certain scatter [21], they differ little among SDS products. Here we chose to employ the up-to-date CSR release 06 products [22] from April 2002 to April 2019 (172 months), including GRACE-FO data from June 2018 to April 2019. We follow the usual post-processing steps; since the GRACE monthly solutions refer to the center of mass, the offset between the center of mass and the center of the frame of the Earth, represented by the degree-1 coefficients of spherical harmonics [23], must be added back when the mass flux solutions are defined in the reference frame. We use degree-1 coefficients derived from satellite laser ranging (SLR) observations [24] to correct this offset. C₂₀ and C₃₀ terms, which are poorly observed by GRACE and GRACE-FO, in particular, with only a single accelerometer, are replaced with SLR observations by the technical note 14 (TN-14) document recommended by Loomis et al. [25].

The response of the solid Earth to the glacial loading and unloading leads to glacial isostatic adjustment (GIA), and this signal should be removed from surface mass variations to enable an interpretation in terms of hydrology. We corrected our mass flux solution with the three-dimensional (3D) GIA model [26] and removed the GIA signal by 3.6 Gt/year from mass flux solution in the TP (Figure 1). With various ice sheet models and viscoelastic Earth models, other GIA models ranging from 0 mm/year to 2.5 mm/year account for mass variations of 0.2 Gt/year to 4.9 Gt/year in the TP, which indicate large uncertainties due to the complexity of the geodynamic process [27].

2.2. Mascon Modeling

In mascon approaches it is common to first transform the spherical harmonic coefficients into a set of pseudo-observables over the region of interest; typically, gravity disturbances at satellite altitude on a regular grid. When we account for the Earth’s elastic yielding, the radial gravity disturbance can be expressed in terms of fully normalized Stokes coefficients as follows [23],

$$\delta g=-\frac{\mathrm{GM}}{r^2}{\displaystyle \sum }_{l=2}^{l_{max}}{\left(\frac{a}{r}\right)}^l\frac{(l+1)}{(1+k_l^{\prime })}{\displaystyle \sum }_{m=0}^l{P}_{lm}\left(\sin \phi \right)(\mathsf{\Delta }{\mathrm{C}}_{lm}\cos m\lambda +\mathsf{\Delta }{\mathrm{S}}_{lm}\sin m\lambda )$$

(1)

where r, λ, and φ are the spherical coordinates of the pseudo observations (i.e., the radius (r = a + 500 km), longitude, and latitude, respectively); G indicates the gravitational constant; a is the mean radius of the Earth; l and m indicate the degree and order of the gravity field model; the maximum degree l_max is 60; $k_l^{\prime }$ (provided by Wahr et al. [11]) is the load Love number of degree l; P_lm stands for the fully normalized associated Legendre functions; $\mathsf{\Delta }{\mathrm{C}}_{lm}$ and $\mathsf{\Delta }{\mathrm{S}}_{lm}$ are the Stokes coefficients after the mean gravity field is removed.

According to Newton’s law of gravity, assuming n pseudo observations at satellite attitude caused by t mass-points in the study area, the total radial gravitational disturbance at each location j is as follows [18],

$$\delta g_i=\mathrm{G}{\displaystyle \sum }_{j=1}^t{\delta m}_j\frac{(r-a\cos {\psi }_{i,j})}{{(a^2+r^2-2ar\cos {\psi }_{i,j})}^{3/2}},i=1,2,\dots ,n,j=1,2,\dots ,t$$

(2)

where${\delta m}_j$ (j = 1,2,…,t) is the unknown mass variations to be estimated and ${\psi }_{i,j}$ is the spherical distance.

Together with Equations (1) and (2), the radial gravitational disturbance in the pseudo observation points at satellite altitude connects the spherical harmonic coefficients and the mass variation due to the ground mass-points, and we reformulated the observation equation here as

$$\mathsf{y}=\mathsf{Ax}+\mathsf{e},{\mathsf{A}}_{i,j}=\mathrm{G}\frac{(r-a\cos {\psi }_{i,j})}{{(a^2+r^2-2ar\cos {\psi }_{i,j})}^{3/2}},i=1,2,\dots ,n,j=1,2,\dots ,t,$$

(3)

where y is an n-vector of pseudo observations of the radial gravitational disturbance calculated by spherical harmonic coefficients; A is an n × t design matrix (n > t) and denotes an overdetermined system of the equation; x is a t-vector of unknown ground mass-points to be estimated (i.e., ${\delta m}_j$); e denotes the n-vector of random errors with zero mean and variance of unit weight ${\sigma }_0^2$. Via the law of error propagation, the weight matrix P is computed with P = (BDB^T)⁻¹, in which D is the covariance matrix from the CSR release 06 products; B is the coefficient matrix of projecting the spherical harmonics to the pseudo observation vector with its ith row elements corresponding to spherical harmonics $\mathsf{\Delta }{\mathrm{C}}_{lm}$ and$\mathsf{\Delta }{\mathrm{S}}_{lm}$which are written as [28],

$$\mathsf{B}(i,j_{\mathsf{\Delta }{C}_{lm}})=-\frac{\mathrm{GM}}{r^2}{(\frac{a}{r})}^l\frac{l+1}{{1+k}_l^{\prime }}{P}_{lm}(\sin {\phi }_i){\cos m\lambda }_i,$$

(4)

$$\mathsf{B}(i,j_{\mathsf{\Delta }{S}_{lm}})=-\frac{\mathrm{GM}}{r^2}{(\frac{a}{r})}^l\frac{l+1}{{1+k}_l^{\prime }}{P}_{lm}\left(\sin {\phi }_i\right){\sin m\lambda }_i,$$

(5)

where $({\lambda }_i,{\phi }_i,r)$ is the spherical coordinate of the ith pseudo observation; $j_{\mathsf{\Delta }{C}_{lm}}$ and $j_{\mathsf{\Delta }{S}_{lm}}$ are the column numbers corresponding to spherical harmonics$\mathsf{\Delta }{\mathrm{C}}_{lm}$ and$\mathsf{\Delta }{\mathrm{S}}_{lm}$, respectively. According to Ran et al. [29], introducing the weight matrix $\mathsf{P}$can improve the quality of estimated mass variations. Tikhonov regularization, with the L-curve method applied to choose the regularization parameter, is widely used in GRACE-derived regional and global mascon solutions [18,30].

Within our study area of 63°E-105°E, 25°N-46°N, we constructed a mascon grid at ground level, at 1° resolution with 946 mascons (Figure 1). In the figure, the center of each mascon (1° × 1°) stands for the coordinate of the corresponding ground mascon to be estimated. The mass flux solution for the entire TP consists of the individual mass variations of each mascon inside the boundary of the TP (black line in Figure 1). The pseudo observations of radial gravitational disturbance cover an area of 60°E-108°E, 22°N-49°N of 0.8° resolution and 1° resolution, respectively.

2.3. Combined Use of TSVD and Tikhonov Regularization

To recover the solution x, we applied the least-squares adjustment to minimize the square sum of e. The design matrix A (n × t) can be expressed by singular value decomposition (SVD) as

$$\mathsf{A}={\mathsf{USV}}^{\mathrm{T}},$$

(6)

in which the columns of U = [ u₁, u₂, …, u_n] are the left singular vectors of A; the diagonal matrix S = diag[s₁, s₂, …, s_t] is the singular value of A; the columns of V = [v₁, v₂, …, v_t] are the right singular vectors of A. The least-squares solution ${\mathsf{x}}_{\mathrm{LS}}$ is,

$${\mathsf{x}}_{\mathrm{LS}}={\displaystyle \sum }_{j=1}^{\mathrm{t}}\frac{{\mathsf{u}}_j^T\mathsf{y}}{{\mathsf{s}}_j}{\mathsf{v}}_j.$$

(7)

As the index increases, the singular value ${\mathsf{s}}_j^{-1}$ of A increases as well and small perturbations in the observations can cause significant perturbations in the solution. Due to the sampling and geometry of pseudo observation, we found a condition number for the design matrix A of 6.9·10¹⁷ (with the geodesic grid of 1° × 1°), indicating an ill-conditioned normal equation N = A^TPA. There are two regularization methods presented: TSVD and Tikhonov regularization. By eliminating small singular values, the TSVD solution x_k is in the form of [31],

$${\mathsf{x}}_{\mathrm{k}}={\displaystyle \sum }_{j=1}^k\frac{{\mathsf{u}}_j^T\mathsf{y}}{{\mathsf{s}}_j}{\mathsf{v}}_j,$$

(8)

where 1 < k < t. If the truncation parameter k is chosen too large, the condition number of A remains large; however, a small k value leads to losing a large part of the signals. Tikhonov regularization is commonly applied to stabilize ill-posed problems, it is well-known that the solution minimizes the cost function

$$\mathsf{\Phi }\left(\mathsf{x}\right)={\left(\mathsf{Ax}-\mathsf{y}\right)}^{\mathrm{T}}\mathsf{P}\left(\mathsf{Ax}-\mathsf{y}\right)+\mu {\mathsf{x}}^{\mathrm{T}}\mathsf{Rx},$$

(9)

where $\mu$ indicates the regularization parameter ($\mu$ > 0) and R is the regularization matrix. The solution to Equation (9) is expressed as,

$${\mathsf{x}}_{\mu }={\left({\mathsf{A}}^{\mathrm{T}}\mathsf{PA}+\mu \mathsf{R}\right)}^{-1}{\mathsf{A}}^{\mathrm{T}}\mathsf{Py},$$

(10)

which shows that Tikhonov regularization dampens every component of the design matrix.

Here we construct a new TSVD and Tikhonov regularization by applying Tikhonov regularization to the truncated normal equation. As the regularization matrix R is chosen as an identity matrix I and the design matrix A is truncated to the first k singular values, the solution ${\mathsf{x}}_{\mu k}$is,

$${\mathsf{x}}_{\mu k}=\mathsf{V}{\left({\mathsf{S}}_{\mathrm{k}}{}^{\mathrm{T}}{\mathsf{PS}}_{\mathrm{k}}+\mu \mathsf{I}\right)}^{-1}{\mathsf{S}}_{\mathrm{k}}{}^{\mathrm{T}}{\mathsf{U}}^{\mathrm{T}}\mathsf{Py},{\mathsf{S}}_{\mathrm{k}}={\mathrm{diag}[\mathrm{s}}_1{,\mathrm{s}}_2{,\dots ,\mathrm{s}}_k,0,\dots ,0].$$

(11)

Consisting of the variance-covariance matrix of $x_{\mu k}$ and the bias correction of both Tikhonov regularization and the truncated part [32,33], the corresponding mean squared error (MSE) matrix is,

$$\mathrm{MSE}\left(x_{\mu k}\right)={\mathsf{\sigma }}_0^2{\mathsf{Q}}_{\mu k}{\mathsf{VS}}_{\mathrm{k}}{}^{\mathrm{T}}{\mathsf{PS}}_{\mathrm{k}}{\mathsf{V}}^{\mathrm{T}}{\mathsf{Q}}_{\mu k}+{\mu }^2{\mathsf{Q}}_{\mu k}\tilde{\mathsf{x}}{\tilde{\mathsf{x}}}^{\mathrm{T}}{\mathsf{Q}}_{\mu k}+{\tilde{\mathsf{x}}}^{\mathrm{T}}{\mathsf{V}}_{\mathrm{k}}{\mathsf{V}}_{\mathrm{K}}^{\mathrm{T}}\tilde{\mathsf{x}},{\mathsf{V}}_{\mathrm{k}}=[0,\dots ,0,{\mathsf{v}}_{\mathrm{k}+1},\dots ,{\mathsf{v}}_{\mathrm{k}}]$$

(12)

in which $\tilde{\mathsf{x}}$ is the true value, and ${\mathsf{Q}}_{\mu k}=\mathsf{V}{\left({\mathsf{S}}_{\mathrm{k}}{}^{\mathrm{T}}{\mathsf{PS}}_{\mathrm{k}}+\mu \mathsf{I}\right)}^{-1}{\mathsf{V}}^{\mathrm{T}}$. The variance matrix of $x_{\mu k}$is ${\mathsf{\sigma }}_0^2{\mathsf{Q}}_{\mu k}{\mathsf{VS}}_{\mathrm{k}}{}^{\mathrm{T}}{\mathsf{PS}}_{\mathrm{k}}{\mathsf{V}}^{\mathrm{T}}{\mathsf{Q}}_{\mu k}$. The ${\mu }^2{\mathsf{Q}}_{\mu k}\tilde{\mathsf{x}}{\tilde{\mathsf{x}}}^{\mathrm{T}}{\mathsf{Q}}_{\mu k}$ and ${\tilde{\mathsf{x}}}^{\mathrm{T}}{\mathsf{V}}_{\mathrm{k}}{\mathsf{V}}_{\mathrm{K}}^{\mathrm{T}}\tilde{\mathsf{x}},$correspond to bias impacts of Tikhonov regularization and TSVD, respectively. For more details, please refer to [32,33]. To estimate the unknown ${\sigma }_0$, the following equation is introduced [34],

$${\hat{\sigma }}_0^2=\frac{{\hat{\mathsf{e}}}^{\mathrm{T}}\mathsf{P}\hat{\mathsf{e}}-{\mu }^2{\tilde{\mathsf{x}}}^{\mathrm{T}}({\mathsf{Q}}_{\mu k}-\mu {\mathsf{Q}}_{\mu k}^2)\tilde{\mathsf{x}}}{(n-t)+{\mu }^2\mathrm{tr}({\mathsf{Q}}_{\mu k}^2)}.$$

(13)

in which $\hat{\mathsf{e}}$ is the estimation of the error vector. Since the true value $\tilde{\mathsf{x}}$ remains unknown, we replaced it with the estimation $x_{\mu k}$ by the combined method.

Comparing Equation (10) with Equation (11), it is obvious that this combined Tikhonov regularization adds a small regularization parameter $\mu$to the truncated ${\mathrm{s}}_j$($1<j\le k)$ only. The choice of the truncated number k controls the amount of useful information derived from the normal equation and the Tikhonov regularization further dampens the errors of higher degree and order.

When considering the location of the inflection point of the eigenvalue curve for the case of 1° resolution, we selected the truncation parameter k as 270 (about 2° resolution) for the combined method; we then applied Tikhonov regularization to the truncated normal equation and derived the optimal regularization parameter $\mu$ based on the criterion of minimizing the traced MSE (for details, see [33,34]); finally, we calculated the mass flux solution with Equation (11). Other methods to determine the regularization parameter such as the discrepancy principle [35], the generalized cross-validation method [36], or the L-curve criterion [18,30] could be investigated.

Figure 2a shows the eigenvalues of the normal equation with the 1° resolution with different methods in January 2004. It should be noted that the choice of regularization parameter and truncation value may bias the mass variations derived from the combined method, TSVD, or Tikhonov regularization. To compare the combined solution with other solutions, we did not truncate the normal equation with Tikhonov regularization only (Figure 2a, blue dotted line) and derived the optimal regularization parameter $\mu$ based on the criterion of minimizing the traced MSE. Note that the regularization parameters of the combined method and Tikhonov regularization are separately calculated for each month. The truncating parameter k of TSVD solutions is chosen according to the MSE criterion in each month [30]. We further compared the condition numbers by the combined method, and Tikhonov regularization, respectively in Figure 2b. The condition numbers by the combined methods were reduced slightly more compared to those by Tikhonov regularization during 2002 to 2019. Figure 2c presents how traced MSE of the combined method varies with the regularization parameter in January 2004. The regularization parameter of 1.6 × 10⁻²⁹ is selected by the minimum traced MSE.

Limited by the GRACE data, the spatial resolution of the solution is about 2°. Therefore all 1° grids are overparameterized; due to rounding error, the eigenvalues lower than 10⁻³³ drop very slowly, and the correspondent eigenvectors form linearly dependent parameters. Truncating them to about a 2° resolution is better than direct parameterizing with 2° grids, since the combinations with larger eigenvalues are left. However, if we truncate more terms to get a stable solution, we may lose the signals corresponding to small eigenvalues. The terms for the eigenvalues larger than 10⁻³³ are corresponding to the spatial resolution of GRACE data; Tikhonov regularization is therefore used to derive a stable solution.

2.4. Leakage Correction

The leakage errors in GRACE observations increase proportionally with the strength of the signal source and attenuate quickly with the increasing distance between the source and the mass points. Similar to other studies [28], we expanded the study area by 3°, covering the area of 60°E-108°E, 22°N-49°N, and calculated these additional mass variations outside the study area with mascon modeling as well. To eliminate the leakage error, we removed the radial gravitational disturbance derived from these nearby additional mass variations from our pseudo observations.

3. Results and Discussion

3.1. MSE Roots

Here we introduced MSE, which is one of the most common measures to describe how methods perform in parameter estimation. Estimators (i.e., mascons in our case) with a smaller MSE, indicate that they are expected to have a closer distance to the true value. To compare the performance of each method in the parameter space, we calculated the MSE roots of 946 mascons in the TP using the combined method with Equation (12), those by TSVD according to Xu [32], and those by Tikhonov regularization according to Shen and colleagues [33]. Figure 3 shows the MSE roots in EWH (cm) of all mascons for each month. Except for a few months, the combined method has a better performance than Tikhonov regularization. The maximum MSE roots of 946 mascons were successfully reduced from 8.31 cm in the case of TSVD to 6.74 cm in the case of the combined method. The maximum MSE root by Tikhonov regularization is equal to 7.56 cm. The MSE roots of 3.08 cm by the combined method is smaller than those of 4.23 cm and 3.53 cm by TSVD and Tikhonov regularization, respectively (

Table 1. MSE roots for three methods in EWH (cm).

MSE Roots	Maximum	Minimum	Mean
Tikhonov + TSVD	6.74	1.29	3.08
TSVD	8.31	1.70	4.23
Tikhonov	7.56	1.53	3.53

Note: TSVD = truncated singular value decomposition.

). In other words, the improvements are 27.2% and 12.7% for the combined method over TSVD and Tikhonov regularization.

3.2. Total Mass Variations

As expected, the time series of mass flux solution in the TP (i.e., within the black boundary in Figure 1) from April 2002 to April 2019 exhibited strong annual changes with a decreasing secular trend at a rate of −5.6 ± 4.2 Gt/year, −7.9 ± 3.9 Gt/year, −6.8 ± 5.2 Gt/year, and −8.6 ± 5.8 Gt/year with combined TSVD and Tikhonov regularization, TSVD only, Tikhonov regularization only, and P4M6 + 400 km Gaussian filtering [10,11,12] applying to mass points exactly those used in the combined method, respectively (Figure 4). Note that the TSVD solution suffers unrealistic south-north strips (see Figure 5d), indicating that using TSVD alone is not suitable to recover mass flux solution in the TP.

To investigate the secular trend and the seasonal mass variations, we established a fitting function to estimate the linear rate, (semi-) annual variations, as well as to remove a 161-day S₂ alias that would remove an error due to incorrect modeling of the S2-tide, as follows,

$${\mathrm{EWH}\left(t\right)=a·t+b+A}_{1}cos\left(\frac{2\mathsf{\pi }}{T_1}{t-\theta }_1\right)+A_{2}cos\left(\frac{2\mathsf{\pi }}{T_2}{t-\mathsf{\theta }}_2\right){+A}_3\cos \left(\frac{2\mathsf{\pi }}{T_3}{t-\theta }_3\right)+\mathsf{\Delta },{\mathrm{T}}_1=1,{\mathrm{T}}_2=0.5,{\mathrm{T}}_3=\frac{161}{365.25},$$

(14)

where, A₁, A₂, A₃, ${\theta }_1$, ${\theta }_2$, ${\theta }_3$ stand for the annual, semi-annual, and 161-day amplitudes and phases; t is the time tag in years; b is bias parameter; $\mathsf{\Delta }$ is the residuals. We realize that Equation (14) represents a mathematical model that enables an unbiased trend estimate, rather than a tool for assessing solution errors in the presence of real geophysical signals which do not necessarily follow this model. However, we argue that (1) it is indeed our primary aim to derive the mass trends, and (2) several modeling studies have confirmed that a trend plus (semi-) annual model describes snow accumulation, groundwater, and surface water change in the TP region quite well [37,38].

The similarity of the annual phase retrieved by all methods in

Table 2. Linear trend and annual signals of the total mass variation in the TP.

Method	Trend (Gt/year)	Annual		RMSE (cm)	RMS Ratio
		Amplitude(cm)	Phase (°)
Tikhonov + TSVD	−5.6 ± 4.2	2.8 ± 0.5	226.8 ± 14.4	1.9	1.21
TSVD	−8.9 ± 5.9	2.2 ± 1.9	262.6 ± 34.3	1.6	0.45
Tikhonov	−6.8 ± 5.2	2.3 ± 0.5	220.2 ± 26.4	2.2	0.78
P4M6 + 400 km	−8.6 ± 5.8	2.3 ± 0.6	223.1 ± 23.5	1.8	1.16

Note: significant at 2-σ level for trends, annual amplitudes, and phases. RMSE = root mean squared error; RMS = root mean squared.

suggests that the mass storage of the Tibetan Plateau reaches a maximum in September, likely related to the large amount of precipitation in the TP between June and September which is brought by the southwest monsoon of the Indian Ocean [39,40]. Assume for the moment the fitting curve represents the true signal, we further derive an estimate for the root mean squared error (RMSE), the standard deviation of the residuals, and root mean squared (RMS) ratio defined as,

$$\mathrm{RMS}\mathrm{Ratio}=\frac{\mathrm{RMS}\left(\mathrm{fitting}\mathrm{curve}\right)}{\mathrm{RMS}\left(\mathrm{residuals}\right)},$$

(15)

where RMS(fitting curve) refers to the latitude weighted RMS of the fitting curve; RMS(residuals) stands for RMS of residuals $\mathsf{\Delta }$. The RMSE of the combined method is 1.9 cm, close to the RMSE of filtering, while that of the Tikhonov regularization is 2.2 cm (Table 2). The RMS ratio of 1.21 by the combined method is higher than those of 0.45, 0.78, and 1.16 by TSVD, Tikhonov regularization, and filtering, respectively, indicating that the combined method can retrieve annual and semi-annual signals at a higher confidence level. The improvements of the combined method relative to the Tikhonov regularization is 14% for RMSE, and 55% for RMS ratio.

When evaluated over the same period, the secular trend of total mass variations in the TP that we find with the combined method at 1° resolution is among the range of previous studies with different datasets (

Table 3. Linear trend in the TP of different datasets and previous studies.

Method	Time Intervals	GRACE Data	Trend (Gt/year)	Trend of Combined Method (Gt/year)
TSVD + Tikhonov	Apr 2002–Apr 2019	GFZ Release 06	−5.9 ± 4.3	−5.6 ± 4.2
Jacob et al. [8]	Jan 2003–Dec 2010	CSR Release 04	−4 ± 20	−2.3 ± 5.7
Yi and Sun [35]	Jan 2003–Dec 2012	CSR Release 05	−7.8 ± 5.7	−4.6 ± 5.3
Zou et al. [37]	Aug 2002–Dec 2016	CSR Release 05	−6.2 ± 1.7	−7.6 ± 3.8

). The use of GeoForschungsZentrum (GFZ) release 06 data agrees with the results of the Center for Space Research (CSR) release 06 within 5%.

3.3. Mass Variations Distribution

Figure 5 presents the spatial distribution of the secular trend, derived from the mass variations in TP: with the combined method, with the Tikhonov regularization, with P4M6 + 400 km Gaussian filter, with TSVD, with the filter minus those of the combined method, and with the Tikhonov regularization minus those of the combined method, respectively. The filtering results show signals relatively similar with the combined method; however, the TSVD-only method still exhibits significant south-north stripes. The trend map in Figure 5a–c shows the mass accumulation in the Inner TP at a maximum trend of 2.4 cm/year with the combined method. Two other significant signal sources are the Indus Basin located in North India and the Brahmaputra Basin located in the central and eastern Himalayas. The minimum secular trend from the mass flux solution derived from the combined method is found at −11.2 cm/year in the Indus Basin and −5.4 cm/year in the Brahmaputra Basin, while with Tikhonov regularization these trends remain as −8.9 cm/year and −4.2 cm/year, respectively. We speculate that in addition to the dampening by Tikhonov regularization, the relatively weaker signals could also be blamed by the low spatial resolution because a larger mass point covers more area and the mass variations in a large area will be smoothed if the strong signal is averaged together with weak ones in the same mass point. Note that the trend map is arranged from −2 to 2 cm in Figure 5e,f. Figure 5e shows the trends of the filtering solution minus those of the combined method, indicating that the result of the combined method presents no stripe pattern, just like the two-step filtering technique which is designed to remove the south-north stripes. Figure 5f gives the trends of the Tikhonov only method minus those of the combined method, showing south-north stripes. Since Figure 5e shows no evident south-north stripes, it is the result of Tikhonov regularization that contains stripe patterns.

The geophysical processes causing mass variations are complicated in the TP. The strong mass loss signal in north India is likely caused by the anthropogenic groundwater depletion, confirmed by groundwater-level monitoring data [41]. Global Positioning System observes long-term uplift rates as 0.5–0.7 cm/year in the TP, which might be isostatically compensated by an increasing mass deficiency at depth [42]. Gardner et al. [43] suggested rapid thinning rates of mountain glaciers in different sub-regions based on glacier elevation variations. The recharge pathway of the melt water remains unclear, and tends to sink into the ground, whose local storage capacity is limited due to permafrost [8].

We then visualized the RMSE of each mass point in the TP to compare with those from the Tikhonov regularization, and from the filtering results (Figure 6). Note that the RMSE map is arranged from 0–30 cm and 0–50 cm in Figure 6c,d, respectively. The latitude weighed RMSE of the TP with the combined method was the smallest among the three methods (9.6 cm vs. 10.4 cm for Tikhonov, 12.2 cm for filtering, and 23.0 cm for TSVD), which is consistent with the improved RMSE and RMS ratio derived from the total mass variations. The largest difference of the RMSE is located in the southeast of the Indus Basin, reaching 22.9 cm from the combined method, 25.3 cm from the Tikhonov regularization, and 29.7 cm when using the filtering method. Along the Karakoram Mountains and the Himalayas Mountains, the larger RMSE indicates that the Tikhonov regularization and the filtering methods perform worse in retrieving secular and seasonal signals than the combined method.

4. Summary

We introduced a combined TSVD and Tikhonov regularization method to retrieve mass flux solutions at 1° resolution in the TP. This combined method truncates the terms beyond the native resolution of GRACE/GRACE-FO data and dampens the errors in higher degree and order components by Tikhonov regularization. Of course, the number of degrees of freedom in the truncated normal equation is approximately equal to those directly parameterized as 2°. The improvements for the combined method over TSVD and Tikhonov regularization are 27.2% and 12.7% in terms of MSE roots from 2002 to 2019. Furthermore, using TSVD or Tikhonov regularization alone leaves corresponding mass flux solutions in the TP with unrealistic south-north stripes, while the result of the combined method shows no evidence of stripes.