## 2 Materials and Methods
### Data Processing
#### 2.1.2 GRACE Data and Processing  
In this paper, we use the Green's function convolution method to calculate the influence of regional grid point loading on the stations. Assuming that the mass load on the Earth's surface expressed in terms of the equivalent water height \(h_{w}(\theta,\lambda)\) is known, then the surface fluid mass load can be expanded into the summation form of the spherical associated Legendre function [23]:  
\[h_{w}=\sum_{n=0}^{N_{\max}}\sum_{m=0}^{n}\left(\Delta C_{nm}^{h}\cos m\lambda +\Delta S_{nm}^{h}\sin m\lambda\right)\ \overline{P}_{nm}(\cos\theta) \tag{1}\]  
where \(\theta\), \(\lambda\) denote the residual latitude and longitude east of the load point; \(\Delta C_{nm}^{h}\), \(\Delta S_{nm}^{h}\) are the spherical harmonic coefficients of the global surface load equivalent water height expansion; \(\overline{P}_{nm}\) is the \(4\pi\) normalized associated Legendre function; \(n\) and \(m\) denote the order and number, respectively; \(N_{\max}\) is the maximum order of truncation, which is related to the spatial resolution.  
Assume that \(\rho_{e}\) is the average density of the earth, \(\rho_{w}\) is the density of water, \(R\) is the average radius of the earth, \((\theta,\ \lambda)\) and \((\theta^{\prime},\ \lambda^{\prime})\) denote the residual latitude and longitude east of the point to be sought and the load point, respectively, and \(P_{n}\) is the Legendre function; \(l^{\prime}{}_{n}\) and \(l^{\prime}{}_{n}\) are the load Left number. According to the elastic load theory, the equation of the displacement of the surface in 3 directions caused by the mass load in the local Cartesian coordinates coordinate system:  
\[U=\sum_{n=0}^{N_{\max}}\frac{3\rho_{w}}{\rho_{e}}\times\frac{l^{\prime}{}_{n }}{2n+1}\sum_{m=0}^{n}\left(\Delta C_{nm}^{h}\cos m\lambda+\Delta S_{nm}^{h} \sin m\lambda\right)\ \overline{P}_{nm}(\cos\theta) \tag{2}\]  
\[N=\sum_{n=0}^{N_{\max}}-\frac{3\rho_{w}}{\rho_{e}}\times\frac{l^{\prime}{}_{n }}{2n+1}\sum_{m=0}^{n}\left(\Delta C_{nm}^{h}\cos m\lambda+\Delta S_{nm}^{h} \sin m\lambda\right)\ \frac{\partial}{\partial\theta}\overline{P}_{nm}(\cos\theta) \tag{3}\]  
\[E=\sum_{n=0}^{N_{\max}}\frac{3\rho_{w}}{\rho_{e}}\times\frac{l^{\prime}{}_{n }}{2n+1}\sum_{m=0}^{n}m\left(\Delta S_{nm}^{h}\cos m\lambda-\Delta C_{nm}^{h} \sin m\lambda\right)\ \frac{\partial}{\partial\theta}\overline{P}_{nm}(\cos\theta) \tag{4}\]  
We used the load Love number provided by [PERSON] et al. [24] and the RL05 equivalent water height data product with a resolution of \(0.5^{\circ}\times 0.5^{\circ}\times 10\) days provided by the Centre National d'Etudes Spatiales/Groupe de Recherche de Geodesie Spatiale (CNES/GRGS) [25] to calculate the crustal displacement of gravity load that could facilitate comparisons with GNSS data. Atmospheric and non-tidal ocean loading models were corrected whileprocessing these data. To extract and analyze data periodic terms using the WT method and ensure the continuity of the signal frequency domain, we used the cubic-polynomial interpolation method to provide the missing data. Seasonal analysis was not performed for long periods of missing data between GRACE and GRACE-FO during 2017-2018.