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

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.


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. 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 (Baker 1984).

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 hours, 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 semidiurnal 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  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 seconds, ZWD every 30 minutes and its horizontal gradients every 60 minutes, 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.

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 globallydistributed 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 Figure 1 and listed in Appendix A along with the data availability and spans used) which fulfilled the criteria now described.
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 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. Therefore we only considered GNSS stations away from these areas, and selected 49 stations for GNSS processing which all fulfilled the criterion { , = 2 , 2 , 2 , 1 , 1 , 1 , 1 } < 1 as well as fulfilling the GNSS data criteria described above. As stated above, we used an elastic inter-model agreement for the predicted M2 OTL height displacements at OHI2, TOW2, TRO1, VARS, and WARK (labelled in Figure 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.

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 Figure 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 hour 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 Figure 1), we found optimum values of 1.0 /√ℎ and 3.0 /√ 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 Figure 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 = (0.02, 0.01) 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 Figures 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 hours, there is no significant degradation in using 24 hour session lengths with concatenated 24 hour orbits and clocks rather than 30 hour GPS processing session lengths with 30 hour orbit and clocks.
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 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 Figure 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 (24h period), which we would expect to persist in all solutions based on 24-hour data segments.

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 26559 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 25471 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 Figure 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°-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 {2 ( )} −1 , following Gendt et al. (2003).
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 GPSonly 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 Figure 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 Figure 6 (middle and bottom panels), the amplitude modulation of the K2 and K1 constituents on the GLONASS satellite ground track repetition signal (K 1/8 ) causes peaks which are symmetrically distributed around K2 and K1, but which are not present in the equivalent plots for GPS shown in Figure 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.    Although the noise level is generally higher, most of the peaks at frequencies n*K1 in the GLONASSonly 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⅛ 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 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 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.

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 Figure 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 section 2. 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 Figures 9-15.
As depicted in Figure 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 Figure 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 Figure 7. Also in accordance with the PSD GPS+GLONASS noise reductions over GPS-only and GLONASS-only, the combined GPS+GLONASS data provides the smallest residuals for M2: for the 44 bettermodelled 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, Figure 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 Figure 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 Figure 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 Figure 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 ( Figure 14) which is similarly close to 24 hours 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 Figure 6 that a more robust kinematic PPP solution would arise for the GLONASS-only solutions at higher latitude stations. Therefore in Figures   9 to 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 Figure 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 GLONASSderived 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.

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 GLONASSonly 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. 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.