Introduction

The earth rotation parameters (ERPs) belong to the fundamental transformation parameters linking the International Celestial Reference Frame (ICRF) and the International Terrestrial Reference Frame (ITRF) (Petit and Luzum 2010). The ERPs are conventionally composed of the X and Y pole coordinates (polar motion, PM), their rates, UT1-UTC, and its rate called length-of-day (LoD). Since the satellite orbits realize an inertial system, processing of data collected by the globally distributed network of terrestrial stations allows for the estimation of ERPs.

The Global Navigational Satellite System (GNSS) is one of the major satellite techniques which delivers reliable information about ERPs (Dow et al. 2009; Ferland and Piraszewski 2009; Mireault et al. 1999; Steigenberger et al. 2006). Altamimi et al. (2016) indicate that the GNSS-based PM series dominate the combined PM within ITRF, which is due to the specific characteristics of the GNSS technique. Essential advantages of GNSS in the context of PM estimation are (1) continuous signal registration, (2) a globally distributed network of ground sites, and (3) a wide and stable constellation of satellites (Ray et al. 2017). This concerns, in particular, the whole multi-GNSS constellation with almost 120 satellites in space, including GPS, GLONASS, Galileo, BeiDou, and QZSS at different altitudes and different revolution periods (Montenbruck et al. 2017; Weiss et al. 2017). In contrast to PM, the GNSS technique is not capable of delivering independent information about UT1 as it requires the ability to separate the axial rotation of the earth relative to the celestial frame and rotation of the satellite orbit nodes in the inertial frame (Senior et al. 2010; Thaller et al. 2007). Nonetheless, the LoD, which is a negative time derivative of UT1-UTC expressed over 24 h, is typically estimated based on satellite techniques (Bizouard et al. 2018). The LoD estimation is sensitive to any uncertainties in orbit modeling, especially those related to the Keplerian orbit rotation parameters (Rothacher et al. 1999).

The ERPs, as well as all the other GNSS-based products, are vulnerable to the spurious effects which originate from orbit modeling issues and GNSS constellation characteristics. Direct solar radiation pressure (SRP) is the main non-gravitational source of accelerations acting on GNSS satellites (Bury et al. 2019; Prange et al. 2019; Rodriguez-Solano et al. 2012). Nowadays, direct SRP is one of the most challenging aspects of the multi-GNSS orbit determination. Arnold et al. (2015) tested various cases of SRP modeling using different sets of empirical parameters of the Empirical CODE Orbit Model (ECOM, Beutler et al. 1994). The latter decomposes the accelerations acting on satellites in the DYB frame, i.e., D—pointing toward the sun, Y—along the solar panel rotation axis, and B—perpendicular to D and Y axes, completing the right-handed orthogonal frame. The authors indicated that some artificial signals at the harmonics of a draconitic year, which is the repeatability period of a particular GNSS constellation with respect to the sun, occur depending on the set of estimated ECOM parameters. A similar conclusion was drawn by Rodriguez-Solano et al. (2014), who indicated that using the adjustable box-wing model for SRP modeling significantly reduces the GPS draconitic signals in the PM rates. Lutz et al. (2016) pointed out that the arc length also has a perceptible impact on the ERPs, especially PM rates. The rates based on the 3-day arc solutions are generally superior to those based on the 1-day arc solutions.

Meindl (2011) showed the differences in ERPs derived from the independent GPS-only and GLONASS-only solutions for 2009–2011. He found that both solutions are fraught with different systematic offsets with respect to the International Earth Rotation and Reference Systems Service C04 (IERS-C04) series. The scatter of the PM estimates was also larger by a factor of two for GLONASS compared to GPS. Motivated by Meindl et al. (2011, 2013), Scaramuzza et al. (2018) studied the impact of the number of GNSS orbital planes on the geodynamic parameters, including ERPs. The study shows that 3 instead of 6 orbital planes in the constellation may lead to spurious signals in PM, especially at the harmonic of 3 cycles per draconitic year (cpy). A combination of different GNSS, regardless of the orientation of their nodes, reduces the differences between system-specific parameters. Unfortunately, the test study did not cover the analyses of Galileo, which, similarly to GLONASS, comprises 3 orbital planes.

Since December 2018, the Galileo system constellation is composed of 24 usable satellites and is nearly fully operational together with the legacy GPS and GLONASS systems. Moreover, the detailed metadata containing optical and geometrical properties for the In Orbit Validation (IOV) and Fully Operational Capability (FOC) satellites have been released since 2017 by the European GNSS Service Centre. The metadata brought new perspectives for advanced SRP modeling based on hybrid or semi-empirical SRP models (Bury et al. 2019; Li et al. 2019). Finally, Galileo is certainly a perfect system to check its suitability for ERP estimation and to confirm the hitherto results based on GPS and GLONASS systems.

Motivation

Since 2012, the first multi-GNSS ERP products are supplied to the GNSS community, based on the operationally delivered MGEX products (Multi-GNSS Pilot Project of International GNSS Service, Montenbruck et al. 2015, 2017). Six MGEX Analysis Centers (ACs) have been routinely delivering daily ERP estimates based on multi-GNSS precise orbit determination (POD) processing: Center for Orbit Determination in Europe (CODE), Wuhan University (WUM), GFZ German Research Centre for Geosciences (GFZ), Japan Aerospace Exploration Agency (JAXA), French Space Agency (CNES/CLS), and the European Space Agency (ESA), generating similar products available only on ESA servers (Deng et al. 2016; Guo et al. 2016; Guo et al. 2017; Prange et al. 2017; Katsigianni et al. 2019; Enderle et al. 2019). Figure 1 shows the comparison of the time series of PM and LoD, which are delivered by the ACs, in reference to the IERS-C04-14 values (Bizouard et al. 2018). One can notice an improvement in ERP precision, which has been achieved since the beginning of the MGEX Pilot Project, especially for the PM.

Fig. 1
figure 1

Differences of ERPs delivered by the ACs in reference to the IERS-C04-14 series

Full size image

The POD strategies of the ACs have undergone several changes over the last few years, because of the experimental nature of the MGEX products, including updates in background models, changes in the approach to SRP modeling, as well as the growing number of satellites and tracking stations. CODE has been providing the longest and continuous series of multi-GNSS solutions since 2012 (Fig. 2). Before 2015, when the reduced ECOM model was applied for the SRP modeling (Springer 1999), the variability of the pole coordinate residuals was at the level of 60, 80 µas, for the X and Y pole coordinates, respectively. The dominant signal was also visible at 3 cpy for both pole coordinates. The 3 cpy signal was explicitly related to inefficiencies in the GLONASS orbit modeling (Arnold et al. 2015; Dach et al. 2019; Lutz et al. 2016). After the change of the SRP modeling strategy for CODE MGEX products in 2015 from the reduced ECOM to the so-called extended ECOM2 (Arnold et al. 2015; Sośnica et al. 2018), the standard deviation of the pole coordinates decreased to the level of 50 µas, and even 40 µas when the number of the Galileo satellites exceeded 17 at the beginning of 2017. Since December 2018, 24 Galileo satellites are included in the solutions of all ACs. Furthermore, since January 2019, ESA and WUM supplement their products with BeiDou-3 satellites. The internal consistency between the ACs, as well as the ERP quality as measured by the offset and standard deviation with respect to the IERS-C04-14 is visibly improving over time. The appearance of new Galileo and BeiDou satellites in the solutions improved the variety of the observation geometries and, therefore, decorrelated the pole coordinates from other estimated parameters, which was also supported by the enhancement of the background models used for data processing. In line with theory, ERPs should benefit from the diversity of satellites with different orbital characteristics as well as with the increase in independent data (Freymueller 2017).

Fig. 2
figure 2

Number of satellites used in the multi-GNSS POD processing by the particular ACs. The gaps in data reflect the unavailability of the products

Full size image

The comparison of the ERP products delivered by the different ACs shows that the quality of ERP may depend on the adopted processing strategy. On the other hand, based on the operational MGEX products only, we cannot state clearly what is the impact of each constituent separately. Therefore, the first objective is to study how the particular satellite systems, including GPS, GLONASS, and Galileo, influence the estimates in the combined multi-GNSS solution. Second, we verify the differences between the ERP estimates obtained from GPS and GLONASS, and Galileo constellations. For that purpose, 3 years of GNSS observations from a globally distributed network of stations have been processed in different test cases. The next ITRF realization, as well as the International GNSS Service (IGS, Johnston et al. 2017) reprocessing, repro3, will include the Galileo system for the first time. Thus, the potential impact of the Galileo system on future IGS products seems to be crucial in terms of processing strategies in IGS.

Methodology

The processing strategy was adopted from the CODE MGEX solutions (Prange et al. 2017). The processing details and the overview of the background models are summarized in Table 1. The computations were performed using a development version of the Bernese GNSS Software (Dach et al. 2015). Observations were collected from the network of approximately 100 stations, all tracking GPS, GLONASS, and Galileo satellites. Thus, we do not expect any differences between GPS, GLONASS, and Galileo that could originate from different ground network distribution (Zajdel et al. 2019).

Table 1 Description of the processing strategy
Full size table

The X, Y pole coordinates and UT1-UTC are estimated with 24 h temporal resolution. The ERP values and their rates are estimated at noons in 1-day solutions, to be in line with the IGS-related products. The ERP representation consisting of the offset and rate can be transformed into two offsets at the beginning and the end of a day by considering in addition sub-daily ERP corrections. Satellite techniques are not capable of determining the daily rotation of the earth over a longer period because of the correlation with the rates of orbital parameters describing the orientation of orbital planes, which is the ascending node and the argument of latitude (Rothacher et al. 2001). However, the singularity is removed by fixing one of the UT1-UTC parameters, (e.g., the offset at the beginning of the orbital arc in the case of the presented processing strategy), to the corresponding value from the reference IERS-C04-14 time series and freely estimating the rate (or the offset at the end of the day). This is a valid approach that is implemented in the Bernese GNSS Software. The reference series of UT1-UTC in IERS-C04-14 is naturally the one obtained from the Very Long Baseline Interferometry (VLBI) series (Nothnagel et al. 2017).

We estimate ERPs based on both 1- and ‘overlapped’ 3-day orbital arc solutions. The 3-day arc solutions are based on the stacked 1-day normal equation systems (NEQs) of three consecutive days separated by 1-day intervals, which means that the middle day of one solution and the first day of the subsequent solution overlap. In 1-day arc solutions, the ERPs are parametrized by the offset and rate for each day independently, as in line with the current IGS requirements. In 3-day arc solutions, the ERPs are parametrized as piece-wise linear with common nodal points at midnights; thus, the 3-day arc solutions are represented by four midnight offsets for each ERP component (Thaller et al. 2007; Lutz et al. 2016). One common set of station coordinates and satellite orbits is estimated for 3 days, while the ERPs are estimated with 24 h temporal resolution. Such an approach ensures that the ERPs are continuous at the boundaries of the middle day, but allows for instantaneous changes in the rates (Lutz et al. 2016). For the analyses, we extract only the set of ERPs which refer to the middle day of the 3-day-long arc. Two offsets, at the beginning and the end of the middle day of the arc, can be transformed back to the offset-and-rate representation at the middle, 12-h epoch, for comparison with rates derived from 1-day solutions or the IERS-C04 series. The estimated orbital arcs are continuous over the established 1- or 3-day windows; however, the instantaneous velocity changes represented as pseudo-stochastic pulses are allowed at all noon and midnight epochs with particular constraints (Table 1).

In total, 6 solutions were designed to assess the individual impact of the GPS, GLONASS, and Galileo constellation on the ERP estimates (Table 2). The computations of system-specific ERPs, i.e., estimating ERPs separately for GPS, GLONASS, and Galileo, are based on the same NEQs as the standard combined solution. It is worth noting that the system-specific solutions share some parameters such as troposphere and station coordinates but leave geocenter coordinates and ERPs as system-specific. The advantages of such a parametrization are described in detail by Scaramuzza et al. (2018). By using the same NEQs as a basis for each solution, we ensure the consistency in the terrestrial reference frame realization for all the solutions. Furthermore, the resulting ERPs and geocenter coordinates directly reflect the system-specific signals, which arise from the changes in observational geometry in different GNSS constellations. The ERPs and geocenter estimates may also include some of the modeling deficiencies, which would be transferred into other parameters, e.g., station coordinates, in the case of independent single-GNSS processing (Scaramuzza et al. 2018).

Table 2 Description of the solutions
Full size table

We compare three approaches of SRP modeling in the Galileo orbit solution. The first approach assumes the empirical ECOM2 model with 7 parameters (3 constant accelerations in DYB directions, sine and cosine once-per-revolution terms in B, and sine and cosine twice-per-revolution terms in D). The second and third approach combines a physically interpreted box-wing and empirical ECOM model with a reduced set of parameters (Bury et al. 2019). The composed box-wing model considers the direct SRP, earth albedo, and infrared radiation, based on the Galileo metadata released in the years 2017–2019 by the European GNSS Service Centre. The hybrid strategy of combining both specific terms of empirical and physical models is beneficial for the orbit quality (Bury et al. 2020; Li et al. 2019). The selection of the set of empirical parameters, which should be chosen in the hybrid solution, is extensively studied by Bury et al. (2019). The authors indicated that in terms of the orbit quality, two cases are satisfactory: (1) the box-wing model and estimated constant accelerations in the DYB directions, (2) the box-wing model and estimated periodic accelerations in the B direction together with constant DYB accelerations (equivalent to reduced ECOM). Yan et al. (2019) stated that the case (2) is also more suited for the estimation of phase center offsets and variations for BeiDou-3 satellites. However, the impact of SRP modeling on ERPs in Galileo processing has not been tackled. Motivated by Bury et al. (2019), we composed both (1) and (2) Galileo POD strategies (Table 2). After the investigation of both hybrid SRP strategies, we concluded that the case (2) is more advantageous for ERP estimation, especially in the case of the 3-day arcs. Figure 3 shows the spectral analysis of the pole coordinate estimates for all three Galileo-based solutions. The pronounced peaks at the periods of 3.4 and 121 days (3 cpy) are significantly larger for the GBB solution than for the GAB solution. Thus, only the GAB solution is considered in the next sections of this study.

Fig. 3
figure 3

Amplitude spectra of the Galileo-based X and Y pole coordinates with respect to IERS-C04-14 for the a 1-day arc solutions and b 3-day arc solutions. The vertical gray lines in this and all the following figures denote the harmonics of a draconitic year

Full size image

Results

In this section, the analyses of GNSS-based ERPs are described. The estimates of the X and Y pole coordinates, their rates, and LoD are compared to the reference values of IERS-C04-14 (Bizouard et al. 2018). In the next sections, the statistical values are provided both for the whole analyzed period and the year 2019 only, because in 2019, 24 Galileo satellites were in operation. The IERS-C04-14 series is provided through the IERS ftp servers in two variants, which are sampled at noon or midnight epochs. The IERS series is a multi-technique solution incorporating VLBI, GNSS,  Satellite Laser Ranging (SLR), and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) data. The combined solution benefits from the merits of all the geodetic techniques, while most of their weaknesses are mitigated. We have to be also aware that IERS-C04 series is smoothed, loosing most of the real PM signal, especially at the 1–2 day window and affecting signals with periods up to 3 days. Bizouard et al. (2018) assessed the general consistency of the GNSS-based IGS results within the IERS-C04-14 combination for the period 2010–2015, which reached an impressive level of 30 µas and 10 µs/day for PM and LoD, respectively. A similar comparison can be made for different periods using the dedicated online tool provided by the Earth Orientation Parameter Product Centre of the IERS. Looking at the relevant period, which corresponds to this analysis, the offset and standard deviation of residuals between IGS-final ERP product and IERS-C04-14 equal 40 ± 60 µas, −18 ± 45 µas, and 1 ± 6 µs/day for the X, Y pole coordinates, and LoD, respectively. In the last part of this section, we discuss ERP misclosures as the internal quality indicator.

Parameter errors of earth rotation parameters

Table 3 summarizes the formal errors of the ERP estimates for particular solutions. The formal errors of the parameters are lower by 30–40% for GPS solutions than for solutions based on GLONASS and Galileo constellations. Such a discrepancy may arise from the different number of satellites in the constellations, which is about 30–32, 21–24, and 14–24 satellites for GPS, GLONASS, and Galileo, respectively. However, in the case of the cumulated multi-GNSS constellation in the GRE solution, the formal errors roughly correspond to the GPS solution. Thus, the formal errors also reflect the mutual inconsistencies of the individual subsystems in the combined solution. While using 3-day arcs, we theoretically expect an improvement of the parameter errors at the level of \(\sqrt 3\) ≈ 1.73, when compared to the 1-day arc solutions because approximately 3 times more observations are used. However, the observed improvement is even higher. The extension of the arc length from 1 to 3 days decreases the formal errors by almost a factor of 3 for GLO, GAL, and GAB solutions, and by a factor of 2.2 for GRE and GPS solutions. Thus, the decrease in formal errors should also be attributed to the reduction of correlations between parameters in the 3-day solutions, because of the continuity of orbits and ERP over the orbital arc (Lutz et al. 2016). We discuss the raw formal error outputs as given by Bernese GNSS Software, without using any scaling factors neither for 1- nor 3-day arc solutions.

Table 3 Mean formal solution errors for the respective solutions
Full size table

The parameter errors from 1-day arc solutions vary in time significantly. Figure 4 shows the time series and the spectral analysis of the parameter errors from the particular solutions delivered in the 1-day arc and 3-day arc cases. The characteristic patterns are reduced almost to zero for the corresponding 3-day solutions. These variations most likely originate from the deficiencies in SRP modeling, since the errors correlate with the sun elevation above orbital planes (Fig. 4a). A similar dependency has already been indicated for the Z component of the geocenter coordinates (Scaramuzza et al. 2018; Zajdel et al. 2019). Interestingly, Bury et al. (2019) pointed out that in the case of the Z component of geocenter coordinates, when the hybrid approach to SRP modeling is applied (equivalent to GAB solution), the variations of the parameter formal errors are significantly reduced. In the case of the ERPs, the characteristic patterns are still visible when cross-referencing GAL and GAB solutions (Fig. 4b). Although Galileo and GLONASS comprise three nominal orbital planes, the periodic variations of the parameter formal errors have different patterns. In the case of LoD, the most pronounced signal is visible at the frequencies of 4 cpy and 2 cpy for Galileo and GLONASS, respectively (Fig. 4b). In the case of the Y pole coordinate, the main signal occurs at 4 and 2 cpy for GLONASS and 6 and 4 cpy for Galileo. For the X pole component, the main signal is visible at even harmonics for both GLO and Galileo-based solutions. The combination of different systems almost completely mitigates the system-specific periodic signals for all of the considered parameters (Fig. 4).

Fig. 4
figure 4

Formal errors of ERPs from the respective 1-day arc solutions (colored lines) and corresponding 3-day arc solutions (black lines). a Time series of LoD formal errors; gray lines denote the sun elevation angles above orbital planes (β). All values in µs/days. b Amplitude spectra of the formal errors of ERPs. The vertical gray lines denote the harmonics of a draconitic year

Full size image

Polar motion

Table 4 shows the mean offset and standard deviation of the PM residuals with respect to the IERS-C04-14. We distinguished the statistics for the whole analyzed period as well as the year 2019 only, when 24 Galileo satellites were available. The mean offsets for the GRE and GPS solutions are similar and equal about 50 and − 30 µas for the X and Y pole coordinates, respectively. Moreover, the offsets for the GPS-based PM estimates are consistent with the expected values as derived from the IGS-final products which are dominated by GPS. Interestingly, the mean offset significantly differs between the constellations. The offsets for 1-day arc solutions equal 49, 107, 31, 23 µas for the X coordinate and − 30, −95, − 12, − 15 µas for the Y coordinate, for GPS, GLO, GAL, and GAB solutions, respectively. The offsets for the corresponding solutions do not change significantly when a 3-day arc case is applied (Table 4 and Fig. 5). However, the statistics from the year 2019 show that the offset changes with the development of respective satellite constellations. In the case of the GLONASS constellation, some of the old GLONASS-M satellites have been exempted from service in 2018 and 2019, i.e., SVN 733 (Q2 2018), 716 (Q4 2018), 717 and 734 (Q3 2019). Dach et al. (2019) revealed conspicuous signatures for certain GLONASS satellites when using nominal horizontal satellite antenna offsets used by IGS. In exchange, new GLONASS-M and M + satellites took their spots including 856, 857 (Q3 2018), and 858 (Q3 2019). For Galileo, it turns out that along with the improvement of the constellation, the PM results are getting closer to these from GPS (Table 4).

Fig. 5
figure 5

Time series of PM differences with respect to IERS-C04-14 in the a 1-day arc solutions and b the 3-day arc solutions. The time series are smoothed using a moving average filter with a 5-day window

Full size image
Table 4 Statistics of the polar motion differences with respect to IERS-C04-14. STD denotes a standard deviation
Full size table

The remaining system-specific offsets in PM may originate from modeling issues, such as the lack of precise antenna phase center offsets and variations for Galileo frequencies applied in the processing. Moreover, Zajdel et al. (2019) indicated that geocenter coordinates, as seen by GPS, GLONASS, and Galileo, are also systematically biased. Thus, the realization of the reference frame may be noticeably translated, as seen by different GNSS constellations. As stressed by Ray et al. (2017), the geocenter translations can produce opportunities for systematic rotational errors, especially for no optimal distribution of GNSS stations. The applied processing strategy, which assumes that the ERPs and geocenter coordinates are estimated for each system separately, while the remaining parameters are common to all solutions, may also somehow systematically bias the results. The GPS and GRE results are least scattered having a standard deviation of approximately 40–50 µas. In contrast, the standard deviation equals 91, 66, 66 µas for the X coordinate and 63, 57, 52 µas for the Y coordinate, for GLO, GAL, and GAB solutions, respectively. In the case of the 1-day arc GLONASS- and Galileo-based series, the standard deviation of residuals is even lower when only the year 2019 is considered. The difference in the quality of both PM components suggests a possible relation with the distribution of GNSS ground stations between the eastern and western hemispheres as the X and Y values refer to the earth-fixed frame. In the case of the Y pole coordinate, the 3-day arc leads to a decrease in the standard deviation of residuals from approximately 10% for GAB solutions up to 20% for GPS solution. However, when concerning the X pole coordinate, the 3-day arc makes the standard deviation of the residuals slightly worse up to 23% for the GAB solution.

Figure 6 illustrates the spectral analysis of the X and Y pole coordinate differences with respect to IERS-C04-14. The X pole coordinate is more affected by spurious signals than the Y pole coordinate. The strong 3 cpy signal with the amplitude up to 41 µas is evident in the time series of the X pole coordinate for the GRE, GAL, GAB, and GLO solutions. Interestingly, the 3 cpy signals inherently affect all 3-plane constellations. Moreover, the combination of different GNSS is insufficient to decorrelate this issue entirely. Thus, the GRE solution is also affected by the spurious 3 cpy signal. The 3 cpy signal in the X pole coordinate series even increases in the 3-day arc case. Then, the 3 cpy signal reaches the amplitude of 78, 30, 24, 21 µas for GLO, GRE, GAL, GAB, respectively. Moreover, the signal with a period close to 52 days (7 cpy) and the amplitudes at the level of 20–30 µas appears for all the 3-plane constellations, especially for 3-day arc solutions. In the 3-day GLO and GAL/GAB solutions, we also see some increased spurious signals for the periods in the range of 20–30 days. The orbit modeling deficiencies might cause all of these issues, i.e., SRP mismodeling for GLONASS and Galileo (Bury et al. 2020; Prange et al. 2017). If some significant orbit mismodeling issues exist, their effect might be amplified with long arcs because the same number of parameters must absorb errors accumulated over a longer period. In the case of GLONASS, Dach et al. (2019) described the issues in orbit modeling for the particular GLONASS satellites.

Fig. 6
figure 6

Amplitude spectra of the estimated X and Y pole coordinates with respect to the IERS-C04-14 series for the a 1-day arc solutions and b the 3-day arc solutions. The vertical gray lines denote the harmonics of a draconitic year

Full size image

Polar motion rates

The PM rates \(\dot{X}\) and \(\dot{Y}\) are not explicitly provided in the official IERS-C04 series. For the comparisons, we calculated the reference values of PM rates as the differences of the IERS-C04-14 pole coordinate values of two subsequent midnight epochs. Lutz et al. (2016) concluded that the extension of the arc length from 1 to 3 days has an implicit impact on the PM rates. The rates based on 3-day solutions should be, in general, much better than those based on the 1-day standard IGS solution. Table 5 summarizes the mean offset and standard deviation of the estimated pole coordinate rate residuals with respect to the IERS-C04-14.

Table 5 Statistics of the polar motion rate differences with respect to IERS-C04-14. STD denotes standard deviation
Full size table

Considering the 1-day arc solutions, the consistency of the PM rates with the reference values is of very poor quality as in line with expectations. The change of the orbital arc length from 1 to 3 days results in the improvement of the consistency with IERS-C04-14. In the case of GLO and GAL solutions, the standard deviation is even 2.5 and 3.3 times better for \(\dot{X}\) and \(\dot{Y}\), respectively, when using 3-day arcs compared to 1-day arcs. Finally, based on the 3-day arc solutions, the standard deviation for \(\dot{X}\)(\(\dot{Y}\)) pole coordinates rates are at the level of 98(94), 99(94), 113(107), 128(118), 131(119) µas/day for GRE, GPS, GLO, GAL, GAB solutions, respectively. In the case of GLONASS and Galileo, the standard deviation of residuals is even lower when only the year 2019 is considered. Finally, the mean offset decreases to the level of a several µas/d when switched from 1- to 3-day arc for all the analyzed solutions.

Figure 7 illustrates the spectral analysis of the PM rates. The artificial signals at the harmonics of the draconitic year are visible for all the system-specific solutions. Besides, one can also see some evident artificial peaks in the GLO, GAL and GAB series at the specific periods up to 20 days. The estimation of ERPs using GNSS is limited since the variations in earth rotation are observed by the dynamic system of satellites, which also rotates in conjunction with the earth. Therefore, we may expect spurious signals at the periods which stem from the combination of both the period of earth rotation (sidereal day) and the revolution period of the constellations (11 h 58 min for GPS, 11 h 16 min for GLONASS, and 14 h 05 min for Galileo). This dependency shall be described as follows:

$$\left| {\frac{1}{{n*f_{\text{S}} + m*f_{\text{E}} }}} \right|, \;{\text{with}}\; n,m = \left\{ { \ldots , - 4, - 3, - 2, - 1,0,1,2,3,4, \ldots } \right\}$$
(1)

where \(f_{\text{E}}\) and \(f_{\text{S}}\) are the frequency of earth rotation and satellite revolution, respectively. Then, we can calculate system-specific signals for the particular n and m values. As expected, some artificial signals are clearly visible at periods close to 2.5 days (n = 2, m = −3), 3.4 days (n = 1, m = −2), 10 days (n = 3, m = −5) for Galileo-based solutions and 2.6 days (n = 3, m = −6), 3.9 days (n = 2, m = −4), 7.9 days (n = 1, m = −2) for GLONASS-based solutions. These signals were also visible in the series of pole coordinates (Figs. 3 and 6) and formal errors of ERP estimates (Fig. 4b) but with a noticeably lower magnitude. In the case of Galileo, these kinds of signals are visible regardless of the SRP approach used. When cross-referencing the corresponding 1- and 3-day arc solutions, we observe that most of the spurious signals are almost completely mitigated. The minor exception is the 3 cpy peak with the amplitude of 24 µas in the series of GLO \(\dot{Y}\). In a circular polar motion of angular velocity ω, we may say that \(\dot{y} = \omega x\). Thus, the noticeable 3 cpy signal in the \(\dot{Y}\) series is a consequence of the spurious 3 cpy signal, which prevails in the corresponding series of the X pole coordinate.

Fig. 7
figure 7

Amplitude spectra of the polar motion rate residuals with respect to IERS-C04-14 for a the 1-day arc solutions and b the 3-day arc solutions. The vertical gray lines mark the harmonics of a draconitic year

Full size image

Length-of-day

According to the adopted parametrization of the ERP estimation, we can provide information about the UT1-UTC drift, i.e., LoD. Because the UT1-UTC values in IERS-C04-14 are highly constrained to the VLBI-based values, the analysis of the LoD series can indicate how rapidly a GNSS-derived UT1-UTC would drift away from the reference values. Figure 8 shows the drift of the accumulated LoD values. In the case of the 1-day arcs, the linear drift for the GPS solution is at the level of −8.2 ms/year, which is 1.7 times larger than that for the GRE solution and even 11 times larger than for the GAL solution (Fig. 8). It is visible that the combined GRE solution is affected by the disruptive impact of the GPS solution. Meindl (2011) concluded that the different patterns of the linear drift for the GPS and GLONASS might originate from the applied approach to the SRP modeling. However, the differences between GAL and GAB solutions are minor. Thus, we infer that the difference in drift is partially caused by the orbital period and the arc length. In the case of GPS, the orbital period is in 2:1 resonance with the earth rotation and corresponds with the arc length. The corresponding resonances for GLONASS and Galileo are much weaker as they are represented by larger numbers and equal to 17:8 and 17:10, respectively. In the case of the 3-day arc solutions, individual system-specific solutions are more consistent with each other. The linear drift remains between − 0.7 and − 1.7 ms/year for GAL and GPS solutions, respectively. Thus, the extension of the orbital arc is beneficial, especially for the GPS solution.

Fig. 8
figure 8

Accumulated LoD values with respect to IERS-C04-14 for a the 1-day arc solutions and b the 3-day arc solutions. The dotted line denotes the linear regression line

Full size image

Figure 9 shows the spectral analysis of the LoD residuals with respect to IERS-C04-14 series. The LoD estimates delivered from the 3-day arc solutions have smaller STD by approximately 30% than those delivered from 1-day arc solutions (Table 6). The 1-day arc solutions are also more affected by spurious system-specific signals. However, the combination of GPS, GLONASS, and Galileo significantly decreases the amplitudes of most of these signals, most likely because of the variety of system characteristics. Same as for the PM rates (Fig. 7), the signals close to 3.4 days for Galileo and 3.9 days for GLONASS are visible in the 1-day solutions. For the 1-day Galileo-based solutions, the spurious signals with the amplitudes up to 8 µs/day are visible at the even harmonics of the draconitic year, i.e., 4, 6, 8, and 10 cpy. The switch from 1- to 3-day arc reduces most of these spurious signals. However, the 3 cpy signal with the amplitude of 7 µs/day remains for both GAL and GAB solutions. Finally, the signals close to 14.2 and 14.8 days are visible in all the solutions regardless of the GNSS constellation, especially when the 3-day arc is used. These signals originate from the errors in the conventional model of the sub-daily ERPs and aliasing of the O1 and M2 tidal terms (Griffiths and Ray 2013; Ray et al. 2017). The latest report given in 2019 by the IERS Working Group on the High-Frequency Earth Orientation Parameters recommends using the Desai and Sibois model with corrected tidal terms (Desai and Sibois 2016) instead of the model which has been recommended so far by the IERS conventions.

Fig. 9
figure 9

Amplitude spectra of the LoD residuals with respect to IERS-C04-14 for a the 1-day arc solutions and b the 3-day arc solutions. The vertical gray lines denote the harmonics of a draconitic year

Full size image
Table 6 Statistics of the LoD differences with respect to IERS-C04-14. STD denotes standard deviation
Full size table

Orbit misclosures and pole misclosures

The ERPs are modeled as linear between the particular midnight epochs within the orbital arc, assuming a conventional sub-daily motion. Thus, we may calculate the pole misclosures between the independent estimates from the successive days analogously as it is performed for the consecutive orbit arcs in the so-called orbit misclosures (Lutz et al. 2016). The quality of ERP misclosures derives directly from the quality of PM rates and LoD. Thus, the spectral analysis of the pole misclosures corresponds well with Fig. 7 and is not shown here. Similar to the PM rates, the quality dramatically improves when the 3-day arcs are generated instead of 1 day orbital arcs (Fig. 10). Because ERPs are the transformation parameters between the terrestrial and the celestial reference frames, ERP misclosures inevitably affect the orbit quality, especially close to the day boundaries.

Fig. 10
figure 10

Misclosures of the pole coordinates. The box ranges from 1st to the 3rd quartile. A vertical line inside the box reflects the median. The whiskers go from each quartile to the minimum and maximum value, excluding outliers, which are outside 1.5 times the interquartile range above the upper quartile and below the lower quartile

Full size image

Orbit misclosures may refer to the inertial (celestial) or the earth-fixed (terrestrial) frame (Lutz et al. 2016). Orbit misclosures, which are expressed in the inertial system, are integrally influenced by the ERP errors, including the PM and UT1-UTC misclosures, which depend on the errors in the a priori series, sub-daily models, and the quality of LoD estimation. Figure 11 illustrates the absolute position differences at the orbit boundaries in both inertial and earth-fixed frames, as well as the differences which arise from the ERPs. For the sake of consistency, the orbits and ERPs delivered in GRE solutions were used in the orbit transformation for all systems. First, because only the rotation is used, the orbit misclosures represented in the earth-fixed and inertial frames differ only in along-track and cross-track directions, while the radial direction remains unchanged. In the earth-fixed frame, the median misclosure is lowest for the GPS satellites compared to GLONASS and Galileo, in both 1- and 3-day arc cases. In the case of Galileo, the change of the SRP model from ECOM2 to the hybrid box-wing + reduced ECOM improves the median orbit misclosures by approximately 22 and 4% for 1-day and 3-day arcs, respectively. The magnitude of the impact of the ERPs on the orbit misclosures depends on the orbit altitude above the earth surface. Thus, the differences are greater for Galileo (23 225 km), than for GPS (20 200 km) and GLONASS (19 132 km) (Fig. 11c). In theory, a change in rotation of 100 µas corresponds to 12, 12.5 and 14.5 mm position errors at the altitudes of GLONASS, GPS, and Galileo, respectively. We may deduce that the orbit misclosures in the inertial frame are even better for GLONASS than for GPS and Galileo.

Fig. 11
figure 11

Orbit misclosures in a earth-fixed frame; b inertial frame; c impact of the ERPs on the orbit misclosures as seen by the difference of the orbit misclosures in both frames

Full size image

Conclusions

We have studied the ERP estimates, including X, Y pole coordinates, rates of the pole coordinates, and LoD, obtained from different cases of GNSS processing. Three main aspects have been addressed:

  • What are the differences in ERPs obtained from GPS, GLONASS, and Galileo solutions?

  • What is the benefit of using a combined multi-GNSS solution rather than system-specific solutions?

  • Is it advantageous for ERP estimation when a 3-day orbital arc is used instead of a 1-day arc?

All the PM estimates, which are obtained from different system-specific solutions, are systematically biased among themselves. We may notice that the offset changes with the development of respective satellite constellations, i.e., exchange of the old GLONASS-M satellites or completing the Galileo constellation. The remaining offset may originate from the system-specific issues, e.g., the lack of precise phase center offsets and variations of the receiver and satellite antennas for Galileo frequencies in our processing scheme. The reprocessing of Galileo data with absolute Galileo satellite and receiver antenna PCOs/PCVs chamber and robot calibrations remains a topic for future research, whenever the complete set of absolute receiver antenna calibration PCOs/PCVs become available for Galileo frequencies. Furthermore, the handling of inter-system offsets in ERP estimates should be also further investigated. The GPS-based PM estimates are of superior quality compared to those derived from GLONASS or Galileo. In the case of standard 1-day arcs, the scatter of the GPS-based PM residuals with respect to IERS-C04-14 are 1.6, and 1.2 times better for the X component, and 1.5 and 1.3 times better for the Y component, compared to GLONASS and Galileo, respectively. Considering oppositely, the drift of the UT1-UTC is up to 11 times larger for the GPS-based solutions than for the Galileo-based solutions, because of the strong orbital resonance.

The spurious signals inherently influence the Galileo-based and GLONASS-based ERPs at frequencies that arise from the resonance between the satellite revolution period and the earth rotation, e.g., 3.4 days for Galileo and 3.9 days for GLONASS. These signals dominate the series of pole coordinate rates and LoD. Moreover, the harmonics of the draconitic years, e.g., 183, 122, 91, 57 days, and aliasing periods of sub-daily ERP models, e.g., 14.2 and 14.8 days, are detectable in system-specific ERP estimates. Despite that all the system-specific solutions are affected by the spurious signals, the combination of different GNSS mitigates most of the uncertainties and improves the ERP results. Therefore, we recommend that the combined multi-GNSS solution should be used instead of system-specific solutions for ERP estimation purposes. Since the GLONASS solution delivered the worst results, the combination of GPS + Galileo needs further investigation. Moreover, the BeiDou constellation reached a number of 48 satellites in space in 2019. Using the combined constellation of MEO, IGSO, and GEO satellites may also decorrelate some of the systematic signals, which are still visible in the ERP estimates.

The extension of the orbital arc in the GNSS processing from 1- to 3- days is superior for the quality of the ERPs, especially for pole coordinate rates and LoD. The quality of the ERP rates also has a direct impact on the internal quality of the orbit in the inertial frame. Nonetheless, some of the artificial signals are amplified in the PM series of GLONASS-based and Galileo-based solutions. The amplitude of the 3 cpy signal is higher by a factor of 1.9 in the series of the GLONASS-based X pole coordinate when a 3-day arc solution is applied instead of 1-day. These issues may constitute a side effect of the deficiencies in orbit modeling for GLONASS and Galileo (Bury et al. 2020; Dach et al. 2019). If some significant orbit mismodeling exists, its effect might be amplified when using long arcs because the same number of parameters has to absorb errors accumulated over a long period. Finally, we infer that the 3-day arc is beneficial for GNSS results provided that the orbit is correctly modeled.

Motivated by Bury et al. (2019), who found an improvement in the precise Galileo orbit when a hybrid approach for SRP modeling is applied, we followed the same methodology and check the impact of SRP modeling on the ERP products. The GAB solution is slightly less scattered than the standard GAL solutions. The greatest improvement of approximately 10% is seen for the Y pole coordinate when cross-referencing GAB and GAL solutions. The amplitudes of the draconitic signals are also reduced, however not fully mitigated. Therefore, more research is needed to understand this matter and to further enhance the orbit modeling for Galileo.