Benefits of combining GPS and GLONASS for measuring ocean tide loading displacement

Abstract

GPS has been used to estimate ocean tide loading (OTL) height displacement amplitudes to accuracies of within 0.5 mm at the M2 frequency, but such estimation has been problematic at luni-solar K2 and K1 frequencies because they coincide with the GPS orbital period and revisit period, leading to repeating multipath and satellite orbit errors. We therefore investigate the potential of using the GLONASS constellation (with orbital period 11.26 h and true site revisit period of 8 sidereal days distinct from K2 and K1) for OTL displacement estimation, analysing 3–7 years of GPS and GLONASS data from 49 globally distributed stations. Using the PANDA software in kinematic precise point positioning mode with float ambiguities, we demonstrate that GLONASS can estimate OTL height displacement at the M2, N2, O1 and Q1 lunar frequencies with similar accuracy to GPS: 95th percentile agreements of 0.6–1.3 mm between estimated and FES2014b ocean tide model displacements. At the K2 and K1 luni-solar frequencies, 95th percentile agreements between GPS estimates and model values of 3.9–4.4 mm improved to 2.0–2.8 mm using GLONASS-only solutions. A combined GPS+GLONASS float solution improves accuracy of the lunar OTL constituents and P1 (but not significantly for K1 or K2) compared with a single-constellation solution and results in hourly-to-weekly spectral noise very similar to a GPS ambiguity-fixed solution, but without needing uncalibrated phase delay information. GLONASS estimates are more accurate at higher compared with lower latitudes because of improved satellite visibility, although this can be countered by using a lower elevation cut-off angle.

1 Introduction

Geodetic measurements, for example from Global Navigation Satellite Systems (GNSS), Very Long Baseline Interferometry (VLBI), Satellite Laser Ranging (SLR) and Doppler Orbitography Integrated by Satellite (DORIS), are sensitive to ocean tide loading (OTL) deformation of the solid Earth which is caused by the periodic change in ocean mass distribution arising from the gravitational attractions of the moon and Sun. The IERS Conventions (Petit and Luzum 2010) provide utilities to correct geodetic measurements for this OTL deformation, requiring as input OTL displacement coefficients at the dominant tidal periods (including those listed in Table 1). These are generated by convolving a global model of the ocean tides with a loading Green’s function, which is dependent on the material properties of the Earth’s interior. Because any errors in these OTL displacement coefficients will propagate to the normally estimated geodetic parameters and degrade, for example, resultant coordinate time series used for reference frame definition and the monitoring of millimetre-level land movements, it is important that accurate models are used for their derivation. One way in which the accuracy of these Earth and numerical ocean tide models can be verified is by independent geodetic analysis in which the OTL displacements are the parameters of interest.

Table 1 Principal semi-diurnal and diurnal tidal constituents of the tidal potential (after Kudryavtsev 2004)

VLBI data were first shown to be able to estimate OTL displacement by Schuh and Moehlmann (1989) and then Sovers (1994), who included harmonic parameters at the dominant tidal frequencies in the primary least squares estimation. Schenewerk et al. (2001) showed this was also possible with global solutions of double differenced Global Positioning System (GPS) data but to an accuracy of ~ 5 mm for 90% of the sites studied, whereas Allinson et al. (2004) used precise point positioning (PPP) for at least 90 days of GPS data to obtain M2 OTL displacement agreements with geophysical models within ~ 1 mm. Thereafter King et al. (2005) used the PPP GPS method of Allinson et al. (2004) to obtain OTL displacements to validate ocean tide models around Antarctica. Thomas et al. (2007) compared VLBI and PPP GPS analyses (each using several years of data) and concluded similar millimetre-level agreement for GPS and VLBI when compared with OTL computed from existing Earth and ocean tide models, for the majority of tidal constituents. An alternative approach was followed by Khan and Tscherning (2001) and Melachroinos et al. (2008), who undertook harmonic analysis of GPS coordinate time series to estimate the OTL displacement, obtaining observed versus model differences of several millimetres but using only 7–15 weeks of GPS data. Penna et al. (2015) refined this time series analysis technique by determining the optimum tropospheric and coordinate process noise through comparisons with radiosonde tropospheric delays and synthetic harmonic ground displacements. This led to the estimation of OTL displacement using GPS to an accuracy of around 0.4 mm when using time series from 2.5 years of data, improving to about 0.2 mm with 4 years or more of data. While ocean tide model errors have historically been assumed to be the limiting accuracy factor in the modelling of OTL displacement (e.g. Bos and Baker 2005), recent advances in ocean tide modelling (e.g. Stammer et al. 2014) have led to GPS-estimated OTL displacements being used to not only validate and identify deficiencies in ocean tide models, but also to measure the elastic and anelastic properties of the Earth’s interior (e.g. Ito and Simons 2011; Bos et al. 2015).

Studies to date on probing the Earth’s interior properties at tidal frequencies using GPS have mostly considered the M2 constituent only (e.g. Bos et al. 2015; Yuan and Chao 2012; Martens et al. 2016), and validation of ocean tide models using GPS-estimated OTL displacements has proved especially problematic at the K2 and K1 frequencies (e.g. Allinson et al. 2004; King et al. 2005; Thomas et al. 2007). This is because K2 and K1 coincide with the GPS orbital period and sidereal geometry repeat period, respectively, so any orbit errors and multipath effects degrade the OTL displacement estimates even over time spans of several years (e.g. Thomas et al. 2007). Post-processing periodic error mitigation techniques, e.g. sidereal filtering for multipath, would inadvertently remove part of the OTL and hence cannot help overcome the GPS problem. The completion of the GLONASS satellite constellation replenishment in 2010, the subsequent upgrade of networks of GNSS receivers worldwide such that there are now over five years of both GLONASS and GPS dual frequency observations widely available, together with IGS Analysis Centres generating high accuracy GPS and GLONASS satellite orbits and high-rate clocks, now facilitate the estimation of OTL displacement using GLONASS. This is particularly desirable, as the GLONASS orbital period of 11.26 h (~ 2.131 cycle per day) and the sidereal geometry repeat period of 8 days (0.125 cycle per day) are distinct from any major tidal frequencies, so K2 and K1 OTL displacement estimation becomes potentially feasible. This complements the promise shown by GLONASS for longer period crustal deformation studies, with Abraha et al. (2018) demonstrating that GLONASS can result in reductions compared with GPS in artificial longer period signals arising from the propagation of unmodelled semi-diurnal and diurnal tidal displacements, because of the different geometry repeat period.

This paper investigates how well OTL displacement may be estimated using GLONASS observations, in particular for the GPS-problematic K2 and K1 frequencies. We also investigate whether combining GPS and GLONASS observations can lead to more accurate OTL displacement estimation than when using either GPS or GLONASS observations alone. A globally distributed set of GPS+GLONASS continuous receiver data spanning at least three years at carefully selected stations is used, with validation undertaken by comparison with forward geophysical model OTL displacement values. We focus on the height component, as these OTL displacements are typically three times the size of the horizontal components (Baker 1984).

2 OTL displacement estimation using multi-GNSS kinematic PPP

As described by Penna et al. (2015), OTL displacement can be estimated by the GNSS precise point positioning (PPP) technique in two ways, which they termed harmonic estimation and the kinematic approach. In kinematic PPP (which we use here, following Penna et al. 2015), satellite positions and clock offsets are held fixed, and a variety of systematic errors, including antenna phase centre variations (PCV), phase windup, atmospheric propagation effects and tidal displacements, are corrected. Then, parameters of interest are estimated which include time-varying 3D station coordinates, receiver clock offsets (at each data collection epoch), unmodelled time-varying tropospheric delays and phase biases for each satellite during each phase-connected arc. Thereafter, the station coordinates in each processing session, e.g. 24 h, are concatenated to form a time series and then screened to remove blunders. In addition, a low pass filter in the form of a window average may be used to eliminate time series noise with periods much shorter than the diurnal and semi-diurnal tidal bands. Finally, by least squares spectral analysis of the time series for each desired coordinate component and tidal constituent, the amplitude and phase lag of the tidal displacement signals are estimated.

Here, for the kinematic PPP data processing, we use the Position and Navigation Data Analyst (PANDA) software (Liu and Ge 2003), as it not only has a proven capability in kinematic PPP combined multi-GNSS processing (e.g. Penna et al. 2018) but also allows the processing of either GPS or GLONASS data separately. To fix the satellite positions and their clock offsets at each data collection epoch, we use the ESA final (operational) products as they have the longest continuous record of high-rate (30-second) GPS+GLONASS satellite clock availability (2010 onwards) of all the IGS Analysis Centres, and they are of high quality throughout this interval, ~ 1 cm weighted root mean square difference from the IGS combined solution (Villiger and Dach 2018). We use EOP information provided by IERS bulletins as fixed values for each daily kinematic PPP batch. We apply IGS receiver and satellite antenna phase centre variation models, and we use the ionosphere-free combination of dual frequency data to mitigate ionospheric effects. The predictable parts of the tropospheric delay and tidal displacement (including the perturbation due to the free core nutation) are removed from GNSS observations according to the IERS Conventions 2010 (Petit and Luzum 2010), using the Saastamoinen (1972) formula and the Global Mapping Function (GMF) (Boehm et al. 2006) to reduce the hydrostatic and wet tropospheric delays. Due to their simpler modelling approach (e.g. Mathews et al. 1997), Earth body tide calculations are typically performed within PPP software packages. However, OTL displacement computations, which require information for the ocean tidal height and coastline geometry (e.g. Farrell 1972; Baker 1984), require a separate computation procedure. We input the FES2014b ocean tide model (Carrère et al. 2016), and an elastic Earth response Green’s function computed from the Preliminary Reference Earth Model (PREM) (Dziewonski and Anderson 1981) by Wang et al. (2020), to the NLOADF (SPOTL) software (Agnew 1997; 2012) to compute a priori OTL displacements. These are then applied in the GNSS processing via the hardisp program of the IERS Conventions 2010. Thus, the tidal displacement signals we estimate are residuals to this a priori model. We focused our study on the M2, N2, K2, K1, O1, P1 and Q1 constituents which typically have the largest semi-diurnal and diurnal OTL displacements (Table 1). We disregarded S2 OTL despite its typically large magnitude, as GNSS observations are also affected by S2 atmospheric loading displacement (e.g. Tregoning and van Dam 2005) and these two physical signals cannot be separated in the frequency domain. Because ESA satellite orbit/clock information is provided in the centre of GNSS network (CN) reference frame, which is a realisation of the Earth’s centre of figure (CF) frame (e.g. Dong et al. 2003), we compute the predicted OTL displacement with respect to the centre of mass of the solid Earth (CE), which closely resembles CF.

We adopt a dynamic model for the estimated time-varying kinematic PPP parameters consisting of white process noise for receiver clock offsets, and random walk process noise for the station coordinates and tropospheric zenith wet delay (ZWD) and its northward and eastward horizontal gradients. As described below, we use the method of Penna et al. (2015) to tune appropriate process noise values. The unknown parameters in the GPS-only solutions, namely 3D station coordinates and receiver clock corrections every 30 s, ZWD every 30 min and its horizontal gradients every 60 min, and a real-valued phase bias for each phase-connected arc, are estimated in a recursive least squares adjustment, over 24-hour sessions (chosen to minimise any additional errors from day break effects when concatenating 24-hour ESA orbits and clocks). The GLONASS-only and GPS+GLONASS solutions were parameterised as for GPS except: for GLONASS, also a time-constant inter-frequency bias is estimated per satellite (except for a reference satellite); for GPS+GLONASS, also a time-constant inter-system time and inter-frequency bias is estimated per satellite. As precise satellite orbit products for GPS and GLONASS are both computed in the International Terrestrial Reference Frame (ITRF), a coordinate frame transformation between them is not required in combined GPS+GLONASS PPP. Gross outliers were removed from the resulting 30-second detrended height coordinate time series if they were more than ten times the median absolute deviation from the median, before coordinate averaging in 30-minute bins, which were then used to estimate harmonic displacements at specific, defined tidal frequencies using least squares.

2.1 GNSS data selection

The GNSS stations used to assess the benefit of GLONASS for OTL displacement estimation were selected according to (i) GNSS data availability and quality, and (ii) the suitability of using forward geophysical models for the validation of the estimated OTL displacements. As the quality of a kinematic PPP solution will be directly dependent on the quality of satellite orbits and high-rate clocks, as well as how well ionospheric and tropospheric effects are mitigated, we selected a globally distributed set of GNSS stations to assess the impact of these effects on the OTL displacement estimation. In total, 49 globally distributed GNSS stations were selected (shown in Fig. 1 and listed in “Appendix A” along with the data availability and spans used) which fulfilled the criteria now described.

Fig. 1

RMS agreement of the magnitudes of the vector differences for the predicted OTL height displacement (in mm) per cell of a 0.25° global grid for the M2 (top) and K1 (bottom) constituents based on seven recent ocean tide models (FES2014b, GOT4.10c, TPXO8-Atlas, NAO.99b, HAMTIDE11a, DTU10 and EOT11a), with GNSS stations used in this paper shown as dark blue dots. The colour scale saturates at 2 mm RMS (maximum RMS for both M2 and K1 is 3 mm). The five stations with the largest (0.7–0.8 mm) inter-model disagreement (OHI2, TOW2, TRO1, VARS and WARK) are labelled, as are the five stations co-located with radiosonde observations (CAS1, CHUR, HOB2, TIXI and UFPR)

The accuracy of a GNSS-estimated tidal displacement is a function of data completeness within each daily PPP session, and the entire data processing window size. Penna et al. (2015) found that if at least 2.5 years of data are used, a harmonic displacement in the semi-diurnal tidal band may be estimated to within about 0.4 mm. They also found that at least 70% data coverage is needed over the given time span. Therefore to be conservative in selecting our data set, we used globally distributed stations which had 90% annual coverage for at least three consecutive years between 2012.0 and 2019.0. Daily station data files were only considered as candidates if there were at least 20 h of GPS and GLONASS continuous data and if the GPS analyses of Blewitt et al. (2018) for the station per day resulted in sub-3 cm values both for the RMS of the daily post-fit residuals from all satellites and for the formal error of the estimated daily 3D coordinates (ftp://data-out.unavco.org/pub/products/unr_qa). For validation of the GLONASS-estimated OTL displacements, GPS-derived OTL displacements using established methodology could be used for most constituents, but for the K2 and K1 constituents which are expected to be problematic for GPS, we must validate using OTL displacements computed by forward modelling. Therefore, after assuring data completeness as described above, we further restricted our choice of GNSS stations to locations where precise and accurate tidal displacement modelling is possible. This must in principle include the modelling of the Earth body tide, but referring to Yuan et al. (2013), we expect sub-millimetre uncertainty for this at any station. Therefore we are concerned only with the accuracy of predicted OTL displacement, which is a function of errors in each of the ocean tide model, the Green’s function incorporating the Earth model, and the computational strategy for convolving these. Penna et al. (2008) found sub-millimetre agreement for the convolution integral computation using different OTL software packages, even in the worst case of coastal stations, and agreement better than 0.2 mm for stations more than 150 km inland. On the other hand, Bos et al. (2015) reported 0.2–0.4 mm disagreement between GPS observations and the predicted M2 OTL height displacement using a Green’s function that accounted for anelasticity effects, commensurate with the effects of the established GPS observation error and uncertainties in the ocean tide models that they used. This suggests that computational and Earth model errors can be reduced to negligible amounts provided a suitable Green’s function and ocean tide model are used. Hence, we believe ocean tide model error remains the main source of potential uncertainty for OTL displacement prediction.

To determine the GNSS station locations at which ocean tide model errors are minimised, OTL displacements based on eight global ocean tide models, FES2014b (Carrère et al. 2016), GOT4.10c (Ray 2013), TPXO8-Atlas (Egbert and Erofeeva 2002), NAO.99b (Matsumoto et al. 2000), HAMTIDE11a (Taguchi et al. 2014), DTU10 (Cheng and Andersen 2011) and EOT11a (Savcenko and Bosch 2012), were computed using NLOADF, and the phasor differences from the mean of the displacements for each cell of a 0.25°×0.25° global grid were generated. The RMS magnitudes of these phasor (vector) differences for the modelled height component for the M2 and K1 constituents are shown in Fig. 1, and similar maps for the smaller constituents N2, K2, O1, P1 and Q1 are provided in “Appendix B”. The largest inter-model discrepancies of about 3 mm arise around the Weddell Sea, the Ross Sea, Baffin Bay, Baffin Island (outside the TOPEX/POSEIDON and Jason altimetry satellite data coverage and where the ice grounding zone is poorly determined) and the Philippines, while there are further more localised areas where the discrepancies are about 1 mm, such as in the Arctic Ocean, northern Australia, the Gulf of Alaska and the north coast of Brazil. Therefore, we only considered GNSS stations away from these areas and selected 49 stations for GNSS processing which all fulfilled the criterion ${\text{Max}}\left\{ {{\text{RMS}}_{i} ,\;i = M_{2} , N_{2} ,K_{2} , K_{1} , O_{1} , P_{1} , Q_{1} } \right\} < 1\;{\text{mm}}$ as well as fulfilling the GNSS data criteria described above. As stated above, we used an elastic Green’s function based on PREM for all of our computations, and at our 49 stations, the M2 constituent height component displacements differ by only 0.16 mm RMS from when the anelastic Green’s function of Wang et al. (2020) is used. Ocean tide model variations caused 0.7–0.8 mm RMS inter-model agreement for the predicted M2 OTL height displacements at OHI2, TOW2, TRO1, VARS and WARK (labelled in Fig. 1), which is mostly caused by 0.8–1.5 mm discrepancies arising from the NAO.99b model. If this model is excluded, the RMS inter-model agreement per station is reduced to 0.4–0.7 mm, but these stations are still the worst-performing. All other stations’ predicted OTL displacements agree better than 0.5 mm regardless of ocean tide model choice.

2.2 PANDA software validation

As we are not aware of any previous publications using PANDA kinematic PPP to estimate OTL displacements, we initially assessed its GPS-only capability using two tests. First, we introduced a synthetic harmonic displacement signal and assessed how well it may be recovered using our PANDA kinematic PPP estimated GPS height time series. Second, the power spectral density (PSD) from GPS kinematic PPP height time series from PANDA were compared with those using the GNSS-Inferred Position System (GIPSY) software, which is regarded as valid because GIPSY height time series have been shown by Penna et al. (2015), Bos et al. (2015) and Martens et al. (2016) to estimate OTL displacement to an accuracy of better than 0.5 mm.

All data from all stations marked on Fig. 1 (and listed in “Appendix A”) were processed using PANDA in GPS-only mode with a 10° elevation angle cut-off, and a 5 mm amplitude (phase assigned as zero at J2000) synthetic harmonic height displacement with 13.96 h period was applied to the data in order to test the harmonic displacement measurement accuracy and precision. This follows the validation methodology of Penna et al. (2015), except here we implemented this by changing the satellite instantaneous position rather than the nominal reference coordinate of the ground station. At each data epoch, we generated a height displacement signal in the GNSS station’s local topocentric frame. Then, the station’s approximate latitude and longitude were used to construct the matrices to convert from the topocentric frame to the geocentric Earth fixed frame of the orbits. After converting the synthetic signal 3D coordinates (with zero values for the east-west and north-south components) to the IGS orbit coordinate frame in this way, the displacements were applied to the satellite positions. Similar to Penna et al. (2015), we then varied the process noise values of the station coordinates and the ZWD, to minimise the synthetic signal recovery error, estimated height repeatability, RMS of the observation post-fit residuals and RMS discrepancy between GPS-estimated tropospheric delay and that estimated from nearby radiosonde data where available. Based on analysis of five of the stations in different parts of the world (CAS1, CHUR, HOB2, TIXI and UFPR, labelled on Fig. 1), we found optimum values of $1.0\;{\text{mm}}/\sqrt {\text{h}}$ and $3.0\;{\text{mm}}/\sqrt {\text{s}}$ for the ZWD and coordinate process noise, respectively, and hence these values are used for the GNSS kinematic PPP data processing throughout the rest of this paper. We evaluated the recovery error in the spectral domain instead of the time domain. In Fig. 2a, the phasor differences between the true synthetic signal and its estimated values at all stations are shown. As this figure indicates, the residual vectors are randomly distributed with a very small mean $R_{\text{m}} = \left( {0.02, 0.01} \right)$  mm. Therefore, we applied the Rayleigh distribution (e.g. Maymon 2018) for the statistical assessment of the synthetic signal recovery error magnitude, and the best fit probability distribution function (PDF) and its equivalent cumulative distribution function (CDF) are shown in Figs. 2b and 2c, respectively. For more than 95% of the tested GNSS stations, the synthetic signal was recovered with an error less than 0.5 mm in magnitude. This is approximately equivalent to the 0.2–0.4 mm RMS reported by Penna et al. (2015) for GIPSY, but uses more stations (49 rather than 21) which are distributed globally, not just in western Europe. The PANDA solution uses float not fixed carrier phase ambiguities. The similarity between the PANDA synthetic signal displacement recoveries and those of Penna et al. (2015) also suggest that for tidal constituents with periods clearly distinct from 12 or 24 h, there is no significant degradation in using 24 h session lengths with concatenated 24 h orbits and clocks rather than 30 h GPS processing session lengths with 30 h orbit and clocks.

Fig. 2

a Signal recovery error phasors from the introduction of a 13.96 h harmonic height displacement for the PANDA GPS-only height solutions, b normalised probability distribution function (PDF) histogram of their vector magnitudes, and c cumulative probability distribution (CDF) over all 49 GNSS stations listed in Appendix A. The smooth red curves in (b) and (c) are for the fitted distribution functions

To compare directly with the PANDA GPS height time series, GPS data over the same time span from all 49 stations were processed using GIPSY v6.4 in kinematic PPP mode, with the processing method following that described in Bos et al. (2015). The key differences between the PANDA and GIPSY processing are that in GIPSY: the VMF1 mapping function was used; the data were processed in 30-hour sessions and the central 24 h of estimated coordinates extracted and concatenated; JPL repro3 fiducial satellite orbits and 30-second clocks computed in the IGb14 reference frame were held fixed; and the nominal FES2014b / PREM Green’s function OTL displacements applied were computed in the CM frame to ensure compatibility with the JPL orbits and clocks. The resulting mean (stacked) PSD plot for the ambiguity-float GIPSY height time series for all 49 stations is shown in Fig. 3, and superimposed on it is that from the PANDA GPS-only processing. It can be seen that they are very similar, with the PANDA results showing slightly (11–18%) more noise PSD averaged across the non-tidal bands 0.2–0.8 cycles per day (cpd), 1.2–1.8 cpd and 2.2–2.8 cpd. This confirms that PANDA GPS-only processing gives commensurate kinematic float PPP results to GIPSY. This similarity exists despite these solutions using different orbit and clock products which have different reference frames (the ESA products are operational and initially in the IGS08 frame but switched to IGb14 at GPS week 1934, whereas the JPL ones are repro3 products in the IGb14 frame) and may be subject to changes in ESA processing strategy over time for the operational products. However, this similarity further substantiates the findings of Penna et al. (2015) who noted that OTL displacement estimates are not sensitive to reference frame changes.

Fig. 3

Mean stacked power spectral density (PSD) for the GPS-derived height time series for the 49 globally distributed GNSS stations processed using both PANDA (ambiguity-float) and GIPSY (ambiguity-float and ambiguity-fixed). All solutions used a 10° elevation cut-off angle. The shaded bandwidths (0.2–0.8 cpd, 1.2–1.8 cpd and 2.2–2.8 cpd) are used for the noise PSD comparisons. The lower panes show enlargements of the diurnal and semi-diurnal frequency bands

Previous studies (e.g. Penna et al. 2015; Bos et al. 2015; Martens et al. 2016) have used ambiguity-fixed GPS kinematic PPP within GIPSY (Bertiger et al. 2010) as the most robust solution for the GPS-derived OTL displacement, so this will be taken as the reference solution for comparison of the ambiguity-float PANDA GPS-only, GLONASS-only and combined GPS+GLONASS solutions later in this paper. To illustrate the effect of ambiguity fixing, Fig. 3 also compares the stacked PSDs of our ambiguity-fixed and ambiguity-float GIPSY GPS solutions. Ambiguity fixing leads to a reduction in noise across the entire frequency range (35–45% smaller noise PSD in the three non-tidal bands mentioned above), although this reduction is marginal at the highest frequencies. We will in the next section investigate to what extent the addition of GLONASS data can mitigate the lack of ambiguity resolution in our PANDA solutions, and constellation-related GPS errors. A notable feature of all solutions shown in Fig. 3 is the frequency comb of increased noise at frequency multiples of K1 (23h56m period) and K2 (11h58m period), arising from errors in GPS which are sidereally repeating (station-satellite geometry and multipath) and orbitally repeating (satellite orbits and clocks), respectively. These errors, resulting from the 11h58m orbital period of GPS satellites, are a principal motivating factor for including GLONASS in our analysis, although they are accompanied by some daybreak noise, centred on frequency multiples of S1 (24 h period), which we would expect to persist in all solutions based on 24-hour data segments.

3 GLONASS data contribution to OTL displacement measurement

The quality of kinematic PPP solutions is very sensitive to the number of satellites and their geometric distribution at each epoch (e.g. Li et al. 2015). The GPS constellation consists of at least 24 satellites distributed in six near-circular orbits of approximate radius 26,559 km, inclined at 55° to the equatorial plane, with a 60° longitude separation between their ascending nodes. The GLONASS constellation also consists of 24 operational satellites, but they are distributed evenly across three near-circular orbits with approximate radius 25,471 km, inclination angle 65° and longitude separation of 120° for the ascending nodes. These differences in satellite constellation change the temporal and spatial variation in GNSS satellites’ availability and viewing geometry, and the consequent dilution of precision (DOP) in different locations; thus, kinematic PPP performance is affected (Pan et al. 2017). In particular, to estimate independent coordinates and receiver clock terms at each epoch within a phase-connected data arc, a minimum of four satellites is required for a single-constellation solution, or five satellites for a dual-constellation solution where the GPS–GLONASS system time offset also needs to be estimated. Epochs when this minimum is approached, or when the geometric dilution of precision is high, may not achieve reliable outlier identification, and hence, the position estimates may be unreliable (especially as 30-minute tropospheric parameters and constant ambiguity parameters are also estimated in our solutions).

Using the initial elevation cut-off angle of 10°, we noted particularly poor performance of some GLONASS-only PPP solutions, which we investigated as follows. We used the TEQC program (Estey and Meertens 1999) to inspect the RINEX observation files of all stations from 00:00 UTC on 15 January to 00:00 UTC on 21 January 2016, a sample time span during which all stations recorded all 30-second data epochs with no receiver tracking outages. The average daily percentage of epochs for which at least seven GPS or seven GLONASS satellites were recorded is shown in Fig. 4. It can be seen that when a 10° mask angle is used, all stations obtain data from at least seven GPS satellites at virtually all epochs, whereas for GLONASS data this success rate varies with station latitude, from around 50% for latitudes between 20° and 30° rising to at least 95% at latitudes of 50° and above. Figure 4 also indicates that for stations with latitudes less than 50°, reducing the mask angle to 5° can significantly increase the percentage of epochs with ample GLONASS observations. Although satellites at lower elevation angles will have lower quality observations because of increased atmospheric propagation errors and multipath, this is mitigated by elevation angle dependent data weighting. In PANDA, for any observation collected at an elevation angle (E) less than 30°, the pre-defined standard error is scaled by $\left\{ {2\sin \left( E \right)} \right\}^{ - 1}$, following Gendt et al. (2003).

Fig. 4

Mean percentage of epochs with at least seven recorded satellites, as a function of station absolute latitude, for six consecutive days in January 2016 for all 49 stations, for different elevation cut-off angles

Hence, we classified stations into two groups based on their latitude: stations within 50° of the equator, and those at higher latitudes, to evaluate the impact of the data mask angle on kinematic PPP performance. After running kinematic PPP solutions for all stations with 5° and 10° elevation cut-off angles for GPS-only as well as GLONASS-only data, the mean PSDs of the estimated height time series for each region were computed. Figure 5 demonstrates slightly lower performance for the GPS-only kinematic PPP solution for stations in the equatorial band, compared with the high-latitude group. For GPS, mean vertical DOP improves slightly at lower latitudes, but we hypothesise that this is offset by greater atmospheric delay variability which impacts the position estimates. Reducing the elevation cut-off angle improves the time series precision very slightly in both regions, which can be explained by the typically increased number of recorded GPS measurements and reduced DOP at each data collection epoch. In Fig. 6, we present the mean stacked PSDs for the estimated height time series from GLONASS-only data, which also show larger noise for the lower-latitude group. However, in this case there is much smaller latitudinal variation in mean DOP, and we hypothesise that the effects of atmospheric delay variability are more extensively compounded because of the smaller number of satellites typically observed. As can be seen in Fig. 6 (middle and bottom panels), the amplitude modulation of the K2 and K1 constituents on the GLONASS satellite ground track repetition signal $\left( {{\text{K}}_{ 1 / 8} } \right)$ causes peaks which are symmetrically distributed around K2 and K1, but which are not present in the equivalent plots for GPS shown in Fig. 5. Figure 6 also indicates that a reduction in data processing elevation cut-off angle enhances GLONASS-only kinematic PPP performance more for lower than for higher latitude stations. Therefore because of this improved precision with a 5° instead of 10° elevation cut-off angle for both GPS and GLONASS constellations and across all latitude bands, we use a 5˚ cut-off angle for all PANDA GPS-only, GLONASS-only and combined GPS+GLONASS data processing for the remainder of this paper.

Fig. 5

Mean stacked PSD of the height time series from PANDA GPS-only kinematic PPP ambiguity-float solutions with different elevation cut-off angles. Stations with absolute latitude (ϕ) greater than 50° are in the left panel, lower-latitude stations are in the centre panel, and the right panel is for the entire data set

Fig. 6

Similar to Fig. 5 but for GLONASS-only data

In Fig. 7, the mean stacked height time series PSDs of the GPS ambiguity-fixed solutions from GIPSY, and those for the GPS, GLONASS and combined GPS+GLONASS ambiguity-float solutions from PANDA are compared. The GLONASS-only solution has generally greater noise than GPS-only, likely because of fewer available GLONASS satellites especially in mid-latitude areas, and lower quality GLONASS satellite clock/orbit products (e.g. Prange et al. 2017). However, by combining GPS and GLONASS data in a float solution, the noise level of the estimated height time series is considerably reduced, and it shows generally similar or even smaller noise compared with the GPS-only ambiguity-fixed solution in GIPSY (0–8% reduction in noise PSD across the 0.2–0.8 cpd, 1.2–1.8 cpd and 2.2–2.8 cpd non-tidal bands). This demonstrates the benefit of incorporating GLONASS data if an ambiguity-fixing algorithm is not implemented in a PPP software package, or when either uncalibrated phase delay (UPD) information or a dense regional network, which are required for PPP ambiguity fixing (e.g. Bertiger et al. 2010: Geng et al. 2011), are not available.

Fig. 7

Mean stacked height time series PSDs from GIPSY ambiguity-fixed GPS-only solutions and GPS-only, GLONASS-only and combined GPS+GLONASS ambiguity-float solutions in PANDA. For all PANDA solutions, a 5° elevation cut-off angle is used. The shaded bandwidths (0.2–0.8 cpd, 1.2–1.8 cpd and 2.2–2.8 cpd) are used for the noise comparison

Although the noise level is generally higher, most of the peaks at frequencies n*K1 in the GLONASS-only PSD are smaller in absolute terms than those in any of the GPS-only solutions. This is because the 11h16m orbital period of GLONASS satellites does not combine with the sidereal rotation of the Earth to create an exact station-satellite geometry repeat as it does for GPS, so sidereally repeating errors such as multipath are randomised and much reduced on average in a GLONASS-only solution. However, small errors remain because there does exist a weak approximate geometry repeat arising from the interaction between the 2.125 GLONASS satellite orbits per sidereal day and the equal separation of eight satellites per GLONASS orbital plane. This means that after one sidereal day the satellite geometry as seen from a station will repeat, although different satellites will be involved. These small peaks can be seen in the GLONASS spectrum, with larger peaks at 3*K1 (K3) and 9*K1 (K9) caused by the 120° longitude separation of the three GLONASS orbital planes (see Daly (1988) for a related discussion of GLONASS viewing geometry repeat). Also, the GLONASS solution shows slightly increased noise at period ${\text{K}}_{ 1 / 8}$ and its frequency multiples, caused by the true GLONASS geometry repeat interval of 8 sidereal days. The combined ambiguity-float GPS+GLONASS PANDA solution is still contaminated by the sidereally repeating errors arising principally from GPS and an overtone of the abovementioned ${\text{K}}_{ 1 / 8}$ artefact, but whereas overall noise levels are similar, the magnitude of all n*K1 peaks is reduced compared with any of the GPS-only solutions.

4 Comparison between GNSS-derived and modelled OTL displacements

We inspect OTL height displacements for the M2, N2, K2, K1, O1, P1 and Q1 constituents obtained from GPS-only, GLONASS-only and combined GPS+GLONASS solutions at each of the 49 stations. The vector differences between the predicted (modelled) and GNSS-derived OTL displacements are shown in Fig. 8, and their statistics summarised in Table 2. The largest M2 residuals (of about 1.2 mm even for the combined GPS+GLONASS solution) are for the stations TOW2, TRO1, VARS and WARK, for which 0.7–0.8 mm inter-model disagreement for the predicted M2 OTL height displacement was noted in Sect. 2.1. Figure 8 demonstrates that the vector differences for all constituents are distributed randomly around zero with a mean well below 0.5 mm, which again leads us to use the Rayleigh distribution for their statistical analysis. The estimated OTL height displacement residuals with their best-fitted Rayleigh CDF are presented in Figs. 9, 10, 11, 12, 13, 14 and 15.

Fig. 8

Vector differences between GNSS-derived and modelled OTL height displacement for all 49 stations for M2, N2, K2, K1, O1, P1 and Q1. In each panel, the mean of all vector differences $\left( {R_{\text{m}} } \right)$ is provided. Note that K2, K1 and P1 are plotted with a different scale to the other constituents. For M2, phasors are highlighted in cyan for stations OHI2, TOW2, TRO1, VARS and WARK, which show larger disagreement among ocean tide models

Table 2 95th percentile of the magnitude of the vector differences between GNSS-derived and modelled OTL height displacement for the 49 stations, for GPS-only, GLONASS-only and GPS+GLONASS solutions

Fig. 9

Magnitude of vector differences between GNSS-derived and modelled M2 OTL height displacement. In the lower panels, the observed cumulative distribution function (CDF) with its fitted counterpart (based on the Rayleigh probability distribution function (PDF)) is shown and the 95th percentile of the fitted CDF is labelled

Fig. 10

Similar to Fig. 9 but for the N2 constituent

Fig. 11

Similar to Fig. 9 but for the K2 constituent (note the different scale matching K1 and P1)

Fig. 12

Similar to Fig. 9 but for the K1 constituent (note the different scale matching K2 and P1)

Fig. 13

Similar to Fig. 9 but for the O1 constituent

Fig. 14

Similar to Fig. 9 but for the P1 constituent (note the different scale matching K2 and K1)

Fig. 15

Similar to Fig. 9 but for the Q1 constituent

As depicted in Fig. 9 for the M2 constituent, the estimated OTL height displacement residuals with GPS-only and GLONASS-only measurements are smaller than 1.2 mm and 1.3 mm, respectively, at about 95% of the processed stations. When excluding the five stations at which the largest M2 disagreements arose (OHI2, TOW2, TRO1, VARS and WARK), these 95% limits are slightly reduced to 1.0 mm and 1.2 mm (not shown in Fig. 9). This indicates the near-similarity in capability of GPS-only and GLONASS-only data to estimate OTL displacement for M2. The slight improvement in vector difference residuals by a factor of 1.2 with GPS rather than GLONASS is also commensurate with the PSD differences around the M2 frequency shown in Fig. 7. Also in accordance with the PSD GPS+GLONASS noise reductions over GPS-only and GLONASS-only, the combined GPS+GLONASS data provide the smallest residuals for M2: for the 44 better-modelled stations the 95th percentile is reduced to 0.9 mm for this estimate and the mean magnitude of these residuals is 0.4 mm, commensurate with the ambiguity-fixed GPS results of Bos et al. (2015). In comparison, Fig. 10 shows that for N2 OTL height displacement, which is only marginally affected by ocean tide model uncertainty, the estimated residual with combined GPS+GLONASS is smaller than 0.3 mm at 95% of the 49 stations, compared with 0.5 mm and 0.6 mm for GPS-only and GLONASS-only, respectively. Similar behaviour for the estimated O1 height residual can be seen in Fig. 13 and, as for N2 and M2, the improvements in the residuals with GPS+GLONASS are commensurate with the PSD reductions over GPS-only and GLONASS-only shown in Fig. 7. We suggest that these results, for constituents whose OTL modelling uncertainty is low, are indicative of the inherent GNSS measurement error budget at frequencies well separated from the sidereal and satellite orbit and geometry repeat periods. Poorer agreement at M2 is at least partly due to the greater OTL modelling uncertainty, but might also indicate systematic lunar-origin errors in satellite orbit and clock or Earth body tide modelling.

Figure 11 clearly demonstrates the problem of measuring OTL displacement at the K2 frequency from GPS data. The 95th percentile of the estimated K2 height residuals estimated by GPS is 4.4 mm, which is much larger than any uncertainty in OTL modelling and more than two times larger than its counterpart estimated by GLONASS (2.0 mm). Hence, the ability of GLONASS to partially overcome GPS problems in measuring K2 tidal displacement is confirmed. However, the lack of GLONASS agreement to within the level of OTL modelling uncertainty that is indicated by the results for N2 and O1 implies that systematic errors remain, which we suggest may be due to overtones of sidereally repeating errors such as multipath arising from the approximate geometry repeat of the GLONASS constellation. Figure 11 also indicates that the increase in satellite availability and better DOP in the combined GPS+GLONASS kinematic PPP can compensate GPS-specific error in the estimated K2 tidal displacement to some extent, but the latter error dominates and so a combined solution (95th percentile residual 2.4 mm) is not as accurate as GLONASS-only. For the K1 tidal constituent shown in Fig. 12, the 95th percentiles of the GLONASS-derived and combined GPS+GLONASS-estimated K1 residuals are 2.8 mm and 2.6 mm, respectively, roughly two-thirds of the GPS-derived value of 3.9 mm. In comparison to K2, these larger discrepancies might imply further systematic errors in addition to the fundamental sidereally repeating geometry related errors. Such errors may arise from the 24-hour data segments used in processing and/or orbit integration, as evidenced by the larger discrepancies also noted for the P1 constituent (Fig. 14) which is similarly close to 24 h in period. The discrepancies at P1 are identical for GPS-only and GLONASS-only solutions (95th percentile 2.4 mm), indicating that they are not related to orbital or geometry repeat period, but reduce to a 95th percentile of 1.3 mm for the combined solution as expected in accordance with the decreased overall noise level.

It was anticipated from the noise reductions shown in Fig. 6 that a more robust kinematic PPP solution would arise for the GLONASS-only solutions at higher latitude stations. Therefore in Figs. 9, 10, 11, 12, 13, 14 and 15, we grouped the residuals into two latitude bands, with smaller GLONASS-only M2, N2, O1, P1 and Q1 residuals seen for the higher latitude band than the lower. To quantify this, we computed the 95th percentiles for the estimated OTL displacement residuals by latitude band and plot these in Fig. 16 for each of the GPS-only, GLONASS-only and GPS+GLONASS solutions. It can be seen that the M2, N2 and O1 OTL displacements can be measured by GLONASS data with similar accuracy to the GPS observations for high-latitude stations, whereas the accuracy of the GLONASS-derived M2, N2, O1, P1 and Q1 estimates is reduced by around 0.2–1.9 mm for the low latitude stations. For K2 and K1, the station latitude effect cannot be seen because the K2 and K1 error sources discussed above dominate.

Fig. 16

95th percentile for the GNSS minus model M2, N2, K2, K1, O1, P1 and Q1 OTL height displacement residuals grouped by station latitude. The cyan bars shown for M2 were computed after excluding the poorly modelled stations mentioned in Fig. 8 caption

5 Discussion and conclusions

We have validated PANDA’s robustness as a kinematic PPP displacement estimation software by comparing the spectral characteristics of its height time series noise to those from GIPSY, at hourly-to-weekly periods. We used a network of globally distributed GNSS stations fulfilling daily and annual data completeness, located in regions with sub-millimetre consistency in predicted OTL height displacement when computed with different ocean tide models. We found that a low (5° instead of 10°) elevation cut-off angle mask was especially beneficial for processing lower-latitude GLONASS-only solutions and had small but positive impact in other situations. Our investigation of GPS-only, GLONASS-only and combined GPS+GLONASS observation processing in kinematic PPP mode demonstrates three main benefits of incorporating GLONASS data, with particular relevance for measuring OTL displacement.

First, combined GPS+GLONASS kinematic PPP with float ambiguity estimation is as precise as GPS-only fixed ambiguity PPP. However, the former is more flexible, as it is independent of UPD corrections. Even with available UPD information, the ambiguity-fixing success rate will decrease when estimated float ambiguities are not precise enough (Teunissen 2017), for instance in the situation of poor DOP, extreme ionospheric activity or short phase-connected arcs.

Second, in addition to noise reduction in the combined GPS+GLONASS kinematic PPP compared with single-constellation solutions, it is verified that GLONASS-only float solutions are able to measure the M2, N2, O1 and Q1 constituents of OTL height displacement with almost similar accuracy to GPS-only float solutions. With GLONASS-only, 95% of tidal displacements agreed with forward geophysical models to within 0.6–1.3 mm for the M2, N2, O1 and Q1 constituents, compared with 0.4–1.2 mm for GPS-only. Hence, GLONASS-derived M2, N2, O1 and Q1 OTL displacements can be used as a check for GPS-derived ones and vice versa. Furthermore, OTL displacement estimation from a float solution using combined GPS+GLONASS observations can be as robust as a GPS-only ambiguity-fixed solution for these constituents.

Third, we have demonstrated the improved ability of GLONASS data to resolve OTL height displacements at the luni-solar semi-diurnal and diurnal periods (K2 and K1) which are not reliably measurable by GPS on account of the latter’s orbital period and station-satellite geometry repeat interval. We found very distinct improvement for purely GLONASS-derived K2 tidal displacement compared to its GPS-derived counterpart (2.0 mm rather than 4.4 mm 95th percentile residual values), and also improved compared with the combined GPS+GLONASS estimate (2.4 mm 95th percentile) which appears to be dominated by GPS-related errors. For the K1 constituent, the GLONASS-only and combined solutions are of comparable quality. However, even the best solutions at K1 and K2 do not agree with the modelled OTL at the level achieved for M2, N2, O1 and Q1. This disagreement may be caused by sidereally repeating errors which also exist for GLONASS because of its approximate sidereal station-satellite geometry repeat, and errors arising from the use of 24-hour data segments which also affect the nearby P1 constituent.

When several years of complete constellation data together with corresponding high accuracy satellite orbits and high-rate clocks are available for the Galileo and BeiDou systems (which have further differences in orbital and geometry repeat periods), a combined GPS+GLONASS+Galileo+BeiDou kinematic PPP solution may achieve a further reduction in the K1 and K2 residuals. For example, Abraha et al. (2018) showed that even limited Galileo data when added to GPS+GLONASS data reduced the amplitude of spurious propagated tidal signals in GPS coordinate time series. Hence, a full, four-constellation kinematic PPP solution using longer data segments could in future provide the potential to be utilised for the refinement of solid-Earth Green’s functions and numerical ocean tide models, including for the K1, K2 and P1 constituents.

Appendix

1.1 A: GNSS station and observation information

See Fig. 17, Table 3.

Fig. 17

Processed GNSS data time span with percentage data availability per day at each station

Table 3 GNSS station name, coordinates and receiver operator or data provider

1.2 B: OTL height displacement inter-model agreements

See Fig. 18, 19.

Fig. 18

RMS agreement of the magnitudes of the vector differences for the predicted OTL height displacement (in mm) per cell of a 0.25° global grid for the semi-diurnal N2 (top) and K2 (bottom) constituents based on seven recent ocean tide models (FES2014b, GOT4.10c, TPXO8-Atlas, NAO.99b, HAMTIDE11a, DTU10 and EOT11a), with GNSS stations used in this paper shown as dark blue dots

Fig. 19

RMS agreement of the magnitudes of the vector differences for the predicted OTL height displacement (in mm) per cell of a 0.25° global grid for the diurnal O1 (top), P1 (middle) and Q1 (bottom) constituents based on seven recent ocean tide models (FES2014b, GOT4.10c, TPXO8-Atlas, NAO.99b, HAMTIDE11a, DTU10 and EOT11a), with GNSS stations used in this paper shown as dark blue dots. The colour scale saturates at 2 mm RMS (maximum RMS for O1 is 3 mm)