Geopotential Determination Based on Precise Point Positioning Time Comparison: A Case Study Using Simulated Observation

Based on the general relativity theory, the geopotential difference can be determined by comparing the change in time difference between precise clocks using the precise point positioning (PPP) time transfer technique, referred to as the relativistic PPP time comparison approach. We focused on high-precision time comparison between two high-performance clocks for determining the geopotential difference using this approach, and conducted the simulation experiments to validate this approach. In the experiment, we consider three cases for evaluating the performance of this approach with different types of atomic clocks, namely, the fractional frequency stabilities of the clocks equipped at three selected ground stations (BRUX, OPMT, and PTBB) are <inline-formula> <tex-math notation="LaTeX">$4.0 \times 10^{-13} / \sqrt {\tau }$ </tex-math></inline-formula> (Case 1), <inline-formula> <tex-math notation="LaTeX">$2.3 \times 10^{-14} / \sqrt {\tau }$ </tex-math></inline-formula> (Case 2), and <inline-formula> <tex-math notation="LaTeX">$2.8 \times 10^{-15} / \sqrt {\tau }$ </tex-math></inline-formula> (Case 3) at averaging time <inline-formula> <tex-math notation="LaTeX">${\tau }$ </tex-math></inline-formula>, respectively, and the accuracy of these clock have been evaluated to be <inline-formula> <tex-math notation="LaTeX">$5.3 \times 10^{-16}$ </tex-math></inline-formula>, <inline-formula> <tex-math notation="LaTeX">$7.8 \times 10^{-17}$ </tex-math></inline-formula>, and <inline-formula> <tex-math notation="LaTeX">$8.6 \times 10^{-18}$ </tex-math></inline-formula>. Two main conclusions can be drawn from the experimental results. First, high-performance clocks can significantly improve the precision for GNSS PPP time transfer. Compared to Case 1, the long-term stabilities of the time link OPMT-BRUX as well as PTBB-BRUX are improved in Cases 2 and 3. Second, the geopotential difference between any two stations can be determined at the decimeter level, and the accuracy of geopotential difference is consistent with the stabilities of the time links in Cases 1–3. In Case 3, the determined geopotential differences between OPMT and BRUX deviate from the EIGEN-6C4 model values by −0.64, m<sup>2</sup>/<sup>2</sup> with an uncertainty of 1.11 m<sup>2</sup>/s<sup>2</sup>, whereas the deviation error between PTBB and BRUX is 0.76 m<sup>2</sup>/s<sup>2</sup> with an uncertainty of 1.79 m<sup>2</sup>/s<sup>2</sup>. The results of this study suggest that a one decimeter-level geopotential difference between two arbitrary stations can be determined based on this approach.


I. INTRODUCTION
The geopotential plays a very important role in geodesy and has broad applications in various fields. The classic method of determining geopotential difference is based on leveling with additional gravimetry [1], which has drawbacks in that it is The associate editor coordinating the review of this manuscript and approving it for publication was Halil Ersin Soken .
labor-intensive and geographically limited [2]. To overcome these drawbacks, Bjerhammar [3] proposed that geopotential difference could be determined using a clock transportation comparison (CTC) approach [4] based on the general relativity theory (GRT) [5]. The CTC technique is based on continuously comparing the change in time difference between a fixed clock with a transportable clock [6]. The key problem lies in the difficulties associated with transporting VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ high-precision clocks [7] because it is very difficult to control ideal operating conditions (temperature and humidity) during the transportation process [8].
Alternative approaches of determining the geopotential difference using clocks require connecting two clocks with optical fiber/coaxial cables or satellites to transmit frequency signals or time signals between two stations [9]- [13]. Chou et al. [14] observed time dilation by comparing two separate optical atomic clocks connected by a 75 m optical fiber, and they detected a change in height of 37 ± 15 cm compared to the actual height of 33 cm. Atomic clocks with a stability of 1.0 × 10 −18 may be capable of sensing height variations of one centimeter [15]. A time transfer simulation experiment was carried out using two atomic clocks with stabilities of 1.0 × 10 −18 connected by a coaxial cable, and the results of the simulation experiment achieved an accuracy of 0.16 m 2 /s 2 (equivalent to 1.6 cm in height) [16]. Although relative accuracies of up to the 10 −19 level can be achieved for comparisons of time and frequency transfer using optical fibers [17], [18], sufficient underground fibers are required to connect two arbitrary stations, which limits its applications, for instance over ocean and mountainous areas.
With the rapid development of the global navigation satellite system (GNSS), receiver clock offsets can be estimated as independent parameters for each measurement epoch [19], which could have potential applications in time-frequency science, providing a good opportunity to determine the geopotential difference. According to the studies of Shen et al. [8], [11], [20], the geopotential difference between a ground station and a satellite can be calculated using the satellite frequency signal transmission based on the Doppler-canceling technique or tri-frequency combination technique. Unlike time and frequency transfer using fibers on the ground, time and frequency transfer between ground stations and satellites poses much more challenging problems, such as signal propagation delay and frequency shift, orbital error, and Earth rotation. If a time and frequency transfer accuracy of less than 1 ns or higher is desired, the aforementioned errors must be eliminated or significantly diminished. As precise satellite orbits and clock products can be generated in the frame of the international GNSS service (IGS) [21], the precise point positioning (PPP) technique is widely applied to compute the time and frequency links with sub-nanosecond accuracy [22], [23]. With its low noise of phase measurements, PPP can provide much higher short-term stability than the GNSS pseudorange-only technique. Petit et al. [24], [25] demonstrated that the frequency transfer accuracy of the integer-PPP technique can reach 1.0 × 10 −16 within a few days for links, and the accuracy is in the low 10 −17 for averaging time above 10 days.
Considering the advantages of the high-precision of GNSS PPP time transfer technique and rapid development of time and frequency science, for instance at present optical atomic clocks achieve a stability of 4.8×10 −17 at 1 s and 6.6×10 −19 over an hour-long measurement [26], in this study, we propose an approach that uses the PPP technique to directly compute clock offsets between two clocks at two arbitrary positions for the determination of geopotential difference, referred to as the relativistic PPP time comparison approach, and the accuracy of this approach depends not only on the both accuracies and stabilities of clocks, but also the time transfer technique itself.
In this paper, we first introduce the relativistic PPP time comparison approach, including the relationship between the change in accumulated time difference and the geopotential difference between two remote clocks at two stations, and the principle of PPP time transfer technique. The clocks used should be previously calibrated at same site. Then the clocks run freely at different stations of our interests without any artificial adjustment. Based on this approach, the geopotential difference could be determined. To validate the proposed approach, simulation experiments are conducted. We simulated GNSS observations from ground stations equipped with ultra-high-performance clocks (say optical clock), and added various measurement errors to the GNSS observation simulation. The change in time difference between two clocks can be estimated from simulated observations using the PPP time transfer technique, and then the geopotential difference between two stations can be determined. Afterward, we described the simulation strategies of GNSS observations, receiver clock offset model, experimental data, and processing strategies. Thereafter, we comprehensively evaluate the performance of this approach in different cases using different types of clocks, and relevant results are presented. Furthermore, we discuss the requirements in this approach and the relationship between the performance of time links and the accuracy of this approach. Finally, we conclude this work in the last section.

II. METHODOLOGY
In this section, we first introduce the relativistic time comparison approach for determining the geopotential and then briefly describe the process of determining the change in time difference between two separate clocks using the PPP time transfer technique.

A. DETERMINATION OF GEOPOTENTIAL DIFFERENCE BETWEEN TWO GROUND STATIONS BASED ON CLOCKS' RUNNING RATES
According to GRT, an accurate clock runs faster at a position with higher geopotential than an identical clock at a position with lower geopotential [4]. Suppose there are two clocks C P and C Q at two different stations P and Q, respectively, and one clock C 0 on the geoid (an equi-geopotential surface nearest to the mean sea level) on which the geopotential constant is noted as W 0 , and the value of W 0 is 62 636 851.71 m 2 /s 2 [27], [28]. After a standard time duration dt 0 , recorded by clock C 0 , clocks C P and C Q will have time durations dt P and dt Q , respectively. Accurate to the level of c −2 , dt P and dt Q could be expressed respectively as [3], [29]: where W 0P = W P − W 0 and W 0Q = W Q − W 0 represent the geopotential numbers between the geoid and the stations P and Q, respectively; W P and W Q denote the geopotentials at stations P and Q; and c is the speed of light in vacuum. In this study, the definition of the geopotential in physical geodesy is applied: it always holds that W ≥ 0, which is different from the definition in physics [8]. Thus, if station P (or Q) is above the geoid, it holds that W P < W 0 (or W Q < W 0 ) [30]. Based on (1) and (2), after standard time duration T , the clock offset between clocks C P and C Q can be described as: where t P and t Q are the clock offset of clocks C P and C Q , respectively. Since the earth's gravity field is weak, the geopotential difference between P and Q can be formulated in approximate form as: where W PQ = W Q − W P ; O c −4 represent the order terms higher than c −2 , which can be omitted in this study. Hence, if GRT holds, we can determine the geopotential difference between the two stations by comparing the change in time difference t PQ after standard time duration T using precise clocks. Inversely, given known two datum points on ground, we may test GRT. Suppose the orthometric height of P (noted as H P ) is given, the orthometric height of Q (noted as H Q ) can be determined, expressed as: whereḡ P andḡ Q are the mean value of the gravity along the plumb line, respectively. Ifḡ i (i = P, Q) can be determined, then the H i can be computed. The mean gravityḡ i can be formulated as:ḡ where g(h) is the actual gravity at the station i which has the height h. Since we cannot determine the mean gravityḡ i precisely, (6) can be approximated as [31]: where g i , in gal, is the gravity measured at the ground station i, which can be measured by absolute gravimeter, and H i in km. The factor 0.0424 refers to the normal density ρ = 2.67 g/cm 3 . According to (5) and (7), we obtain the orthometric height H Q , then the practically useful formula can be expressed as: where W PQ is measured in geopotential units (g.p.u.), and 1 g.p.u.=1000 gal.m. We can observe that the accuracy of the determined W PQ will affect that of H Q , then the determination of geopotential difference has the potential applications in orthometric height determination and in unifying the world vertical height system (WVHS).
The key problem is to determine the change in time difference between two remote clocks using the GNSS satellites microwave signal. In order to eliminate the errors involved in GNSS time transfer, here we use the GPS ionosphere-free combination PPP model (IF-PPP).

B. GPS IONOSPHERE-FREE PPP MODELS
The undifferenced (UD) observation equations of the dual-frequency pseudorange P s r,i and carrier phase L s r,i measurements at one epoch can be respectively modeled as [32]: where indices s, r, and i denote the satellite, receiver, and carrier frequency, respectively; P s r,i and L s r,i denote the pseudorange and carrier phase measurements, respectively, in meters, which need to be simulated at each epoch; ρ s r represents the geometric distance between the phase centers of the satellite s and receiver r antennas; c is the speed of light in vacuum; dt r and dt s denote respectively the clock offsets between the receiver clock and GPS time (GPST) and the satellite clock and GPST; T s r is the tropospheric delay of the signal path in meters; I s r,1 is the slant ionospheric delay on frequency f 1 ; γ s i is the ionospheric factor depending on the frequency r,i is the integer phase ambiguity in cycles; d r,i and d s i are the code delays of receiver and satellite, respectively, in meters; b r,i and b s i are the uncalibrated phase delays (UPDs) at receiver and satellite, respectively in cycles; and ε s r,i and ξ s r,i refer to the multipath effects and unmodeled measurement errors for pseudorange and carrier phase observations, respectively in meters. In addition, some errors, such as the dry component of tropospheric delays, tidal loading, phase wind-up, relativistic effects in the satellite clock, phase center offsets (PCOs), and phase center variations (PCVs), are not mentioned in (9) and (10). These errors can be also corrected by relevant models [27], [33].
The IF-PPP model is most generally used for time and frequency transfer, which can eliminate first-order ionospheric delays [34], expressed as: where α and β are ionosphere-free combination coefficients with and L s r,IF refer to ionosphere-free pseudorange and carrier phase observations, respectively, in meters; dt r and dt s are respectively receiver and satellite clock offsets, which absorb ionosphere-free code hardware delays.

C. DETERMINATION OF GEOPOTENTIAL DIFFERENCE BETWEEN TWO GROUND STATIONS VIA GNSS SATELLITES
As shown in Fig. 1, there are two ground stations P and Q, which can receive the signals emitted from GNSS satellites.
To determine the geopotential difference between two stations using the relativistic PPP time comparison approach, we performed the following procedures. (i) The freely running clocks C P and C Q are used at stations P and Q, respectively; the clocks should be a priori synchronized at same site (say at station P) before the experiment, and their vibration frequencies should not be adjusted during the whole experiment. Hence, the clock's stability is significant in this study. (ii) After accurate synchroniation at station P, the clock C P is fixed at P and clock C Q is transported to station Q, both clocks outputing local 1 pulse per second (1 PPS) signals, and the GNSS receivers linked with precise clocks obtain the observations from GNSS satellites signals. (iii) We employed an open-source program RTKLIB [35] which can perform GNSS IF-PPP processing to compute the time offsets between the clock C P (or C Q ) at a ground station and the GPST, where the GPST is taken as reference. By subtracting one time offset series from another we obtain the change in time difference between the remote clocks C P and C Q on ground, cancelling the ''common'' GPST. Here we see that the GPST is used only as reference (or ''bridge''), having no effect on the time offsets between two remote clocks. (iv) The corresponding geopotential difference between stations P and Q is determined based on (4).

III. DATA SIMULATION STRATEGIES
Although some ground stations are equipped with highperformance clocks, their clocks are frequently adjusted to conform with the Coordinated Universal Time (UTC). Hence, the observations at present IGS stations cannot be used to determine the geopotential. In another aspect, at present the stability of the GNSS PPP frequency transfer is limited to about 1.0 × 10 −15 at 1 day averaging time [24], which is insufficient to validate the relativistic PPP time comparison approach at the decimeter level. These difficulties could be overcome by using free-running clocks with ultra-high FIGURE 1. Principle of determining geopotential difference between stations P and Q via GNSS satellite signal transmission. Two clocks are synchronized accurately at station P, and one clock is transported to station Q after synchronization. Both clocks outputting local 1 PPS signals, and receivers receiving the signals from GNSS satellites to extract the observations. The background is based on the earth relief dataset provided with Generic Mapping Tools [36].
accuracy and stability. Here in this study, we simulated GNSS observations from ground stations equipped with ultra-highperformance clocks (say optical clocks). GNSS observations associated with ground station clocks were simulated to evaluate the performance of this approach, and the flowchart of determining the geopotential difference between two ground stations using this approach is shown in Fig. 2.

A. SIMULATION OF GNSS OBSERVATIONS CONNECTED WITH GROUND STATIONS CLOCKS
In this subsection, we describe the simulation of GNSS observations, which is essentially the inverse process of positioning. As mentioned by Section II-B, all measurement errors, including satellite-dependent, receiver-dependent, and atmospheric propagation errors, were corrected according to existing models [33]. GNSS observations were used to estimate the receiver position and clock offset in the GNSS positioning process. For simulating GNSS observations, if the receiver/satellite position and clock are known values, the satellite-receiver geometric distance can be computed. To simulate the observations as accurately as possible, various measurement errors were calculated using existing models and added to the geometric distance with random noise to generate the observations [37].
As shown in Fig. 2, we focused on simulating all components on the right side of (9) and (10) in the observation simulation process. The satellite-receiver geometric range is calculated by the positions of the GNSS satellite and ground receiver. The positions of the GNSS satellite can be extracted by final GNSS orbit products with a sampling rate of 5 min, acquired from the Center for Orbit Determination in Europe (CODE). Nevertheless, satellite positions need to be interpolated to access the satellite positions at any epoch. In addition, receiver positions can be obtained from weekly solutions of IGS stations. Tropospheric delay is composed of dry and wet components, both of which can be formulated as zenith delay and the corresponding mapping function. In terms of zenith dry delay, which can be calculated by combining the global pressure and temperature (GPT) empirical model with the Global Mapping Function (GMF) [38], [39], according to the Saastamoinen model [40]. Zenith wet delay can be obtained from GNSS PPP solutions. Ionospheric delay is associated with the total electron content (TEC), which can be obtained from the global ionosphere maps (GIMs) of CODE. The initial ambiguity of phase measurement is set to an integer number of cycle constant for each continuous arc, and cycle slips are also introduced at some epochs for certain satellites. The satellite clock offsets are obtained from CODE precise clock products, and the receiver clock offsets are simulated by stochastic differential equations (SDEs) [41], which are described in the next subsection. Moreover, some errors, such as the ocean tide loading displacement, PCOs, and PCVs, need to be included in the observation errors to ensure that the simulated observation data are as close to real data as possible. The ocean tide loading displacement can be calculated by the FES2004 ocean tide model [42]. Antenna products from IGS are employed in PCOs and PCVs simulations for satellites and receivers [27]. To generate realistic data, the observation noises are simulated as a zero-mean-value random measure noises with Gaussian distribution, for which the standard deviation (STD) depends on satellite elevation angle. Specifically, STD decreases with increasing elevation angle. In the zenith direction, the STD of each frequency is set to 0.2 m and 2 mm for code and phase observations, respectively.

B. MATHEMATICAL MODEL OF CLOCK OFFSET
In this subsection, we introduce the mathematical model of the clock offset, which is very useful for clock simulation. Clock noises are typically affected by five random noises [43] generally known as white phase modulation (WPM), flicker phase modulation (FPM), white frequency modulation (WFM), flicker frequency modulation (FFM), and random walk frequency modulation (RWFM). These random noises vary according to the type of clock. Previous studies have reported that WFM and RWFM are the predominant noises in high-precise atomic clocks [41], [44]. Accordingly, the mathematical model of the clock offset can be formulated as [41]: , t ≥ 0 (12) with initial conditions where X 1 (t) represents the clock phase deviation, and X 2 (t) is a part of the clock frequency deviation, which is generally called the random walk component; the constants µ 1 and µ 2 can be interpreted as drift terms for the two Wiener processes, referred to as the deterministic component driving clock errors in general. In particular, µ 1 is related to the constant initial frequency offset, indicated by y 0 ; µ 2 is generally indicated by d, and it is the frequency drift, also termed as aging. In more familiar metrological notation, c 1 + µ 1 = y 0 and µ 2 = d. W 1 (t) is the Wiener process on the clock phase deviation X 1 (t) driven by the white noise of frequency and corresponds to WFM, whereas W 1 (t) denotes the Wiener noise of frequency corresponding to RWFM, which produces an integrated Wiener process on the phase. The constants σ 1 and σ 2 are the diffusion coefficients of the two noise components W 1 and W 2 , respectively, and indicate the intensity of WFM and RWFM. We note that clock offset is dominantly affected by WFM, RWFM, and the deterministic trend terms µ 1 and µ 2 . The relationship between the Allan deviation (ADEV) and the diffusion coefficients of the SDEs can be expressed as [44]: where σ WFM y (τ ) and σ RWFM y (τ ) denote the ADEVs related to the WFM and RWFM at averaging time τ , respectively.
To evaluate the performances of the relativistic PPP time comparison approach for different types of receiver clocks, we propose three cases as summarized in Table 1. The deterministic component µ 1 is considered as the gravitational frequency shift t/T , which can be solved by (4). As the value of frequency drift d does not affect the uncertainty of the frequency drift [44], [45], we can simply assume that d = 0. We simulated the clock in Case 1 with WFM having a fractional frequency stability of σ WFM y (τ ) = 4.0×10 −13 / √ τ and σ RWFM y (τ ) = 4.0 × 10 −19 √ τ , both with averaging time τ , which is typical Cs-fountain performance, and the VOLUME 8, 2020 accuracy of this clock has been evaluated to be 5.3 × 10 −16 [46], [47]. The clocks in Cases 2 and 3 were simulated with WFM having fractional frequency stabilities of σ τ at averaging time τ , and the accuracy of these clocks can reach 7.8 × 10 −17 and 8.6 × 10 −18 , which are typical for optical clocks [48], [49]. It should be noted that the simulation parameters of clocks in Case 1-3 are applied for whole experiment period, without artificially steering outputs of clocks to approximate UTC.

C. EXPERIMENTAL DATA AND PROCESSING STRATEGIES
To evaluate the performance of the PPP time transfer technique and validate the tests of geopotential determination via relativistic PPP time comparison approach, we selected three IGS stations that can provide fixed station coordinates from IGS weekly solutions for ground-based observation simulation. The geographical distribution of the selected stations is shown in Fig. 3, but the GNSS observation data of these stations from IGS are not used in this experiment. Table 2 summarizes the details of the selected IGS stations, including the latitude, longitude, height (here it denotes the height above geoid), and geopotential (W ). The geopotentials at all selected stations were calculated using the new global combined gravity field model EIGEN-6C4. The accuracy of the EIGEN-6C4 model is approximately 10-20 cm in land areas [50], which is sufficient for the precision requirement of this study. The experimental data cover a 30-day period for the day of year (DOY) 061-090 in 2020 (during Modified Julian Date (MJD) 58909-58938).  In addition, the detailed PPP processing strategy is summarized in Table 3, including the estimated parameters and observation models. Using the simulated GPS observations, we can analyze performances of both the PPP time transfer technique and geopotential difference determination with different receiver clock cases. In the data processing, GPS PPP solutions were based on dual-frequency (L1/L2) ionosphere-free combinations for eliminating first-order ionospheric delay. The initial standard deviations for the raw GPS carrier phase and code observations were 2 mm and 0.2 m, respectively. Precise satellite orbit and clock products are provided by CODE at sampling rates of 5 min and 5 s, respectively.

IV. RESULTS
We first evaluated the performances of the clocks in different cases. With the simulated GNSS observation data, we analyzed the performance of the PPP time transfer technique with different cases, and the results of the relativistic PPP time comparison approach are presented.

A. PERFORMANCE EVALUATION OF CLOCK OFFSETS SIMULATION
In order to analyze the performances of the relativistic PPP time comparison approach with different receiver clocks, GNSS observations need to be simulated with different receiver clock offsets, as summarized in Table 2. Based on the SDEs, we simulated the receiver clock offsets over a 30-day period and evaluated the performance. The frequency stabilities of the simulated clock offsets, shown in Fig. 4, included total Allan deviation (TADEV), along with WFM asymptotes for Cases 1-3. The WFM asymptotes for Cases 1-3 follow 4.0 × 10 −13 / √ τ , 2.3 × 10 −14 / √ τ , and 2.8 × 10 −15 / √ τ at every averaging time τ , respectively, which are consistent with simulation parameters of WFM.
To further evaluate the performances of three cases for simulating clock offset, we simulated 500 independent clock offsets for each case based on the SDEs. Fig. 5 presents the 500 clock offsets, and the 1-σ uncertainty (68% confidence level) and 2-σ uncertainty (95% confidence level) of the 500 clock offsets for a 30-day period. For this period, the 1-σ uncertainties of Cases 1-3 were approximately 646.25 ps, 36.97 ps, and 4.44 ps, respectively, which correspond to geopotential difference determination errors of ± 22.41 m 2 /s 2 , ± 1.28 m 2 /s 2 , and ± 0.15 m 2 /s 2 . Therefore, considering Figs. 4 and 5, it can be concluded that the instability of the clocks decreases with increasing averaging time, and the uncertainty of clock offsets increases with experiment  time. Nevertheless, both the instability and uncertainty of the clocks can be improved by increasing the stability of the clocks.

B. ACCURACY ASSESSMENT OF PPP TIME TRANSFER IN DIFFERENT CASES
We evaluated the performance of PPP time transfer in different cases against true values. True values mean that σ WFM y and σ RWFM y are set to 0. The performance of BRUX station in Cases 1-3 between MJD 58909-58938 is presented in Fig. 6. Fig. 6 (a) shows that the clock offset deviations of Case 1 are greater than those of Cases 2 and 3. Moreover, because the accuracy of clock offset is affected by measurement noise, the clock offset deviations of Case 2 are consistent with those of Case 3. Deviations in Case 1 mostly ranged from −0.5 ns to 0.1 ns, but those in Cases 2 and 3 remained within approximately 0.1 ns for the entire period. Fig. 6 (b) presents the frequency stability of BRUX in Cases 1-3. We observed that only the short-term stabilities (within averaging of 1000 s) of Cases 1-3 were very close to each other; the stabilities of Cases 2 and 3 were higher than that of Case 1 for averaging time at and above 1000 s.
To further evaluate the performance of the three cases, we calculated the STD and root mean square (RMS) of clock offsets of Cases 1-3 with respect to the true values for the three stations, and the results are presented in Fig. 7. The accuracy of the clock offsets in terms of STD obtained from Case 1 ranged from 0.1 ns to 0.3 ns, for the three stations, whereas the RMS ranged from 0.2 ns to 0.5 ns. The performance of Case 3 was the best among the three cases, but the STD and RMS of Cases 2 and 3 were very similar for the three stations. The STDs of both Cases 2 and 3 were within approximately 50 ps and the RMS ranged from 40 ps to 100 ps.
We used Cases 1-3 to respectively compute the clock offsets of two time transfer links, for which BRUX was  selected as the reference clock. The performance of time links in Cases 1-3 was evaluated. Fig. 8 shows that the clock offsets of time links obtained from Cases 1-3 with respect to the true values, and Fig. 9 shows the frequency stability of two time links. In addition, Table 4 summarizes the STD and RMS values of clock offsets for each link. Overall, three conclusions can be drawn from the results. First, the clock offsets of both Cases 2 and 3 are fairly consistent with the true values, and the clock offset deviations of Cases 2 and 3 varied from −0.1 ns to 0.1 ns, for both OPMT-BRUX and PTBB-BRUX. In contrast, the deviations of Case 1 of OPMT-BRUX and PTBB-BRUX ranged from −0.5 ns to 0.6 ns and from −0.1 ns to 1 ns, respectively. Second, the frequency stabilities of Cases 2 and 3 for each link were higher than those for Case 1 at averaging time of 1000 s and above. The frequency stabilities of Cases 1-3 were approximately 4.28 × 10 −16 , 4.00 × 10 −17 , and 3.22 × 10 −17 , respectively, at 10-day averaging for OPMT-BRUX. For PTBB-BRUX, these values were approximately 3.73 × 10 −16 , 8.17 × 10 −17 , and 4.64 × 10 −17 . According to previous studies [3], [15], a clock with a stability of 1.0 × 10 −16 may sense one-meter height variations, according to GRT. The deviation errors of geopotential difference were approximately 38.50 m 2 /s 2 , 3.60 m 2 /s 2 , and 2.90 m 2 /s 2 for OPMT-BRUX in Cases 1-3, respectively, after a 10-day test period; for PTBB-BRUX, the errors were approximately 33.56 m 2 /s 2 , 7.34 m 2 /s 2 , and 4.17 m 2 /s 2 . Third, the maximum values of STD and RMS for each link were observed in Case 1, and the minimum values were observed in Case 3. Therefore, we can conclude that the   clock offsets of stations and links obtained from Case 3 are more accurate and stable than those obtained from Case 1, Cases 2 and 3 can provide much higher long-term stability owing to the higher performance clocks, and the frequency stabilities of Cases 2 and 3 facilitate the determination of geopotential difference at the decimeter-level after 10 days of observations.

C. VALIDATION OF RELATIVISTIC PPP TIME COMPARISON APPROACH
To determine the geopotential difference between two stations based on GRT, the change in time difference between the clocks at two stations during the test period need to be calculated. For this purpose, we obtained clock offsets via PPP time transfer with Cases 1-3 and analyzed performance of clock offsets (Section IV-B). We conducted a geopotential difference experiment to test the validity of the relativistic PPP time comparison approach. From the beginning of the experiment, we measured the geopotential difference every day to evaluate the performance of the experimental results with increasing time. The results of the determined geopotential difference deviating from those of the EIGEN-6C4 model are shown in Fig. 10 and Table 5. Error bars represent the 1-σ uncertainty obtained from the TADEV of time links, as shown in Fig. 9. We performed an experiment using 30-day observations, and TADEV was calculated for about one third of the whole period, namely 10-day period at an average time of 10 6 s. Frequency stability after the averaging time of 10 6 s can be obtained by fitting the values of TADEV using the least squares method, and extrapolating the fit to the entire experimental period [15]. The geopotential difference determination errors caused by the frequency stabilities of Cases 1-3, which are larger than the corresponding deviation errors, can also be used to measure the uncertainty of this approach [24]. As shown in Fig. 10, we observed that the deviation errors of geopotential differences between OPMT and BRUX and between PTBB and BRUX decrease with increasing time in all cases. In general, the largest errors of geopotential differences as well as the largest uncertainties were obtained at the beginning of the experiment because of errors in the observation data and insufficient data for estimation. In addition, the deviation errors became stable and close to zero for Cases 2 and 3 after ten days of experimental time, but the errors for Case 1 continued to fluctuate due to the limited clock accuracy and measurement noise. Table 5 presents the results of the experiment for 1-day, 5-day, 10day, 20-day, and 30-day test periods. From these results, three conclusions can be directly drawn. First, because each ground station is equipped with a clock with the same noise level and the same parameters were applied for simulating observation data, the deviation errors and uncertainty of each link are similar in each case. Second, the deviation error decrease as clock accuracy increases. Case 3 shows the highest performance, whereas Case 1 has the largest deviation error as well as uncertainty. In Case 3, the deviation error of measured geopotential differences between OPMT and BRUX was −0.64 m 2 /s 2 with an uncertainty of 1.11 m 2 /s 2 , and the deviation error between PTBB and BRUX was 0.76 m 2 /s 2 with an uncertainty of 1.79 m 2 /s 2 . It is noteworthy that although the deviation error using the 20-day observations is smaller than that using the 30-day observations, the uncertainty of the former is higher (see Table 5). Third, the accuracy of geopotential difference achieved using the relativistic PPP time comparison approach is consistent with the stabilities of the time links in Cases 1-3. Overall, the ground station equipped with a clock of higher accuracy not only determines the geopotential difference between two stations at the 1 m 2 /s 2 level, but also improves the uncertainty of the deviation error. In other words, relativistic PPP time comparison for determining the geopotential difference between two stations is feasible.

V. DISCUSSION AND CONCLUSION
This study focuses on high-performance time comparison between two arbitrary ground stations for determining the geopotential difference using the relativistic PPP time comparison approach, and comprehensively validating this approach with different receiver clocks. For this purpose, two key requirements need to be satisfied. First, all stations of interests should be equipped with high-performance clocks to maintain an accuracy and stable frequency standard without artificially adjustments. Second, the clocks used at different stations should be a priori synchronized at one same site. However, it should be acknowledged that the frequency stabilities of the clocks at IGS stations at present are limited to about 1.0 × 10 −15 at 1-day averaging time. Moreover, a priori synchronization information of two clocks at different IGS stations is not available, and the uncertainty of synchronization affects the accuracy of the experimental results.
The key point is that the frequency of the clocks at IGS stations is steered to stay aligned on UTC, so it can certainly not be used to determine the geopotential. Consequently, it is difficult to validate the relativistic PPP time comparison approach at the decimeter level using real GNSS observations at present IGS stations. To solve this problem, we need GNSS observations connected with ultra-high-performance atomic clocks which freely run, without being steered to any standard reference. As a prelude, in this study we simulate GNSS observations at ground stations equipped with ultra-highperformance atomic clocks without any kinds of adjustments.
We choose three stations and simulate GNSS observations at these stations equipped with atomic clocks with high accuracies and stabilities to assess the performance of the relativistic PPP time comparison approach over several durations of up to 30 days. The validation is performed following a three-step procedure: (i) simulating GNSS observations, (ii) evaluating the performance of PPP time transfer, and (iii) determining the geopotential difference between two stations. In the first step, we simulated the receiver clock offsets in the three cases based on SDEs, and employed them in the observation simulation. We also analyzed the performance of the clock offsets (Section IV-A); the accuracy and uncertainty of clock offsets are the part of errors involved in determining geopotential difference. In the second step, the PPP technique is used to compute the clock offsets and the performance is evaluated (Section IV-B). Combining Figs. 6-9 and Table 4, the performance of clock offsets in Case 3 is found to be better than that in Case 1. Regarding frequency stability, results of Cases 2 and 3 improve the long-term stability of OPMT-BRUX as well as PTBB-BRUX. Comparative analyses show that the precision of PPP time transfer can be improved by using more accurate clocks, but the improvement is also limited by various noises in observations. Compared to Case 2, the performance of Case 3 is only slightly improved. Based on the performance evaluation of PPP time transfer, it is determined the geopotential difference between two stations during the 30-day period in the last step. With the requirement of maintaining continuous operation of the clocks, deviation errors and uncertainties of time links decrease with increasing experiment time, providing a geopotential difference of high accuracy. We may conclude that the accuracy of this approach is consistent with the stabilities of the clocks used for each case; that is, this approach is feasible for determining the geopotential difference between two stations.
The results of this study show that the relativistic PPP time comparison approach has the potential to achieve decimeter-level accuracy. As for this approach, it is not only a new development trend of geopotential determination, but also has important application in various fields, for instance geodynamics, astronomy, and aerospace. However, this approach is mainly limited by the performance of atomic clock, the accuracy of GNSS satellite orbit and clock products as well as GNSS measurement noise. Although GNSS orbit and clock products and other real data were employed to simulation experiments, and carefully considered the performance of atomic clock as well as the measurement errors correction, some abnormal situations (e.g. unhealthy of GNSS satellites, atmosphere anomaly, and ground environment anomaly) may still occur in actual observations. Hence, the performance of this approach in practical applications will be further studied. With the rapid development of GNSS, the accuracy of orbit and clock products as well as the models of various measurement errors will be improved. The determination of geopotential difference between two stations at the centimeter level will be further investigated. The formulation of this study could be also applied to testing GRT.
ZIYU SHEN (Member, IEEE) is currently a Lecturer with the School of Resources, Environmental Science and Engineering, Hubei University of Science and Technology, Xianning, China. His current research interests include relativistic geodesy and satellite frequency transfer technique.
WEI XU received the master's degree from the Anhui University of Science and Technology, in 2018. He is currently pursuing the Ph.D. degree with Wuhan University. His current research interests include multi-frequency GNSS data processing and its applications.