## 2 Materials and Methods
### Geoid Height Computation from GGMs  
For global gravity field analysis, a spherical harmonic series can represent the geoid heights (\(N\)). The longer wavelength geoid heights (\(N\)) contributed by the geopotential models were computed at position (\(\varphi_{P}\), \(\lambda_{P}\)) using Equation (1) [20]:  
\[N(\varphi_{P},\ \lambda_{P})=R\sum_{l=2}^{L}\sum_{m=0}^{l}\left[\overline{C}_{ lm}cosm\lambda_{P}+\overline{S}_{lm}sim\lambda_{P}\right]\overline{P}_{lm}( sin\varphi_{P}) \tag{1}\]  
where \(R\) is the Earth's mean radius; \(\overline{C}_{lm}\) and \(\overline{S}_{lm}\) are the fully normalized spherical harmonic coefficients; \(\overline{P}_{lm}(\sin\varphi_{P})\) is the fully normalized associated Legendre function; and \(l\) and \(m\) are the degree and order, respectively. The finite series is usually truncated at the maximum degree of expansion \(l=L\). The series coefficients allow the determination of geoid heights using Equation (1). Table 2 summarizes the main characteristics of the three high-degree GGMs that were utilized in this study. In Table 2, \(S\) stands for satellite (e.g., GRACE, GOCE, and LAGEOS), \(A\) for altimetry, and \(G\) for ground data (e.g., terrestrial, shipborne, and airborne measurements). The GGMs that were used in this study, as shown in Table 2, were provided by ICGEM in the zero-tide system with respect to a geometrically fixed GRS80 reference ellipsoid, and were calculated from GRACE and GOCE gravitational satellite observations, satellite-altimetry-derived gravity data, land gravity data, etc. For the geoid height accuracy evaluation in this study, high-degree GGMs developed from the combined contributions of various types of satellite tracking data, gravity data, and satellite altimetry data were considered, as shown in Table 2. The geoid heights (\(N_{GGM}\)) obtained through spherical harmonic synthesis from high-degree GGMs were interpolated into the 1182 GNSS/leveling UCPs to assess the accuracy of geoid heights derived from the three different GGMs with maximum degree 2190. The GNSS/leveling geoid heights (\(N_{GNSS/leveling}\)) at the 1182 UCPs were computed by subtracting the known orthometric height (\(H_{leveling}\)) calculated from the precise leveling between the BM and the \begin{table}
\begin{tabular}{c c c c} \hline \hline
**Geopotential Model** & **Year** & **Max. Degree** & **Data Type** \\ \hline GECO & 2016 & 2190 & S (GOCE), EGM2008 \\ EIGEN-6C4 & 2014 & 2190 & S (GOCE, GRACE, LAGEOS), A, G \\ EGM2008 & 2008 & 2190 & S (GRACE), A, G \\ \hline \hline \end{tabular}
\end{table} 
Table 2: Characteristics of the GGMs used in this study. UCP from the known ellipsoidal height (\(h_{GNSS}\)) obtained through GNSS observations, as given in Equation (2). \[N_{GNSS/\text{\emph{leveling}}}=h_{GNSS}-H_{\text{\emph{leveling}}}. \tag{2}\]  
The GNSS-derived geoid heights are invariably different from the values interpolated from a gravimetrically derived geoid model, and are influenced by the datum inconsistencies, biases, and errors associated with the independently derived ellipsoidal and orthometric height data [21; 22; 23; 24]. An external evaluation of the quality of a gravimetric geoid model or GGM can be performed by comparing its interpolated values \(\left(N_{i}^{Model}\right)\) at the 1182 UCPs with the corresponding GNSS/leveling-derived geoid heights \(\left(N_{i}^{GNSS/\text{\emph{leveling}}}\right)\).