Determination of precise Galileo orbits using combined GNSS and SLR observations

Galileo satellites are equipped with laser retroreflector arrays for satellite laser ranging (SLR). In this study, we develop a methodology for the GNSS-SLR combination at the normal equation level with three different weighting strategies and evaluate the impact of laser observations on the determined Galileo orbits. We provide the optimum weighting scheme for precise orbit determination employing the co-location onboard Galileo. The combined GNSS-SLR solution diminishes the semimajor axis formal error by up to 62%, as well as reduces the dependency between values of formal errors and the elevation of the Sun above the orbital plane—the β angle. In the combined solution, the standard deviation of the SLR residuals decreases from 36.1 to 29.6 mm for Galileo-IOV satellites and |β|> 60°, when compared to GNSS-only solutions. Moreover, the bias of the Length-of-Day parameter is 20% lower for the combined solution when compared to the microwave one. As a result, the combination of GNSS and SLR observations provides promising results for future co-locations onboard the Galileo satellites for the orbit determination, realization of the terrestrial reference frames, and deriving geodetic parameters.


Introduction
Modern satellites of the global navigation satellite systems (GNSS) are commonly equipped with laser retroreflector arrays (LRA) for satellite laser ranging (SLR). As a result, two independent space techniques are co-located onboard navigation satellites. This allows investigating the orbital quality using SLR residuals as well as performing an independent GNSS orbit solution based solely on laser observations (Pavlis 1995). Moreover, GNSS orbits can be determined using the combination of two types of observations (Thaller et al. 2011;Hackel et al. 2015).
The International Laser Ranging Service (ILRS, Pearlman et al. 2019) provides recommendations for the SLR stations in the form of the priority list of targets to be tracked. Nowadays, not only the dedicated geodetic satellites are tracked, but also plenty of other scientific missions including GNSS spacecraft. In 2013, the ILRS established a special study group, called Laser Ranging to GNSS s/c Experiment (LARGE) to organize special GNSS tracking campaigns. In 2018, two GNSS tracking campaigns were organized, between February 15 and May 15, and from August 1 to October 31. The results of the second campaign are especially pronounced as eight Galileo satellites were recommended for tracking by the ILRS stations (see Fig. 1).
The research on the combination of microwave GNSS and SLR observations to GNSS satellites was mostly focused on the realization of the terrestrial reference frame (TRF) (Thaller et al. , 2015Glaser et al. 2015;Bruni et al. 2018). Urschl et al. (2007) assessed the formal errors of the GPS orbit semimajor axis based on sparse SLR data and simulated continuous SLR observations to GPS and GLONASS satellites. When the weights for the SLR and GNSS observations are the same, no improvement of the a priori formal error of the semimajor axis was achieved (Urschl et al. 2007). When the weight for SLR increases Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s1029 1-020-01045 -3) contains supplementary material, which is available to authorized users.
by an order of magnitude, the formal error improves; however, the number of four SLR sites is insufficient to determine an accurate orbit solution. Schönemann et al. (2007) and Hidalgo et al. (2008) incorporated SLR observations to determine the very first Galileo In-Orbit Validation Element (GIOVE-A and -B) orbits using combined SLR + GNSS solutions. In both studies, the combined solution was characterized by an improvement of the orbit overlap. However, in both cases, the number of GNSS stations tracking Galileo satellites was limited to 13. In such a situation, any increase in the number of stations would be beneficial in terms of the improvement in the geometry of observations. As a result, the question arises, what would be the contribution of the SLR observations to the GNSS precise orbit determination (POD), if the GNSS data were provided by approximately 100 globally distributed stations, and the satellites were tracked by over a dozen of the SLR stations, as it is the case nowadays. Moreover, both studies applied higher, by an order of magnitude, weights for the SLR observations to GNSS satellites than for microwave GNSS observations. That is understandable because the orbit modeling issues were not resolved (Urschl et al. 2005), and the accuracy of early results for the GIOVE satellites was at the decimeter level (Steigenberger et al. 2011). Currently, the accuracy of GNSS and SLR techniques is significantly better, i.e., the a posteriori error of the multi-GNSS double-difference carrier phase solutions is at the level of 1.5 mm . Therefore, the weighting strategy should be reconsidered on that basis. Hackel et al. (2015) conducted the Galileo In-Orbit Validation (IOV) orbit determination using combined GNSS and SLR observations. The internal consistency of the combined solution deteriorated in terms of, e.g., day boundary misclosures, compared to the GNSS-only solution. The solution consistency was improved by shifting the offset of SLR retroreflector by 5 cm, which should rather be accounted for by the estimation of the SLR signature effect which depends on the laser incidence angle, and type of the LRA mounted on a satellite (Sośnica et al. 2018).
The goal of this study is to develop the best strategy for the precise Galileo orbit determination using combined GNSS and SLR-to-GNSS observations. We check the influence of SLR observations to GNSS satellites combined with microwave GNSS (GPS + GLONASS + Galileo) data collected by 100 GNSS stations and a dozen of SLR stations. We compare the combined solutions to microwave-based GNSS orbits whose accuracy is already at the centimeterlevel . We also test different weighting strategies to develop the optimal POD strategy using combined GNSS and SLR solutions.

Methodology
We conduct a series of Galileo POD test strategies to verify the impact of SLR observations on the combined SLR-GNSS solutions for 2017.0-2019.0. As a reference solution, we determine purely microwave-based Galileo orbits (M1) based on ionospheric-free linear combinations of doubledifferenced GNSS phase observations collected by globally distributed stations, which follows the strategy described by Bury et al. (2019b) and Zajdel et al. (2019aZajdel et al. ( , 2020. For the absorption of the direct solar radiation pressure (SRP), for all solutions, we use the hybrid approach based on the box-wing model for the Galileo satellites ) and the empirical parameters consistent with the 5-parameter Empirical CODE Orbit Model (ECOM, Springer et al. 1999).
The 1-day normal equations (NEQs) for the microwave solutions are the basis for the combining SLR and GNSS solutions. The outliers in GNSS data are rejected in a preprocessing step. The GNSS-based NEQs are calculated using GPS, GLONASS, and Galileo observations provided by approximately 100 multi-GNSS stations distributed homogeneously worldwide. For SLR observations to GNSS satellites, NEQs are generated consistently with the strategy described by Bury et al. (2019a). Contrary to GNSS, the NEQs for SLR are not inverted at the preparation level due to their singularity, but instead, the observation residuals are rejected based on fixed orbits and station coordinates. We remove the SLR observations with residuals exceeding the threshold of 90 mm. In contrast to the GNSS NEQs that are calculated using double-differenced phase observations, the SLR NEQs are calculated based on direct, absolute range measurements. The number of SLR observations is overwhelmingly lower as compared to GNSS, i.e., the median value of the GNSS observations for one day is about 882,809 for GNSS as compared to the 371 SLR observations for GLONASS and Galileo together. Eventually, the combination of the two types of observations demands to cope with different models dedicated to the particular technique. It is essential that both GNSS and SLR NEQs must be prepared using the same background models for orbit integration to avoid any sources of inconsistencies (see Table 1).
The combination of GNSS and SLR data is done at the NEQ level using the independent 1-day NEQs for GNSS and SLR with different weighting strategies. The a posteriori sigma of unit weight is at the level of 1.5 mm for the microwave Galileo solutions. For the SLR technique, the number of SLR observations to the Galileo satellites is insufficient to invert the SLR-only NEQ and obtain the a posteriori sigma. As a consequence, the accuracy of laser ranging to GNSS is established based on the SLR residuals to GNSS satellites using an online tool for the GNSS orbit validation with SLR (https ://govus .pl/slr, Zajdel et al. 2017). For the best-performing stations, the standard deviation (STD) of the SLR residuals is at the level of 15 mm; hence this value is assumed to be an indicator of the SLR data quality (W1). We also calculate the solution, for which the GNSS and SLR observations are treated equally (W3). Moreover, we provide an intermediate case with the mid-range weights for the SLR observations to test the effects of increasing the impact of SLR data (W2). The weights for GNSS and SLR NEQs in this particular strategy is obtained using the −2 values from Table 2.
The datum is realized based on two sets of the GNSS and SLR stations. A priori GNSS station coordinates are expressed in the IGS14 reference frame, whereas the SLR site coordinates are referred to SLRF2014 all of which are consistent with the realization of the International Terrestrial Reference Frame (ITRF2014, Altamimi et al. 2016). Strugarek et al. (2019) proved that the differences between IGS14 and SLRF2014 are negligible and do not introduce any systematic effects to the combination of SLR and GNSS data onboard Sentinel-3A/B satellites. For the combined network of GNSS and SLR stations, we use minimum constraints, i.e., no-net-rotation (NNR) and nonet-translation (NNT) for both networks, as we simultaneously estimate the earth rotation parameters (ERPs) and geocenter coordinates (Table 3). For the GNSS technique, the minimum constraints are applied to the stations which were verified using residuals after Helmert transformation for each day of the solution , whereas for the SLR, we use the ILRS core stations (https ://cddis .nasa.gov/lw21/docs/2018/splin ter/ASC_CANBE RRA_2018_PRESE NTATI ONS.pdf). From the core station list, we rejected stations which may deteriorate the SLR solution because of station-specific systematics as recommended by Zajdel et al. (2019b). The ILRS core station list is prepared based on LAGEOS tracking, hence not all the core stations contribute to the GNSS tracking. In the best case, 13 SLR stations tracked Galileo satellites during the 1-day solution.
The NEQs, including two types of observations, i.e., the GNSS and SLR, are stacked to one NEQ for each day in program ADDNEQ2 of a modified version of the Bernese GNSS Software . The modifications of the Bernese GNSS Software include the extended handling of network constraining separately for SLR and GNSS stations, Galileo data processing using the a priori box-wing model based on Galileo metadata (GSA 2016), and optimized data processing and screening for SLR observations to GNSS .
Apart from the station coordinates and the GNSS orbit parameters, which include six Kepler elements and five empirical parameters, we estimate the ERPs, the geocenter coordinates, troposphere zenith path delays, troposphere gradients, and SLR range biases (RBs). We calculate the annual mean RBs for each satellite-station pair (Appleby et al. 2016;Sośnica et al. 2019) which are then resubstituted as a priori values and strongly constrained in the processing. RBs treated in such a way do not deteriorate the orbit solution ).

Results
We check how the addition of the SLR observations influences the orbit geometry by analyzing the semimajor axis and its formal errors. The POD accuracy is assessed using the SLR residual analysis. We also analyze, how the SLR observations influence the satellite position, as well as the impact of SLR data on the empirical orbit parameters. Finally, we analyze  the systematic drift of the accumulated Length-of-Day (LoD) parameter, also known as LoD bias (Senior et al. 2010), as it is correlated with the drift of the satellite ascending node and may serve as an additional orbit quality indicator. We also investigate RBs as a prerequisite for a reliable SLR solution and analyze the quality of the 2-day orbit predictions. The two later analyses are attached as supplementary material.

Formal error of the semimajor axis
The analysis of the semimajor axis formal error allows us to assess the impact of the addition of Galileo SLR observations on the solution precision. We check whether the small number of SLR observations, when added to the processing of GNSS observations, diminishes the formal error, and if so, we investigate whether there is any dependency between the number of the SLR observation and the reduction of the semimajor axis formal error for different strategies. In all test cases, when we added SLR observations to the GNSS satellites, the formal errors of the semimajor axis are smaller as compared to the values from the microwave-only solution (Fig. 2). This is expected as the addition of a different type of observations to NEQs improves the observational geometry in the functional model. Even though the number of SLR observations is much lower than that for the GNSS technique, the impact of laser ranging is noticeable. SLR provides direct range measurements to the GNSS satellites, mostly in the radial direction. As a result, the SLR observations strengthen the solution thanks to their advantageous characteristic of direct range measurements, as well as improve the observation geometry.
The impact of the SLR technique increases with the number of observations contributing to the particular 1-day solution. Regardless of the strategy, the more SLR observations, the higher the reduction of the semimajor axis formal error. Moreover, the reduction increases when the SLR observations have greater weights. The addition of 50 SLR observations to the GNSS solution diminishes the mean error by 6, 13, and 25% for the Galileo-FOC, and 5, 15, and 37% for the Galileo-IOV for strategies W1, W2, and W3, respectively.
The formal error of the semimajor axis also depends on the sun-satellite geometry. Figure 3 illustrates the reduction of the formal error of the semimajor axis as a function of the β angle which is the elevation of the Sun above the orbital plane. We chose here the Galileo satellites from plane C that is characterized by the highest maximum β, i.e., 78°. The formal errors depend on β in the microwave solutions in such a way that the largest errors are for maximum values of |β| and for |β| close to zero. When introducing SLR observations to GNSS satellites, the β-angle dependence diminishes for both IOV and FOC satellites. The mean formal error of the semimajor axis diminishes by 5, 12, and 23% for W1, W2, and W3 solutions, respectively, for |β|= 78°.

Analysis of the orbit misclosures
Orbit misclosures evaluate the consistency level between the orbital arcs of consecutive days. The differences between  Table 4 shows the orbit misclosure statistics for all strategies.
Contrary to the analysis of the semimajor axis formal error, the orbit misclosures are smallest for the SLR observation weights from solution W1. For the radial direction, which is measured directly by the SLR technique, the mean value of the orbit misclosures becomes smaller with the increase of the SLR weight for both FOC and IOV satellites. However, the STD and mean value of misclosures is higher for solutions W2 and W3 in both along-track and cross-track directions. Although the formal error of the semimajor axis decreases with the strengthening of the impact of SLR observations, the orbit consistency suffers when the SLR observations have higher weights in W3 compared to W1, which reflects best the accuracy of the two techniques. For the solution W1, in all directions, the STD of the orbit misclosures is lower than for M1, which is reflected in Fig. 4 by the accumulation of more solutions around the zero value. Although in the along-track direction, the absolute mean values of the orbit misclosures are higher in W1 than in M1 and with a negative sign, the orbit solution is more consistent, i.e., the STD is lower for the combined solution. According to Bury et al. (2019a), the least accurately determined GNSS orbit component is the cross-track direction when based solely on SLR. However, using the combined GNSS and SLR solutions, the cross-track component, as well as the other directions, are improved by adding SLR observations to microwave-only solutions. When properly weighted, the STD of the orbit misclosures decreases from 47 to 46, 83 to 75, and 67 to 61 mm in the radial, along-track, and cross-track directions, respectively, from the M1 to the W1 solution.

SLR residual analysis
SLR residual analysis represents an external validation tool for microwave-derived orbits. In this study, however, we calculate both the microwave and hybrid solutions, which are based on both observation types, thus, the SLR validation is not fully independent. Table 5 shows the statistics of SLR residuals for all the Galileo-IOV and FOC satellites. The STD of SLR residuals decreases when the SLR observations are added to the solution, i.e., from 26.6 and 29.3 to 25.0 and 26.8 mm in the case of W1 for FOC and IOV satellites, respectively. When the impact of SLR is increased in solution W2, both the offset and STD of SLR residuals further decrease for all Galileo satellites, i.e., to 24.4 and 25.5 mm for FOC and IOV, respectively. However, when the weight of the SLR observations is equal to GNSS (W3), the STD of SLR residuals is higher than for W2. For strategy W3, the STD of the SLR residuals increases by 8% compared to W2. The decreasing offset for the combined solution may suggest that the solutions based on the two techniques become dominated by SLR (Fig. 5). Although the offset for the Galileo-FOC diminishes to almost zero, the solution W3 seems to be incorrectly dominated by SLR observations. The Galileo-IOV orbits suffer from systematic effects when |β| is higher than 60°. This spurious effect in Galileo IOV orbits was also detected in the experimental microwavebased IGS combined orbits ). The application of the SLR minimizes these β-dependent effects, as the SLR observations directly control the radial direction of the orbit. As a result, the STD of residuals for the IOV satellites while |β|> 60° decreases from 36.3 to 29.6 mm for the solutions M1 and W2. For high β angles, the orbital direction B of the ECOM model coincides with the radial direction that is measured directly by SLR, which means that the effect should also be visible in the B 0 parameter. Figure 6 illustrates the SLR residuals as a function of the sun-satellite-earth geometry, i.e., the β angle and the argument of the latitude of the satellite with respect to the argument of the latitude of the Sun (Δu). Considering SLR observations in POD diminishes the STD of the SLR residuals by approximately 8 and 13% for the FOC and IOV, respectively, when compared to the solution M1. The positive impact of the SLR residuals for the IOV satellites for |β|> 60° is visible in the whole range of Δu. Figure 7 illustrates the percentage share of SLR residuals divided into 1-cm intervals. When adding the SLR observations, the number of the SLR residuals smaller than 1-cm increases for strategies W1 and W2. The highest percentage share of the residuals below 2 cm for the FOC satellites is obtained for solution W2 and comprises 71.7% of all observations. For solutions M1, W1, and W3, the SLR residuals below 2 cm comprise 60.6. 66.0, and 68.5%,  respectively for FOC satellites. Moreover, the lowest number of the SLR residuals exceeding 5 cm is also obtained for strategy W2, i.e., 3.8%. For IOV satellites, the strategy W3 might look more consistent than the solution W2; however, the number of SLR residuals exceeding 5 cm for FOC satellites is higher for W3 than for W2. Solution W2 seems to provide the most reasonable combination of SLR and GNSS observations; therefore, the further analyses will concentrate on W2.

Impact of the SLR observations on empirical orbit parameters
The empirical orbit parameters typically reveal the deficiencies in the applied force models. In the previous section, we found that SLR substantially improves the Galileo-IOV orbit quality for |β|> 60° when the B 0 term coincides with the radial component measured directly by the SLR technique. Figure 8 shows differences between estimated empirical parameter B 0 from strategies M1 and W2 for Galileo-IOV E19 and Galileo-FOC E08. The differences of other empirical orbit parameters indicate a non-systematic flat course with the median of 0.001-0.003 nm/s 2 , thus they are not shown here. Although the SLR residuals do not directly indicate a systematic effect for high |β| angles for the FOC satellites, the differences between B 0 parameters, estimated using the microwave (M1) and combined (W2) solutions, show a dependency on |β|. The overall median difference of B 0 for the Galileo-FOC E08 when |β|> 60° is at the level of 0.16 nm/s 2 , which corresponds to a reduction at the level of 25% toward zero value. The maximum difference of B 0 is 1.86 nm/s 2 for FOC. The majority of B 0 differences for Galileo-IOV E19 are negative with the maximum value of − 2.40 nm/s 2 and the median improvement is at the level of 0.50 nm/s 2 , which corresponds to the reduction of the estimated values of B 0 by 23% toward zero accelerations.
To conclude, for both IOV and FOC satellites, the absolute values of B 0 estimates in the combined solutions are lower, i.e., closer to 0. That is especially visible for high β angles when the direction B is close to being aligned with the orbit radial direction and would fully correspond to it if |β| reaches 90°. The radial direction is directly measured by SLR. Therefore, SLR stabilizes the orbit for the maximum values of the β angle, resulting in lower estimates for the B 0 parameter.

Impact of the SLR observations on the satellite positions
Let us now check the impact of the combined solution on the satellite positions. Although the differences between satellite positions neither specify the absolute orbit accuracy nor indicate a more suitable orbit determination strategy, they may show differences induced by, in this case, the addition of the SLR observation to GNSS satellites in the processing.
We compute the differences between satellite positions calculated from the microwave (M1) and combined (W2) solutions every 15 min. Figure 9 illustrates the differences for Galileo-FOC E08 and Galileo-IOV E19 as a function of the β angle. This analysis corresponds to the results from Fig. 8, i.e., the differences in the satellite positions in the radial direction are consistent with the differences between the estimated parameters B 0 . For both Galileo-FOC and IOV, the differences in the radial direction have the same sign as the differences between B 0 estimates. For the Galileo-FOC, the differences reach up to 3.9 cm with the STD of 1.3 cm for |β|> 60°. For the Galileo-IOV E19, for which the differences are negative, the highest absolute impact from the addition of SLR is at the level of 5.6 cm. Differences are also visible for the eclipsing period for the along-track component, i.e., for |β|< 12.6°, they fall between − 3.4 and 2.2 cm for E08 and − 4.3 and 2.1 cm for E19. The cross-track component does not show a prominent pattern, as it is not determined by SLR as reliably as the two remaining components . The cross-track component is well defined by a dense GNSS network, and hence, SLR has a very moderate impact on the orbit position in this direction, even when employing the highest weights.

Impact of the combined solution on the Length-of-Day parameter
The orbit solution cannot be considered stable unless the remaining estimated global geodetic parameters are determined reliably. The quality of the estimated ERPs is crucial for POD as they comprise the transformation parameters linking the celestial reference frame with the TRF. The rate of UT1-UTC is sensitive to systematic errors included in satellite orbits, thus can be used as another orbit validation tool. The quality of ERPs from multi-GNSS solutions is studied in Zajdel et al. (2020); hence in this study, we focus on the accumulated LoD values as this parameter depends on Kepler orbit parameters, that is, d dt (UT1 − UTC) = −LoD = −̇Ω +cos(i)•̇u k (Rothacher et al. , 2001Ray 2009), where Ω is the right ascension of the ascending node, i is the inclination, u is the argument of latitude of the satellite, and k is the ratio of the universal time to sidereal time. The ERPs are estimated with 24 h resolution as two offsets at the beginning and the end of each day. The first value of UT1-UTC for the midnight epoch is fixed to the a priori value from IERS-C04-14 provided by Very Long Baseline Interferometry (VLBI). LoD is then estimated as the change of UT1-UTC over one day. Figure 10 illustrates the accumulated LoD for all strategies. We check how quickly LoD would drift away from the reference value. The solution which is based solely on the GNSS data indicates the highest drift at the level of − 4.8 ms/year. When adding SLR observations, the drift diminishes. However, when the contribution of SLR has the lowest weight in W1, Fig. 9 Differences between satellite positions determined using the microwave (M1) and combined (W2) solutions decomposed into the radial, along-track, and cross-track directions for Galileo-FOC E08 (left) and Galileo-IOV E19 (right) Fig. 10 Accumulated LoD values with respect to the IERS C04 for all the test cases with fitted linear regression function the drift is at the level of − 4.5 ms/year. On the other hand, when the SLR is treated equally with the GNSS data in W3, the drift diminishes to the level of − 0.9 ms/year. However, according to all the analyses, the strategy W2 provides the best results in terms of the orbit quality. The drift of LoD for W2 is at the level of − 3.6 ms/year which is better by 20% when compared to the microwave solution.

Conclusions
All the Galileo satellites are equipped with LRAs for laser ranging. Thanks to observations provided by the ILRS stations, it is possible to conduct the combined SLR and GNSS orbit determination with an appropriate weighting between the two techniques. In this study, we determined Galileo orbits using combined GNSS and SLR observations and checked the impact of the SLR technique on the Galileo POD by testing different weighting strategies.
We found that the higher the weight of the SLR observations, the greater the reduction of the formal error of the semimajor axis that equals up to 62% for the solution W3 when adding 50 SLR observations. We also found the reduction of the dependency of the formal error of the semimajor axis on the Sun elevation above the orbital plane. The largest improvement is for satellites with the highest β angles when the Sun is almost perpendicular to the orbital plane, and the lowest β angles during the eclipsing periods. However, the orbit misclosure statistics for the combined solution for strategy W3 with equal weights is worse than for the GNSSonly solution and for strategies W1 and W2. For strategy W3 the 3D STD is 131 mm for FOC while it is at a level of 116, 106, and 117 mm for M1, W1 and W2, respectively. For W3, the STD of SLR residuals is highest from all combined solutions, i.e., 26.5 mm for FOC satellites, despite the mean offsets equal to zero, which suggests that this strategy is dominated by the SLR technique. The lowest STD of SLR residuals was obtained for solution W2, i.e., 24.4 mm as compared to 26.6 mm for the GNSS-only solution. Moreover, for the strategy W2, we noticed significant mitigation of the systematic effects for the Galileo-IOV satellites for |β|> 60°, i.e., the STD of SLR residuals for Galileo-IOV for high |β| is at a level of 30.8 and 29.6 mm for W1 and W2, respectively as compared to 36.3 mm for M1.
The differences between the satellite positions calculated using solutions M1 and W2 indicate that the highest differences occur for the radial component, which is directly measured by SLR. The largest differences are at the level of 56 mm for the Galileo-IOV E19 for |β|> 60°. A similar dependency was found for the differences between the estimated empirical orbit parameter B 0 from solution M1 and W2. For both IOV and FOC satellites, the estimated B 0 parameter is lower for the combined solution when |β|> 60° which indicates that the SLR observations improve the Galileo orbits for high |β|-angle periods. Additionally, the drift of the accumulated LoD parameter is diminished by 20% for the solution W2 when compared to the solution M1.
The quality of the combined solution is conditioned by the limitation factors affecting GNSS and SLR techniques. Concerning the GNSS technique, for Galileo, we use phase center offsets (PCO) calculated by Steigenberger et al. (2016) in the frame of Multi-GNSS Experiment (MGEX) as these values are consistent with ITRF2014 and IGS14 scales. The Galileo chamber-calibrated values of PCO and phase center variations (PCV) from the Galileo metadata might be considered to increase the quality of the GNSS solutions. However, it also introduces an offset with respect to the ITRF2014 scale. For the SLR technique, one should analyze in detail the SLR signature effect, which can introduce a systematic effect at the level of 15 mm for different SLR stations (Sośnica et al. 2018) as well as the scale issues in the SLR network related to missing RBs for some stations in the ITRF2014 realization. Further inconsistencies to be solved are the result of GNSS orbit modeling deficiencies, requiring, e.g., proper modeling of thermal effects related to the variable temperature of different satellite surfaces.
To conclude, the impact of the SLR observations in the combined solution may not be very large due to the limited number of SLR data provided by thirteen stations at maximum. The small number of observations enforces the incorporation of greater weights for SLR. A significant increase of the number of SLR observations would be most profitable for the combination because it would naturally increase the influence of the SLR in the combined solution. A higher number of SLR observations would also allow for the application of more restrictive data screening criteria as well as for dividing the SLR observations into two independent sets (1) dedicated to POD and (2) employed for the solution validation.
Currently, the combined GNSS and SLR solutions using weighting schemes W1 and W2 provide accurate Galileo orbits with significantly diminished systematic effects for |β|> 60° owing to the stabilization of the satellite positions by the SLR observations. More SLR observations of GNSS satellites would be beneficial for the identification of remaining inconsistencies and, eventually, for the full exploitation of the SLR-GNSS co-location in space. PI/PRO/2019/1/00004/U/00001. We also thank the Wroclaw Center of Networking and Supercomputing computational grant using MATLAB Software License No: 101979 (https ://www.wcss.wroc.pl).

Data availability
The GNSS and SLR observations used for this study can be freely downloaded from the Crustal Dynamics Data Information System ftp://cddis .nasa.gov/gnss/ and ftp://cddis .nasa.gov/slr/produ cts. The metadata for the Galileo satellites are provided by the GSA at https ://www.gsc-europ a.eu/.
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/.
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://creativecommons.org/licenses/by/4.0/.