Geocenter coordinates derived from multi-GNSS: a look into the role of solar radiation pressure modeling

The Global Navigational Satellite System (GNSS) technique is naturally sensitive to the geocenter motion, similar to all satellite techniques. However, the GNSS-based estimates of the geocenter used to contain more orbital artifacts than the geophysical signals, especially for the Z component of the geocenter coordinates. This contribution conveys a discussion on the impact of solar radiation pressure (SRP) modeling on the geocenter motion estimates. To that end, we process 3 years of GPS, GLONASS, and Galileo observations (2017–2019), collected by a globally distributed network of the ground stations. All possible individual system-specific solutions, as well as combinations of the available constellations, are tested in search of characteristic patterns in geocenter coordinates. We show that the addition of a priori information about the SRP-based forces acting on the satellites using a box-wing model mitigates a great majority of the spurious signals in the spectra of the geocenter coordinates. The amplitude of the 3 cpy (about 121 days) signal for GLONASS has been reduced by a factor of 8.5. Moreover, the amplitude of the spurious 7 cpy (about 52 days) signal has been reduced by a factor of 5.8 and 3.1 for Galileo and GPS, respectively. Conversely, the box-wing solutions indicate increased amplitudes of the annual variations in the geocenter signal. The latter reaches the level of 10–11 mm compared to 4.4 and 6.0 mm from the satellite laser ranging observations of LAGEOS satellites and the corresponding GNSS series applying extended empirical CODE orbit model (ECOM2), respectively. Despite the possible improvement in the GLONASS-based Z component of the geocenter coordinates, we show that some significant power can still be found at periods other than annual. The GPS- and Galileo-based estimates are less affected; thus, a combination of GPS and Galileo leads to the best geocenter estimates.


Introduction
The geocenter motion is conventionally interpreted as the movement of the center of mass of the earth system, including the solid earth, atmosphere, and oceans, with respect to the origin of the reference frame (Wu et al. 2012;Altamimi et al. 2016). This movement is expressed by a three-dimensional vector known as geocenter coordinates (GCC). Although the essence of the geocenter motion is straightforward, the direct observation of its instantaneous location is one of the most demanding applications of highprecision geodetic techniques and has a profound impact on the geophysical interpretation of geodetic measurements (Blewitt 2015;Kosek et al. 2020).

State of the art
The prospects for measuring geocenter motion with GNSS have been widely investigated (Lavallée et al. 2006;Meindl et al. 2013;Rebischung et al. 2014). Regarding the difficulties in the determination of the GCC using GNSS, the stateof-the-art knowledge points to some main aspects, among which two should be emphasized: (1) the reduction in the collinearity issues and (2) the reduction in the modeling deficiencies, especially in the aspect of the orbit model (Rebischung 2014). GNSS can determine the equatorial GCC components with good consistency to the other geodetic Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s1029 1-020-01037 -3) contains supplementary material, which is available to authorized users.  Arnold et al. 2015). On the contrary, the GCC-Z component is more difficult to determine using GNSS. The major source of the GCC-Z contamination is attributed to the orbit modeling deficiencies, especially the direct solar radiation pressure (SRP). With only limited a priori knowledge about the non-conservative forces acting on satellites, we must incorporate additional empirical orbit parameters into the solution, i.e., Empirical CODE Orbit Model (ECOM) or JPL GSPM (Springer et al. 1999;Bar-Sever and Kuang 2004). The errors in the orbit model, as well as the correlations between the estimated parameters, may introduce spurious orbit-related signals in the GNSS-based geophysical parameters, i.e., earth rotation parameters (ERP) and GCC (Meindl et al. 2013;Zajdel et al. 2020). Theoretically, an improvement may come from a decrease in the number of estimated parameters or, equivalently, from putting constraints on some of the estimated parameters. Bury et al. (2019) showed that using the boxwing model for the Galileo satellites allows us to decrease the number of the estimated ECOM parameters, which may improve the geocenter estimates. Moreover, the same set of ECOM parameters is not equally reasoned for the different GNSS satellites (Arnold et al. 2015;Prange et al. 2020a;Sidorov et al. 2020).
The second possible way to improve the GNSS-based GCC is the inclusion of additional decorrelating observations, i.e., combined multi-GNSS processing (Montenbruck et al. 2017) or the joint processing of GNSS observations collected by ground stations and LEO satellites Männel and Rothacher 2017). However, the processing of individual GNSS constellations brings us to different GCC, which complicates the successful multi-GNSS integration (Meindl 2011).

Outline
We investigate the quality of the GCC determination in different test cases of GPS, GLONASS, and Galileo processing for the period of 2017-2020. The time series of 3 years is sufficient to reliably estimate linear, annual terms (Blewitt and Lavallée 2002) and the draconitic errors. Since December 2018, the Galileo constellation is composed of 24 available satellites and is almost complete together with the legacy GPS (30 satellites) and GLONASS (24 satellites) systems. Zajdel et al. (2019a) analyzed the first results of the Galileo-based GCC. The errors for the Galileo-only GCC-Z were found to be lower by a factor of two when compared to the GLONASS-only estimates. Moreover, the Galileo-only GCC showed less artificial signals than the GLONASSonly, especially for the GCC-Z component. We first focus on the different approaches to the SRP modeling and then analyze the correlations between the empirical orbit parameters and GCC. Finally, we check the quality of the GCC estimates using various configurations of GNSS constellations, combining the existing GPS, GLONASS, and Galileo constellations in different variants, and including or excluding Galileo satellites accidentally launched into elliptical orbits (Sośnica et al. 2018). In summary, this work provides a comprehensive evaluation of the role of the SRP modeling in the determination of GCC using GNSS, which has never been done using three operational constellations employing 80 satellites.

Methodology
The computations are conducted using a development version of the Bernese GNSS Software . The processing details and the overview of the background models are described in detail by Zajdel et al. (2020) following the strategy used by CODE MGEX solution . The GCC are estimated simultaneously with the orbit parameters, station coordinates, troposphere parameters, and ERP. Noteworthy, the estimation of the troposphere parameters may influence the observability of the geocenter motion . However, that is not the main thread of this research as we may assume that the impact is the same, regardless of the test case. Finally, we test both 1-and 'overlapped' 3-day orbital arc solutions (Lutz et al. 2016). The longer orbital arc is tailored to reference frame realization and is designed to overcome some of the inherent weaknesses of the GNSS observations. However, the long arc might as well aggravate the solution quality if some significant orbit mismodeling exists ). In the case of the 1-day arcs, the stochastic pulses are estimated every midnight. For the 3-day arc solution, the additional stochastic pulses are also estimated at the day boundaries, which helps with the orbit fit over 3-day arcs (Prange et al. 2020b). Table 1 summarizes the naming of test cases, which differ in the approaches for the SRP modeling. We chose the different variations of the ECOM2 model (Arnold et al. 2015), which is employed by most of the analysis centers of the International GNSS Service (IGS; Johnston et al. 2017) for operational GNSS products, as well as the hybrid approach for SRP modeling as defined by Bury et al. (2020). The ECOM model 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 hybrid models combine the empirical parameterization with a priori information about the non-gravitational accelerations acting on the spacecraft based on the satellite dimensions and optical properties. Starting in 2017, the European GNSS Agency (GSA) has provided detailed metadata concerning the ground measured properties for the Galileo satellites (Bury et al. 2019). A similar set of information has never been published neither by GPS nor by GLONASS authorities for the full constellation. However, the properties of the box-wing model have been derived empirically or based on the assumed surface properties by Rodriguez-Solano et al. (2012). Springer et al. (2014) showed that the latter box-wing model for GPS and GLONASS might improve the GNSS-based products. The most recent and coherent information about the new GPS and GLONASS generations is provided in the frame of the IGS repro3 activities (http://acc.igs.org/repro 3/PROPB OXW.f). By the example of Galileo, the box-wing model accounts for approximately 97% of the direct solar radiation in the D direction (D 0 ) and also for a large part of the periodic terms D 2C and B C . However, the box-wing alone is insufficient without additional correction model parameters. Thus, the hybrid strategy of combining both specific terms of empirical and physical models is beneficial as for the orbit quality (Li et al. 2019) and the Galileo-based ERP .
To study in greater detail both single-system solutions and the combined multi-GNSS solutions, we employ a similar approach to that presented initially by Scaramuzza et al. (2018). One general NEQ is prepared for each day. The parameters, such as troposphere and station coordinates, are set up as common parameters. However, the GCC and ERP are set up as plane-specific. Then, the plane-specific parameters are stacked for any subset of either real or artificial GNSS constellations. This kind of processing allows us to flexibly study different constellation test cases as listed in Table 2, without the necessity of the preprocessing of each test case separately.
Primarily, we prepared the individual one-system solutions to investigate the system-specific signals in the GCC. The individual systems have also been grouped in all possible combinations forming multi-GNSS subsets ( Table 2). So far, the combined GPS + GLONASS time series has brought out the traces of the spurious GLONASS signals rather than contributing to the improvement in the GPSonly series (Meindl 2011). Given that the poor quality of the GLONASS orbit was blamed for this setback, we believe that the new light should be shed on this matter after the great improvements in the orbit quality of the multi-GNSS satellites (Zajdel et al. 2017;Sośnica et al. 2020).

Results
In this section, we discuss the GNSS-based geocenter motion estimates. Firstly, we provide an analysis of formal errors and correlations between the estimated parameters. Secondly, the individual single-system estimates of GCC are characterized. Finally, we discuss the GCC estimated in the multi-GNSS solutions.

Correlations between orbit parameters and geocenter coordinates
Over-parameterization of the orbit model may lead to the degradation of the solution due to correlations between the estimated parameters. The correlation coefficient (ρ) between the simultaneously estimated parameters provides information on the separability or the collinearity of the parameters. No significant correlations are found for the equatorial GCC components; thus, we limit the analysis to the GCC-Z only. As stated by Scaramuzza et al. (2018), the variability of the GCC-Z formal errors depends on the mutual orientation of the orbital planes with respect to the position of the sun (denoted as β angles). Therefore, we expect that the correlations between the parameters also vary in time. To verify this, 2 days have been selected for each constellation to reflect the epoch of the low and high GCC-Z  error. The selected days are highlighted in Fig. 1, which also shows the accompanying formal errors of the GCC-Z estimates. Note that the corresponding solutions, which apply the box-wing model or not (e.g., E2 and BX + E2), are equal in terms of formal errors and correlations between the parameters. Thus, only the E0, E1, and E2 solutions are discussed here as the representative series. For the epochs with the local maxima of the errors, significant correlations occur reaching ρ > 0.6. The correlations differ between satellites; thus, we chose one satellite per orbital plane. For the constellations with 3-nominal orbital planes such as Galileo and GLONASS, the GCC-Z is correlated with the periodic accelerations in the ECOM model, especially D 2C and B C (Fig. 2a). The correlation with the D 0 parameter is also visible, especially for the GLONASS satellite in the R-2 plane. At this specific epoch, the planes R-1 and R-3 have the same β angles (β ≈ 30°), while the R-2 complements the constellation with β ≈ -55°. In the case of Galileo, the correlations between GCC-Z and D 2C , B C , and D 0 are lower than for GLONASS. This applies to D 0 in particular. For the GPS constellation, almost all estimated orbit parameters seem to exhibit correlations 0.3 < ρ<0.6. The reduction in the number of estimated periodic terms (solutions E1 and E0) decreases most of the existing correlations (Fig. 2). The fewer the periodic terms are estimated in the ECOM model, the more separable are the GCC from orbital parameters in the GNSS solutions.
For the epoch, when the formal errors are low, we can assume that the observational geometry is favorable for the geocenter determination. Figure 2b shows that if the geometry of the observations is appropriate, no significant correlations between the orbit parameters and the GCC-Z Only one satellite per each orbital plane is selected. All the other parameters from the full variance-covariance matrix have been omitted; a the correlations at the epochs of maximum errors (Fig. 1); b the correlations at the epochs of minimum errors for the solution E2 (Fig. 1). Prefixes for the orbital planes: G-GPS, R-GLONASS, E-Galileo occur, even despite the estimation of the periodic empirical parameters in the ECOM model.

Signal decomposition of the GCC-Z estimates from the individual GNSS
In the next sections, we show the individual system-specific estimates of the GCC. We narrowed the discussion to the GCC-Z only, which is primarily affected. The complementary plots, which refer to the equatorial GCC, are provided as supplementary materials for the sake of completeness. In the figures contained in this section, the bold line comes from the low-pass-filtered series obtained with a Vondrák filter (Vondrák 1969) with a cutoff frequency corresponding to 40 days. Thereby, the low-pass fit reflects the largest-scale mass redistribution in the earth system visible as the annual and semiannual variations plus all the spurious draconitic signals up to 8th harmonic (about 43 days). Figure 3 illustrates the GCC-Z delivered by the GPS constellation. When the ECOM2 model is used solely (E2), the artificial signal is visible with the period close to the 7 cpy (about 52 days). This pattern has already been indicated by Rodriguez-Solano et al. (2014), who found that the eclipsing period of consecutive GPS orbital planes falls on average every 52 days. Figure 4 shows that the β zero-crossing points for the consecutive orbital planes repeat for GPS not only every 52 days, but also 44 and 123 days. Even if eclipses last at a maximum of 1 h within one GPS satellite orbital revolution, it may degrade the orbital arc and deteriorate the GNSS-based products. Using the box-wing model appears to improve the orbit modeling in the eclipsing season and quieten the spurious signal at the 7th cpy. The amplitudes of the annual signal are equal to 7.9, 8.5, 10.8, and 9.9 mm for the E2, BX + E2, BX + E1, and BX + E0 solutions, respectively. In the case of the box-wing solutions, the amplitude of the annual signal is larger than for the E2 solution. It may be caused by the errors in the box-wing model, which is only empirically adjusted best approximation of the GPS spacecraft due to the lack of official information on satellite surface properties. Due to the limited length of the time series, an annual and draconitic period cannot be separated (Männel and Rothacher 2017). However, the discussed increase in the annual amplitude is only due to the handling of SRP modeling. This is why the increase in the annual amplitude should be associated with the orbital error at the period of the GPS draconitic year (about 351 days).

GCC-Z from GPS
Using the 3-day length arc is beneficial for the GPS-based geocenter estimates, as in line with the expectations (Lutz et al. 2016). The amplitude of the annual signal is lower in general, except for solution BX + E0, for which the factor is 1.7. Apparently, the simplified box-wing model of the GPS satellites plus constant ECOM parameters is not enough to handle all the non-conservative forces properly for the GPS satellites.

GCC-Z from GLONASS
GLONASS satellites are known for introducing extreme spurious signals into the estimates of the geocenter motion (Meindl et al. 2013;Lutz et al. 2016). Figure 5 shows the GCC-Z delivered by the GLONASS constellation. The suspicious 3 cpy signal in the GLONASS-based GCC-Z signal reaches the amplitude of approximately 90 mm in the E2 solution. However, the 3 cpy signal in the GLONASS solution can be significantly reduced when the a priori box-wing is applied. The amplitudes are reduced by a factor of 3.2, 2.1, and even 8.5 for the BX + E2, BX + E1, and BX + E0, respectively. It should be recalled that the box-wing model assumes the simplified shape of the satellite bus. But as we can see, even such a simplification helps with the mitigation of the systematic errors in the GCC-Z estimates. In addition to the 3 cpy dominating spectral line, the spectra also show the pronounced peaks at the odd harmonics of the draconitic GLONASS year, i.e., 7th (about 52 days) and 5th (about 70 days) harmonic. These periods could also be identified as the repeating period of both β zero-crossing and equal β points for GLONASS (Fig. 4). These are mitigated in the solution BX + E0 to a large extent.
The 3-day arc has only a minor impact on the GCC-Z estimates and on both 1-and 3-day arc series roughly correspond with each other. However, the annual signal is increased when the E2 set of ECOM parameters is estimated for GLONASS. The increase reaches the factor of 1.1 and 1.4 for BX + E2 and E2 solutions, respectively.

GCC-Z from Galileo
Galileo reached the number of 24 active satellites at the beginning of 2019 and consists, similarly to GLONASS, of three nominal orbital planes. Figure 6 presents the GCC-Z delivered by the Galileo constellation. For Galileo, when the E1 model is applied, the 3 cpy signal dominates the GCC-Z series, which is the same as we have seen for GLONASS (Meindl et al. 2013). The switch from E1 to E2 remarkably improves the GCC-Z estimates. The amplitude of the 3 cpy signal decreases from 101 to 4 mm. However, a pronounced signal remains with a period close to 50 days and an amplitude of 8 mm. Noteworthy, the interval of 52 days is also addressed to the period between the epochs when the two out of three Galileo orbital planes have the same orientation in the sun-earth-satellite frame (equal β angles) (Fig. 4). The GCC-Z signal significantly changes under the influence of the a priori box-wing model. The GCC-Z as visible in the solutions BX + E1 and BX + E0 resembles each other and contains much less orbital artifacts than the ECOM-only solutions.
The influence of the 3-day arc on the Galileo-based GCC-Z estimates is different than on the GLONASS solutions. For Galileo, the amplitude of the annual signal increases when the E1 or E0 set of ECOM parameters is used; and hence, the twice-per-revolution accelerations in D direction are not estimated. The increase reaches the factor of 1.3, 1.3, and 2.4 for the BX + E0, BX + E1, and E1 solutions, respectively. Table 3 shows the variability of the X, Y, and Z geocenter coordinates with the standard deviation (STD) as a criterion. The two sets of metrics have been provided. The first set refers to the raw estimated signal. For the second set, the low-pass filter fit with the cutoff frequency of 40 days has been subtracted from the original series forming the geocenter residuals. The impact of the short-scale geophysical tidal and non-tidal signals with the repetition periods less than 40 days is minor for the geocenter motion. Thus, the high-pass STD may be treated as the noise indicator. The statistical nature of the time-correlated noise in the time series of GNSS-based GCC is characterized by Ma et al. (2020). Accordingly, the noise in the GCC is best modeled by either white plus power law or white plus generalized Gauss Markov (GGM) noise models. Table 4 gives the same set of metrics as Table 3; however, it refers to the 3-day arc solutions, instead. Table 3 confirms that the change in SRP modeling has only a minor impact on the equatorial GCC. By the example of BX + E1 solutions, STD of the X(Y) geocenter coordinates is equal to 4.8(6.2), 6.3(7.8), and 7.8(9.7) mm for the GPS, GLO, and GAL solutions, respectively. However, STD of the high-pass part of the X(Y) geocenter coordinates is equal to 3.8(4.1), 5.3(5.7), and 6.6(7.4) mm for the GPS, GLO, and GAL BX + E1 solutions, respectively. Thus, after the reduction in the low-pass part, approximately 70-80% of uncertainty remains, which is a huge limitation of the GNSS-based geocenter estimates.

Variability of the geocenter coordinates
Comparison of the corresponding values in Tables 3 and 4 shows that the 3-day solutions for the geocenter are, in essence, a smoothed version of their 1-day equivalent. The high-pass STD is generally reduced by up to 30%, 30%, and even 40% for the GPS, GLO, and GAL solutions, respectively. The higher reduction for the GAL solution is, to a large extent, attributed to the mitigation of the Galileorelated orbital signals close to 3.4 days (Fig. 6) (Zajdel et al. 2019b. The 3-day arc has naturally smoothened the signals with periods shorter than 6 days. Thus, for GLONASS, the orbital signal, with a period close to 8 days, remains in both 1-and 3-day arc solutions. Figure 7 illustrates the time series of the high-pass residuals of the GCC-Z estimates. For Galileo, the increase in the residual variability from about 10 mm up to 60 mm is visible for certain epochs. The intensive noise epochs correspond with the increased formal errors of the GCC-Z estimates for Galileo (Fig. 1). The pattern is mitigated along with the reduction in the periodic accelerations in the ECOM model. Therefore, the theoretical improvement, as indicated by the formal errors and parameter cross-correlations, is also reflected in the high-pass noise of the Galileo-based GCC-Z estimates. The latter is visible in both 1-and 3-day arc solutions. Such patterns are not that distinct for GPS solutions. For GLONASS, the high-pass noise is also periodically amplified with the pattern of the same principle as for Galileo. However, the increase is less evident, compared to Galileo, and ranges from about 20 to 60 mm for the worst E2 solution. Noticeably, the general high-pass noise is larger for GLO solutions than for GPS and GAL solutions. The reduction in the high-pass STD between the worst (E2) and the best (BX + E0) cases reaches up to 44, 67, and 64% for GPS, GLO, and GAL 1-day arc solutions, respectively. In the best test-case for the BX + E0 3-day arc solution, the high-pass STD of the GCC-Z equals 2.6, 3.9, and 6.1 mm for the GPS, GAL, and GLO solutions, respectively. During the considered period, the Galileo constellation had grown from 14 satellites at the beginning of 2017 to 24 available satellites in 2019. The increase in the number of observations and improvement in solution geometry lead to the reduction in high-pass STD. In the case of the BX + E0 1-day arc solution, the high-pass STD decreases from 7.3, 7.9, and 8.9 mm to 5.2 (28%), 7.1 (10%), and 4.3 (52%) mm for the X, Y, and Z components of the GCC, respectively.

Impact of the Galileo satellites in the eccentric orbits on the geocenter coordinates
The Galileo constellation is nominally distributed over three orbital planes separated by 120° in the right ascension of the ascending node. However, due to the launch failure of the two first Galileo FOC satellites (E14 and E18), the fourth highly eccentric orbital plane has been formed (Sośnica et al. 2018). Meindl et al. (2013) assumed that the GCC-Z from Galileo and GLONASS should be affected by similar artifacts due to the 3-plane geometry. The results from the previous sections of this work partly denied this statement. However, one may say that Galileo shall be considered as a 4-plane, rather than a 3-plane constellation. Thus, we decide to check whether the inclusion of the fourth Galileo orbital plane with just two eccentric satellites anyhow affects the GCC-Z estimation. Figure 8 shows that the two additional Galileo satellites have an impact on the GCC-Z estimates, especially in the E2 solution. When the eccentric orbit is included, the amplitudes of the spurious orbital signals at the odd harmonics of the Galileo draconitic year are reduced by factors 1.4, 1.5, and 2.5 for the 7th, 5th, and 3rd harmonics, respectively. In the cases of the BX + E1 and BX + E0 solutions, the impact of the additional Galileo satellites is rather minor. Therefore, the example of the E2 solution demonstrates how the improvement in the observational geometry may help in the GCC-Z estimation when the significant errors in the orbit modeling occur. Figure 9 illustrates the formal errors and the high-pass noise of the GCC-Z estimates in the Galileo solution with and without the satellites on the eccentric orbits for the worst case-E2. The formal errors have dropped by more than a half after the inclusion of the two additional satellites on the eccentric orbits. When only three nominal Galileo planes are considered, the pattern of the GCC-Z formal errors corresponds to those delivered in GLONASS E2 solution (Fig. 1). The standard deviation of the high-pass signal is also reduced from 24 to 18 mm (Fig. 9b). This is mainly due to a decrease in noise close to epochs that correspond to the maxima of the GCC-Z formal errors. In the case of the BX + E1 and BX + E0 solutions, the reduction in formal errors is at the submillimeter level. Fig. 7 Residual high-frequency estimates of the GCC-Z (with adjusted low-pass Vondrák filter with a cutoff frequency corresponding to 40 days) in mm. The colored and gray series refer to the 1-and 3-day arc solutions, respectively. The dashed gray lines denote the sun elevation angles above orbital planes (β) 1 Page 10 of 15

Geocenter coordinates from the multi-GNSS solutions
Finally, the question of whether the combination of the individual GNSS constellations improves the geocenter estimates should be addressed. We examine the combined GCC in the 1-day arc approach only as it should emphasize the potential deficiencies in the particular solutions. Figure 10 shows the GCC-Z as estimated in different combinations of the GNSS constellations. On the first hand, if we include the system in the combined solution, the specific errors, which arise from the SRP modeling applied for this system, affect the combined solution as well. On the other hand, in the combined multi-GNSS solutions, we may notice that the amplitudes of the artificial signals are even a few times lower than in the case of the single-system solutions. All solutions that include the Galileo system and the E2 model are affected by the signal with a period close to 52 days. If GLONASS is included and the periodic ECOM parameters are estimated, the GCC-Z solution is affected by the 3 cpy signal.
The combination of the GLONASS and Galileo constellation leads to the worst result, especially in the E2 and BX + E1 cases. The latter is caused mainly by GLONASS, as the corresponding Galileo-only GCC-Z estimates were of good quality (Fig. 6). However, when reducing the estimated set of the ECOM parameters to the constant accelerations only, we see that most of the orbital signals have been reduced and the annual signal dominates as expected from theory. Figure 11 illustrates the amplitudes and phases of the annual sine waves estimated from all the discussed combined solutions, as well as the representative single-system solutions of BX + E0. The signal of the geocenter motion is nonstationary, and the dominant annual signal is characterized by the apparent variable amplitude ). Therefore, a validation with an independent dataset covering the same period is more justified than the comparison with the state-of-the-art results. The reference series has been prepared in the independent processing of the SLR observations to LAGEOS-1 and LAGEOS-2 satellites (Zajdel et al. 2019a).
For the GNSS-based GCC-X, the annual amplitudes are on average in good agreement with the reference series. In the case of the GCC-Y, the amplitude of the annual signal reaches 4-5 mm, which is more than two times larger than for the SLR series. Noteworthy, such a difference between the SLR and GNSS GCC-Y could be, to a large extent, attributed to the distribution of the GNSS and SLR sites. Zajdel et al. (2019b) showed that the selection of datum defining sites in the SLR processing might change the amplitude of the annual signal of the GCC-Y in the range between 2 and 4 mm. On the other hand, the phases of the SLR and GNSS annual signals differ only slightly by about 10° on average. Thus, it confirms in essence that both techniques sense the same geophysical signal. The phases of the annual signal in the GAL BX + E0 are shifted by approximately 45° with respect to the remaining series. However, one should note that the most up-to-date information about the phase center calibrations of the Galileo satellites and ground antennas is not taken into account within this processing (Villiger et al. 2020). Whether the change of this processing feature would improve or not the GCC estimates should be further investigated.
There are some large differences in the amplitudes and phases of the annual signals in the GCC-Z series (Table 5). First, the solutions, which include the box-wing model, reveal the increase in the annual amplitude by a factor of up to 1.7 with respect to the purely empirical E2 solution. For the latter, the annual amplitude is also generally in the best agreement with the SLR series. An exception is the R + E solution, for which the annual amplitude in the E2 solution is also at the level of 10 mm. The phase of the annual signal in all the GNSS-based solutions is lower by approximately 30° with respect to the SLR series. The variations in phases between the individual test cases are by far more coherent than those reported by Ma et al. (2020), who analyzed the GCC estimates delivered by different IGS ACs in the frame of the IGS repro2. However, many of the IGS contributions to the repro2 applied the ECOM1 model for the GLONASS orbit model, which may substantially deceive the metrics of the annual signal.

Conclusions and discussion
The GCC estimated using single GNSS systems are affected by the artificial system-specific signals. These signals may prevail, to varying extents, over the desirable geophysical signal of the geocenter motion. The degree of quality loss depends on the selected orbital arc length and,  to a greater extent, the approach to SRP modeling. The separability of the GCC in the GNSS processing depends on the number of empirical ECOM parameters being estimated. The fewer the empirical parameters are in the orbit model of the solution, the more independent the GCC-Z parameters. The GCC-Z is correlated not only with the D 0 and B C (with ρ > 0.6) but also with the D 2C empirical parameter. The latter concerns mainly GLONASS and Galileo during the specific epochs when two out of three orbital planes are similarly oriented with respect to the sun (β). Finally, we identify the G + E-based GCC solution as the most reliable among all the other possible single-system and multi-GNSS solutions. The great majority of the artificial GPS-specific and Galileo-specific signals are mitigated, while the annual signal is dominant. Figure 12 presents all the three geocenter components as delivered in the G + E solutions to visualize the achieved results. The E2 solution is more consistent with the SLR in terms of the annual signal. However, the signal close to the 52 days corrupts the GCC-Z estimates. On the other side, the box-wing greatly reduces the harmonic draconitic signals in the BX + E0 and BX + E1 solutions. However, the annual amplitude of the GCC-Z is suspiciously too large compared to those estimates from SLR. The improvement in the box-wing model should be the key point for the clarification on this aspect. Because the draconitic period is close to one solar (synodic) year, the draconitic signal modulates the annual signal and obscures the fingerprints of periodic geocenter motions along the Z-axis. The sufficiently long time series would allow recovering the amplitude and phase of both solar and draconitic year; however, the Galileo series is still too short.
Eventually, the estimation of GCC shall benefit from tighter constraints on the GNSS satellite clocks when employing zero-differences instead of double-differences of GNSS observations. The Galileo clocks are more stable than those available in GPS and GLONASS (Kazmierski et al. 2017). Future improvements in the field of GNSS clock modeling should emerge as an important topic of research regarding the geocenter motion determination in GNSS. A potential release of satellite metadata for GPS and GLO-NASS would allow constructing more accurate box-wing models, which would also be beneficial for the improvement in the geocenter determination, as observed for Galileo.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.