Assessing Height Variations in Qinghai-Tibet Plateau from Time-Varying Gravity Data and Hydrological Model

: Height variations caused by mass change make an important contribution to the tectonic uplift of the Qinghai-Tibet Plateau (QTP). To study the deformation attributable to hydrological loading and real potential tectonic vertical motion, satellite gravity data from the Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On (GRACE-FO) with data from the Global Land Data Assimilation System (GLDAS) and Global Positioning System (GPS) are adopted to estimate height variations in QTP. Based on spherical harmonic function (SHF) and Green’s function (GF), the results show the trend of height variations is unevenly distributed in the spatial domain. The SHF indicated that the rate in the southwest of the QTP is ~1 mm/year, while the northern and eastern show a subtle decreasing trend, which indicates hydrological loading is not the main cause of the uplift observed with GRACE. The maximum annual amplitude of height variations is ~12 mm, reaching the annual maximum around February to March. The average correlation coefﬁcients of SHF, and GF height variations with GPS heights are 0.70 and 0.82, respectively. Based on cross wavelet transform, it is concluded that there are annual signals between the height variations derived from GPS with GRACE (-FO) and GLDAS. Finally, the tectonic vertical motion in the QTP is given by removing the effect of hydrological loading, which shows most GPS stations are uplifted at a rate of 0.06 mm/year–1.97 mm/year.


Introduction
The Qinghai-Tibet Plateau (QTP, 25 • N-40 • N, 74 • E-104 • E) with an average altitude of more than 4000 m, known as the third pole, is a tectonic extrusion and deformation zone as result of the collision between the Indian and Eurasian plates [1]. Accurate determination of height variations over the QTP is important for the study of plate movement and evolution, and ecological changes [2]. The continuous uplift and geophysical processes of the QTP have always been a controversial issue [3]. The traditional method of obtaining height variations mainly relies on precise levelling measurements [4]. With the development of space geodesy technology, GPS network provides high-precision quantitative data of crustal motion by observing the coordinate changes of GPS stations [5][6][7]. The GPS vertical coordinates mainly include linear trends and seasonal variations at annual scales, while the causes of linear trends include tectonic and non-tectonic motions (e.g., mass change) and other factors [6]. Liang et al. [8] estimated the three-dimensional velocity field for the QTP, with the Himalayan Mountains uplifting of 2 mm/year. Moreover, repeated absolute gravity measurements are also sensitive to mass migration and vertical deformation of the surface. Xing et al. [9] combined high-precision absolute gravity measurements and GPS data to quantify the average uplift rate of 1.4 mm/year on the QTP. Van Camp et al. [10] discussed the tectonic interpretation from stations in Tibet using absolute gravity survey.
interpolation is performed to fill in missing data for one or two months during the separate GRACE or GRACE-FO mission. However, it is worth noting that linear interpolation may also introduce the uncertainty, especially when flood or drought events occur [33].
Since the GRACE gravity model is centered on the center of mass, the degree-1 coefficients (C 10 , C 11 , S 11 ) estimated by Satellite Laser Ranging (SLR) are used, thus taking into account the geocenter motion [34]. Meanwhile, C 20 derived from GRACE (-FO) show substantially larger uncertainty, in addition to the failure of an accelerometer on the GRACE satellite since October 2016, which caused a lower accuracy of the C 30 term [35]. The SLRderived C 20 and C 30 are suitable for replacing any problematic GRACE (-FO) estimates. The effect of glacial isostatic adjustment (GIA) effect is corrected by ICE-6G_D model [36]. The DDK3 filtering method is taken to process the time-varying gravity spherical harmonic coefficients to attenuate the striping errors and correlation errors of higher degree terms [37].

GLDAS Noah Hydrological Model
The surface water storage changes from Global Land Data Assimilation System (GLDAS) Noah V2.1 hydrological model [38] are used to estimate the hydrological loading, with a spatial-temporal resolution of 1 month and 1 degree, respectively. The 0-2 m soil moisture, canopy water, and snow are derived from the GLDAS and summed as the total land water storage, then the monthly water storage changes are obtained by subtracting the average value of water storage during the study period.

GPS Data
Daily GPS data is provided by China Earthquake Networks Center, Nation Earthquake Data Center (https://data.earthquake.cn/, accessed on 7 November 2021). The data are continuous observations from the Crustal Movement Observation Network of China (CMONOC). Thirty-four GPS continuous stations located on the QTP are used to study the vertical displacement ( Figure 1). Among them, LHAS, DLHA, and XNIN have been observed since 1999, while the rest of the stations have been recording data since around 2010. The raw GPS observations are processed by the GAMIT/GLOBK software [39] with the orbit and clock bias provided by International GNSS Service (IGS). The tropospheric delays are corrected using the Global Mapping Function (GMF) model and the iono- The raw GPS observations are processed by the GAMIT/GLOBK software [39] with the orbit and clock bias provided by International GNSS Service (IGS). The tropospheric delays are corrected using the Global Mapping Function (GMF) model and the ionospheric effects are weakened based on the ionosphere-free combinations. The antenna phase center is corrected by the IGS antenna model, the solid tide, and the polar tide is corrected using the International Earth Rotation Service 2003 (IERS03) model, and the ocean tide loading effects are corrected by the Finite Element Solution 2004 (FES2004) model. The GNSS Missing Data Interpolation Software (GMIS) [40] is used to fill in the missing GPS data.
The signal of environmental loading in the GPS vertical coordinate mainly includes non-tidal atmospheric loading (NTAL), non-tidal ocean loading (NTOL), and hydrological loading. Since the surface vertical deformation derived from GRACE (-FO) gravity solutions mainly contains hydrological loading, NTAL and NTOL in the GPS need to be removed. The Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2) [41] atmospheric loading model, and the Max Planck Institute Ocean Model (MPIOM06) [42] non-tidal ocean loading model provided by the International Mass Load Service (http://massloading.net/, accessed on 7 November 2021) [43] are taken to process GPS vertical coordinate in this paper. The annual amplitudes of NTAL and NTOL at the GPS stations are shown in Figure 2. Compared with the NTOL, NTAL has a greater impact on the QTP, with the maximum annual amplitude up to 4.1 mm. The amplitude of NTAL is inversely proportional to the altitude, and the edge of the study area is more affected by NTAL than the internal high-altitude region. The maximum annual amplitude of the NTOL is about 0.4 mm, with an east-west distribution, i.e., high in the west and low in the east. In addition, the contribution of the thermoelastic effect to the GPS bedrock is at the millimeter level on the global scale [44]. The bedrock thermal expansion [45] of 34 GPS stations is estimated using 1 km monthly temperature data [46] in this paper. Figure 3 shows the annual amplitude of bedrock thermal expansion, which increases from 0.76 mm to 1.49 mm with increasing latitude. The bedrock thermal expansion is deducted from the GPS vertical displacement to eliminate the effect of temperature change. In addition, the contribution of the thermoelastic effect to the GPS bedrock is at the millimeter level on the global scale [44]. The bedrock thermal expansion [45] of 34 GPS stations is estimated using 1 km monthly temperature data [46] in this paper. Figure 3 shows the annual amplitude of bedrock thermal expansion, which increases from 0.76 mm to 1.49 mm with increasing latitude. The bedrock thermal expansion is deducted from the GPS vertical displacement to eliminate the effect of temperature change.
millimeter level on the global scale [44]. The bedrock thermal expansion stations is estimated using 1 km monthly temperature data [46] in this p shows the annual amplitude of bedrock thermal expansion, which increases to 1.49 mm with increasing latitude. The bedrock thermal expansion is ded GPS vertical displacement to eliminate the effect of temperature change.

Height Variations
Height variations ∆H represent the change in orthometric height (ab in this paper, which are determined as follows:

Height Variations
Height variations ∆H represent the change in orthometric height (above the geoid) in this paper, which are determined as follows: where ∆N and ∆h represent the temporal variation of geoid height and surface vertical displacement, respectively. The gravity field is usually represented in the form of the geoid. Temporal variation of geoid height ∆N can be estimated by the spherical harmonic function [47]: (∆C lm cos mλ + ∆S lm sin mλ)P lm (sin ϕ) (2) where (r, ϕ, λ) are geocentric radius, latitude, and longitude on the ellipsoid, respectively, GM is the constant of earth gravitation, a is the semi-major axis of the reference ellipsoid, l and m are the degree and order of the spherical harmonic coefficients, respectively, l max (= 60) is the maximum degree of the spherical harmonic coefficient, P lm are the fully normalized Legendre functions of degree l and order m, γ represents the normal gravity at the computation point, ∆C lm and ∆S lm derived from GRACE (-FO) is the residual of fully normalized stokes coefficients, from which the average stokes coefficients over the study period have been subtracted. The surface vertical displacement ∆h is determined by two methods in this paper, the spherical harmonic function (SHF), and Green's function (GF), which are described in Sections 3.1.1 and 3.1.2, respectively.

Spherical Harmonic Function
The height variations ∆H SHF are estimated by the GRACE (-FO) time-varying gravity field model. Based on the earth loading deformation and gravity field theory [48], the surface vertical displacement is determined by the spherical harmonic coefficient derived from GRACE (-FO) as follows: where ∆h SHF is the geometric height change (above the ellipsoid) assuming that the gravity changes measured by GRACE are the result of changes in mass loading, ρ w is the water density, ρ ave is the earth's average density, h l is the load Love numbers of degree l, and surface density Stokes coefficient ∆C σ lm and ∆S σ lm defined as: Combining Equations (3) and (4) obtains the vertical deformation due to mass change by the time-varying GRACE (-FO) gravity model: where the load Love numbers h l and k l based on the Preliminary Reference Earth Model (PREM) [49] obtained from Wang et al. [50].

Green's Function
The height variations ∆H GF are determined by Green's function (GF) according to the elastic half-space model [24]. The ∆h GF can be determined as follows [51], which is a simplification of the usual load theory based on Farrell [24]: where ∆EWT is the time-variation of equivalent water thickness derived from the GLDAS hydrologic model, R 0 is estimated by assuming a disk loading corresponding to the radius based on the spatial resolution of GRACE satellite missions' data (R 0 = ca. 167 km), g is the gravitational acceleration, v is the Poisson's ratio, and E is Young's modulus. The Poisson's ratio v is generally set to 0.25 [52]. The type of rock in the QTP is mainly syenite-diopside and granite-biotite; therefore, the E value is estimated to be 50 GPa in this study.

Least Square Fitting
The least squares spectral analysis is used to estimate the long-term trend, annual and semi-annual signals of height variations, which is as follows [53]: where t i is the epoch in unit of year, a and b are constant term and linear long-term trend, respectively, the coefficients (c, d) and (e, f ) describe the annual and semi-annual amplitudes, respectively, ν i is the residual. Time series of annual and semi-annual amplitudes (A ann and A semi−ann ) and phase (ϕ ann and ϕ semi−ann ) can be expressed as:

Height Variations on the QTP
According to the GRACE (-FO) temporal gravity field model, the geoid changes ∆N are obtained by Equation (2), and the vertical surface deformation ∆h estimated by SHF, and GF methods are derived by Equations (5) and (6) ∆H is determined by the ∆N and the ∆h through Equation (1). The time series of ∆N, ∆h, and ∆H are shown in Figure 4. The trend and annual signals based on the least squares fitting are shown in Table 1.   From Figure 4, it can be seen that ∆N, ∆h, and ∆H in the QTP all show strong seasonal fluctuations, most notably annual variations. Moreover, the long-term trend is not obvious, only showing a subtle upward trend, which suggests the hydrological mass of the QTP has maintained a long-term balance during the study period. It is confirmed the hydrological loading is not the main reason for the continuous uplift observed with GRACE and GLDAS. The hydrological loading is not significant in the QTP, with the annual amplitude of about 2.12 ± 0.12 and 1.73 ± 0.04 mm by SHF and GF, respectively ( Figure 4b). The annual phase difference between ∆N and ∆h is about 6 months, revealing the relationship between hydrological mass change, geoid height, and surface vertical displacement. The average ∆H of the QTP varies significantly in different periods, for example, the difference of ∆H SHF between September 2003 and April 2017 reaches 18 mm (Figure 4c). Additionally, ∆H reached the maximum and minimum values around February to March and August to September, respectively, and the maximum difference between the SHF, and GF is about 4 mm in the same period.
The long-term trends and annual amplitudes of ∆H in the QTP based on the least squares fitting are shown in Figures 5 and 6, respectively. The trends of height variations in the study area are characterized by uneven spatial distribution ( Figure 5), roughly bounded by block boundaries, with large differences on both sides. The trend of uplift in the southwestern part of the QTP near the Indian plate is about 1 mm/year, while the northern and eastern shows a subtle decrease (Figure 5a), which agree with the results in Duan et al. [25]. Overall, ∆H GF derived by the GLDAS hydrological model are smaller than the ∆H SHF to some degree.

Comparison of Height Variations at GPS Stations
The height derived directly from GPS is the surface vertical displacement, noted as   [52], which are significantly different from the results of time-varying gravity and hydrological models.

Comparison of Height Variations at GPS Stations
The height derived directly from GPS is the surface vertical displacement, noted as   [52], which are significantly different from the results of time-varying gravity and hydrological models.   Figure 6 suggests there is a good agreement between the annual amplitudes of ∆H estimated by SHF, and GF, which all show a latitudinal distribution, i.e., the annual amplitudes decrease gradually with increasing latitude. The maximum annual amplitudes of ∆H SHF , and ∆H GF are 13.06 mm, and 12.15 mm, respectively, where the maximum values are located in the southern part of the QTP, while the rest of the region is not significant less than 5 mm.

Comparison of Height Variations at GPS Stations
The height derived directly from GPS is the surface vertical displacement, noted as ∆h GPS in this study. It is necessary to subtract the geoid change ∆N from averaging ∆h GPS by each month, which converts ∆h GPS into height variations ∆H GPS based on Equation (1). Table 2 shows the long-term trends of height variations ∆H at GPS stations based on the least squares fitting. The long-term trends derived from GRACE (-FO) and GLDAS are fairly close to each other at some GPS stations (e.g., DLHA, QHLH, QHMY, and SCJL), but Remote Sens. 2022, 14, 4707 9 of 20 differ considerably in most stations. GPS-derived ∆H include not only the hydrological loading but also the effect of tectonic movements [52], which are significantly different from the results of time-varying gravity and hydrological models. The annual amplitude of ∆H at the GPS stations is given in Figure 7, where the length of the arrow represents the magnitude of the annual amplitude and the angle turned clockwise from north represents the annual phase. The annual phase at 34 GPS stations is roughly located in the range of 10 • to 50 • (Figure 7), corresponding to the annual maximum around January to March, and the annual phase of ∆H GPS is earlier than other results. There is a phase difference of about 25 days between the ∆H SHF derived from GRACE (-FO) time-varying gravity model with the ∆H GF estimated by GLDAS hydrological model. One possible cause of the obvious phase difference is the local groundwater-induced porous response, which is not considered in our elastic loading model and may lead to a phase shift of ground deformation [54]. The annual amplitudes of ∆H obtained by different methods are approximately equal in the southern parts of the QTP, while ∆H GPS in western and northern parts are larger than the ∆H SHF or ∆H GF . have approximately equal annual amplitudes. In addition, the different results showed significant seasonal fluctuations, with the most significant annual periodic signals. Vertical displacement uplifts in winter and descends in summer. The annual peak of subsidence occurs in summer, which coincides with the time of year when the summer monsoon brings heavy precipitation and increases the surface mass loading [16]. It can be inferred that the seasonal fluctuations in the QTP are largely influenced by the hydrological loading due to precipitation, with the surface loading increasing during the summer and decreasing during the winter. The time series of vertical displacement with the long-term trend removed at LHAS, QHBM, SCLH, and XZCY stations are shown in Figure 8, which suggests there is the largest amplitude of ∆h GPS , while ∆h SHF , and ∆h GF have approximately equal annual amplitudes. In addition, the different results showed significant seasonal fluctuations, with the most significant annual periodic signals. Vertical displacement uplifts in winter and descends in summer. The annual peak of subsidence occurs in summer, which coincides with the time of year when the summer monsoon brings heavy precipitation and increases the surface mass loading [16]. It can be inferred that the seasonal fluctuations in the QTP are largely influenced by the hydrological loading due to precipitation, with the surface loading increasing during the summer and decreasing during the winter.

Correlation Analysis of Height Variations at GPS Stations
The correlation coefficient is used to evaluate the degree of linear correlation between two time series. Figure 9 presents the correlation coefficients of ∆H GPS time series with ∆H SHF , and ∆H GF , respectively, with the long-term trend removed. ∆H GPS are positively correlated between the ∆H estimated by the GRACE and GLDAS. Compared with the ∆H SHF (Figure 9a), the correlation coefficients of ∆H GF and ∆H GPS (Figure 9b) are larger, with the average correlation coefficients of 0.70, and 0.82 between the SHF/GF and GPS, respectively, which indicates that GPS, GRACE (-FO), and the GLDAS detect some degree of the same geophysical signal while mixing the non-hydrological signal in the GPS height. The correlation strength varies among GPS stations, and the correlation is stronger in the southern part of the QTP near the Indian and Myanmar plates than in the northern part, which may be due to the more significant hydrological loading caused by the whole TWS in the southern part.

Correlation Analysis of Height Variations at GPS Stations
The correlation coefficient is used to evaluate the degree of linear correlation between two time series. Figure 9 presents the correlation coefficients of  Figure  9b) are larger, with the average correlation coefficients of 0.70, and 0.82 between the SHF/GF and GPS, respectively, which indicates that GPS, GRACE (-FO), and the GLDAS detect some degree of the same geophysical signal while mixing the non-hydrological signal in the GPS height. The correlation strength varies among GPS stations, and the correlation is stronger in the southern part of the QTP near the Indian and Myanmar plates than in the northern part, which may be due to the more significant hydrological loading caused by the whole TWS in the southern part.

Wavelet Coherence analysis of Height Variations at GPS Stations
Since correlation analysis is sensitive to phase, but cannot evaluate the consistency of amplitudes. The wavelet coherence analysis [55] could characterize the signal in the spatial and frequency domains, so the continuous wavelet transform (CWT) is used to analyze relationships in the time-frequency space of height variations derived from GPS, GRACE (-FO), and GLDAS. The method of wavelet coherence analysis can be found in

Wavelet Coherence analysis of Height Variations at GPS Stations
Since correlation analysis is sensitive to phase, but cannot evaluate the consistency of amplitudes. The wavelet coherence analysis [55] could characterize the signal in the spatial and frequency domains, so the continuous wavelet transform (CWT) is used to analyze relationships in the time-frequency space of height variations derived from GPS, GRACE (-FO), and GLDAS. The method of wavelet coherence analysis can be found in Grinsted et al. [55].
Station LHAS is used as an example for wavelet analysis in this paper. Figure    XWT reveals the phase relationship in the high energy region common to both data within the time-frequency space. From Figure 11, it can be seen that ∆H GPS exists a period of 8 to 16 months (0.7-1.3 years) with ∆H SHF (∆H GF ) and passes the significance test. Most arrows in this period face to the right, indicating the same phase. In addition, there is a semiannual periodic signal of some energy (2012 and 2014), but the energy is not significant.   WTC compensates for the deficiency of XWT in the low-energy region and estimates the correlation between the two datasets. Figure 12 shows that in addition to the annual period (0.7-1.3 years), there are also semi-annual and 32-month (2.7 year) periods between ∆H derived from GPS and GRACE/GLDAS. The 2.7-year periodic oscillation between ∆H GPS and ∆H SHF passed the significance test in the whole period, indicating that ∆H detected by GRACE (-FO) in the high-frequency part is the main factor of the GPS height variations. The 2.7-year period also exists between GPS and GF from 2002 to 2011, and the specific reasons for this period need further analysis.  Height variations derived from different datasets mainly present the annual periodic oscillation signal. To quantitatively express the relative phase relationship between two signals [56], the average values of WTC-based semblance (Figure 13a) and average relative phase angle (Figure 13b) are estimated. Height variations derived from different datasets mainly present the annual periodic oscillation signal. To quantitatively express the relative phase relationship between two signals [56], the average values of WTC-based semblance (Figure 13a) and average relative phase angle (Figure 13b) are estimated.
The average values of WTC-based semblance for most GPS sites are close to 1 (Figure 13a), which indicates the height variations derived from GRACE (-FO) or GLDAS almost explain GPS height in annual periodic signals, but not completely. In addition, the average WTC-based semblance of GPS and SHF is significantly lower than that of GPS and GF at stations DLHA, QHLH, QHMY, XINI, XZGE, and XZRT, which may be since groundwater loading is not the main factor causing height variations at these stations. The average semblances of the annual periodicity for most of the other stations are~0.98, showing there is a strong geophysical coherence between the height variations derived by GPS and GRACE/GLDAS. After removing NTAL, NTOL, and bedrock thermal expansion, the hydrological mass migration is the main reason for the seasonal height variations of the GPS.
The average relative phase angle of ∆H GPS is earlier than other results, and that estimated by GPS and GF is closer, whose difference is smaller than GPS and SHF (Figure 13b). In this regard, we find that the height variations derived from GLDAS are more consistent with GPS in the range of annual fluctuation.
is a strong geophysical coherence between the height variations derived by GPS and GRACE/GLDAS. After removing NTAL, NTOL, and bedrock thermal expansion, the hydrological mass migration is the main reason for the seasonal height variations of the GPS.
The average relative phase angle of GPS H Δ is earlier than other results, and that estimated by GPS and GF is closer, whose difference is smaller than GPS and SHF ( Figure  13b). In this regard, we find that the height variations derived from GLDAS are more consistent with GPS in the range of annual fluctuation.

Differences in Height Variations Derived from GRACE (-FO) and GLDAS
In this paper, GRACE (-FO) and GLDAS are used to study the vertical displacement and height change in the QTP. From Section 4.1, it can be seen that the SHF and GF methods could reveal the long-term trend and seasonal fluctuations of height variations in the QTP, but there are still some differences between the different results. Figure 4b,c show that the annual amplitudes of the height derived from GRACE (-FO) are larger than those from GLDAS. The reason may be that GLDAS only contributes to 0-2 m soil moisture canopy water, and snow to hydrological mass changes, ignoring the effort of deep subsoil water, groundwater, lakes, wetlands, etc. While GRACE (-FO) considers gravity as hydrological mass after removing atmospheric and non-tidal oceanic influences and integrating various factors such as groundwater and deep soil water, therefore the amplitude of height variations estimated by GRACE (-FO) is larger than that of GLDAS.
We estimate groundwater by TWS derived from GRACE minus soil moisture and snow water equivalent obtained from GLDAS based on regional water balance [13,19]. The groundwater of the QTP accounts for~28% of the TWS. To be clear, the proportion of groundwater is overestimated because we ignored surface water due to data limitations. The effect of groundwater on height variations was not modelled separately because we focused on the whole hydrological loading. In general, it can be concluded that the discrepancy in height variations between GRACE (-FO) and GLDAS is mainly due to the different technical models.
In addition, the GRACE mission had some measuring and processing errors, which may related to truncation; van Dam et al. (2007) argue that the error of the GRACE gravity model is less than the millimeter level [57]. The inversion of height variations by GRACE (-FO) has certain drawbacks, which reflect the mass change at a large scale and lacks spatial resolution for the influence of small-scale material sources. There will be a large bias in areas dominated by groundwater loading because the displacement caused by groundwater is in the opposite direction of surface loading [7]. The discrepancies between GRACE and GLDAS are mainly related to their uncertainty, including initial conditions, climate forcing, land cover input, model structure, etc. [58]. There are still some unmodeled errors, such as seasonal variations in immaterial gravity loading (temperature and magnetic fields), whose effects need to be further explored.

Quantitative Assessment of Hydrological Loading on GPS Height
To quantitatively assess the effectiveness of removing the hydrological loadings in GPS vertical displacement by GRACE (-FO) and GLDAS, the weighted root-mean-square (WRMS) and its reduction could be estimated, which is defined by Equations (3) and (4) in He et al. (2017) [59].
The WRMS reduction reflects the consistency of ∆H GPS and ∆H SHF (∆H GF ) in the annual periodic signals and the closer its value is to 1 indicates that the period term of the two fits is closer to each other. Figure 14 shows the WRMS reduction ratio by removing the GRACE (-FO) and GLDAS-derived height variations from GPS-derived seasonal signals at the QTP 34 GPS stations. The average WRMS reduction between GPS and SHF/GF is 29% (Figure 14a), and 37% (Figure 14b), respectively, which indicates the effect of GPS removed by GF methods is better than that of SHF. The height variations ∆H GF estimated by GLDAS is more consistent with the annual signals of ∆H GPS , therefore the improvement of hydrological loading by GLDAS in GPS height is better.  Comparing Figure 9 with Figure 14, it can be seen that there is a certain spatial correlation between the correlation coefficient and the reduction WRMS , which are significantly larger in the south of the QTP than in other regions. The hydrological loading is the dominant factor influenced by the summer monsoon in the southern, where the GRACE (-FO) or GLDAS seasonal signals with GPS show better consistency. While the north of the study area, which is closer to the interior, is less influenced by the summer monsoon, where the correlation and consistency are significantly lower. Comparing Figure 9 with Figure 14, it can be seen that there is a certain spatial correlation between the correlation coefficient and the WRMS reduction , which are significantly larger in the south of the QTP than in other regions. The hydrological loading is the dominant factor influenced by the summer monsoon in the southern, where the GRACE (-FO) or GLDAS seasonal signals with GPS show better consistency. While the north of the study area, which is closer to the interior, is less influenced by the summer monsoon, where the correlation and consistency are significantly lower.

Tectonic Movement Rates
The height variations derived from GPS include not only surface mass loading but also the effects of internal tectonic movements [5]. The surface loading in GPS heights needs to be deducted to reveal the tectonic movements of the QTP. There are some discussions on comparing or combining vertical displacement from GRACE and GPS to study the water cycle and tectonic movements [26,27,30,31,51,55,56]. However, the use of GRACE to remove the hydrological loading in GPS may be imperfect due to differences in observational patterns. Zhang et al., argued that the loading differences between GRACE and GPS are significant both on a regional average and at a single point, and directly using GRACE to remove the seasonal signal from GPS to study the remaining tectonic signal is challenging [60]. In contrast, the results in Figure 14 from the GRACE (-FO) and GLDAS hydrological models to remove the hydrological loading signal from GPS, respectively, show a higher WRMS reduction by GLDAS.
As shown in Section 5.2, WRMS of the height variations of 34 GPS stations are reduced to different degrees by removing the hydrological loading in GPS height with GRACE (-FO) and GLDAS, and WRMS by the GF is reduced to the greatest extent. Therefore, the hydrological loading in the GPS height variations is removed by the ∆h GF and the remaining signal is considered to be caused by the tectonic motion. Figure 15 shows the vertical rates of tectonic motion of GPS station with the loading effects removed (GPS height subtracts GLDAS-derived hydrological loading based on the GF). Except for the subsidence of XZNQ, XINI, XJYT, YNZD, and SCML, all other stations are uplifting at a rate of 0.06mm/year−1.97 mm/year in the study area ( Figure 15). The uplift rate of the main Main Himalaya Thrust and the Karakorum-Jiali Fault Zone is about 1 mm/year, and individual GPS stations (e.g., XZNQ) show a subsidence component, which may be related to the rebound fault activity or local collapse of the margin. The uplift rate of the northeastern part is not significant, while the stations near the Jinsha River-Honghe Fault Zone uplift at a rate of 0.2-1 mm/year, except for the SCML station, which indicates a significant subsidence rate.
Remote Sens. 2022, 14, x FOR PEER REVIEW 18 of 2 multiple distinct motions, e.g., Jinsha River-Honghe Fault Zone. The northern part of th Jinsha River-Honghe Fault Zone, represented by SCLT and SCXC, shows a unified upli trend, which is caused by active subduction of the Yangtze Craton [26]. While the south ern part, represented by SCML, shows a subsidence trend attributed to the asthenosphere It is the same result as that of Hao et al., [61], who argued that the assumption of uniform extension does not explain the subsidence of southern Sichuan-Yunnan fragment.

Conclusions
The crust of Tibetan Plateau is still uplifting and thickening nowadays which wa revealed by the geological and tectonic results. Based on the GRACE (-FO) time-varyin gravity model, GLDAS, and GPS vertical coordinate time series, the height variation (above the geoid) of the QTP from 2002 to 2020 are obtained in this study. The overa The deformation detected by GPS is mainly affected by the local loading around 100-200 km. Some faults do not show a unified tectonic movement, but are composed of multiple distinct motions, e.g., Jinsha River-Honghe Fault Zone. The northern part of the Jinsha River-Honghe Fault Zone, represented by SCLT and SCXC, shows a unified uplift trend, which is caused by active subduction of the Yangtze Craton [26]. While the southern part, represented by SCML, shows a subsidence trend attributed to the asthenosphere. It is the same result as that of Hao et al. [61], who argued that the assumption of uniform extension does not explain the subsidence of southern Sichuan-Yunnan fragment.

Conclusions
The crust of Tibetan Plateau is still uplifting and thickening nowadays which was revealed by the geological and tectonic results. Based on the GRACE (-FO) time-varying gravity model, GLDAS, and GPS vertical coordinate time series, the height variations (above the geoid) of the QTP from 2002 to 2020 are obtained in this study. The overall height variations caused by hydrological mass on QTP during the study period show a slight increase, which indicates hydrological loading is not the main cause of the uplifting observed with GRACE.
The seasonal term (e.g., annual and semi-annual period) caused by hydrological loading is included in the GPS vertical displacement, which needs to be eliminated to reveal the tectonic motion. In terms of hydrological load on GPS height, GLDAS-derived is more consistent with GPS height than that from GRACE (-FO) based on WRMS reduction . The vertical tectonic motion of GPS stations is estimated with the GLDAS-derived hydrological loading removed. Most stations are uplifted at the rate of 0.06mm/year-1.97 mm/year. To explain the uplift rate of QTP in more detail, more intensive GPS continuous stations need to be analyzed, and more accurate physical models of surface vertical displacement such as hydrological and atmospheric loads and bedrock thermal expansion also need to be considered. Overall, the approach of estimating height variations and tectonic motion in this paper provides a valuable reference for the study of other earthquake-prone areas such as Japan, Chile, etc., and engineering works, e.g., deformation monitoring of constructions.