GLONASS precise orbit determination with identification of malfunctioning spacecraft

Due to the continued development of the GLONASS satellites, precise orbit determination (POD) still poses a series of challenges. This study examines the impact of introducing the analytical tube-wing model for GLONASS-M and the box-wing model for GLONASS-K in a series of hybrid POD strategies that consider both the analytical model and a series of empirical parameters. We assess the perturbing accelerations acting on GLONASS spacecraft based on the analytical model. All GLONASS satellites are equipped with laser retroreflectors for satellite laser ranging (SLR). We apply the SLR observations for the GLONASS POD in a series of GNSS + SLR combined solutions. The application of the box-wing model significantly improves GLONASS orbits, especially for GLONASS-K, reducing the STD of SLR residuals from 92.6 to 27.6 mm. Although the metadata for all GLONASS-M satellites reveal similar construction characteristics, we found differences in empirical accelerations and SLR offsets not only between GLONASS-M and GLONASS-M+ but also within the GLONASS-M+ series. Moreover, we identify satellites with inferior orbit solutions and check if we can improve them using the analytical model and SLR observations. For GLONASS-M SVN730, the STD of the SLR residuals for orbits determined using the empirical solution is 48.7 mm. The STD diminishes to 41.2 and 37.8 mm when introducing the tube-wing model and SLR observations, respectively. As a result, both the application of the SLR observations and the analytical model significantly improve the orbit solution as well as reduce systematic errors affecting orbits of GLONASS satellites.


Introduction
GLONASS observations, together with GPS, are the only GNSS data so far used to generate the official products of the International GNSS Service (IGS, Johnston et al. 2017). Among the IGS products, we can distinguish, e.g., precise orbits of GNSS satellites and the earth rotation parameters (ERP); therefore, GLONASS observations contribute to the precise positioning and the description of the earth system. As a result, it is crucial to properly handle the precise GLO-NASS orbit modeling, especially direct solar radiation pressure (SRP), which is the greatest non-conservative perturbing force acting on the GNSS satellites (Milani et al. 1987). A considerable advantage of the GLONASS satellites is the fact that all spacecraft are equipped with laser retroreflector arrays (LRA) for satellite laser ranging (SLR), the precise space geodetic technique that significantly contributes to the GNSS POD .
Currently, the GLONASS system consists of satellites of type M, K, and M+ spacecraft. The M+ satellites are the updated version of the GLONASS-M spacecraft that broadcast signals at three frequencies, including L3, and are capable of laser-based time transfer. The M+ satellites carry devices for the Inter Satellite Laser Navigation and Communication System, allowing for communication between spacecraft using a laser. The latter introduces challenges for POD as the inter-satellite link devices change surface properties and the shape of the satellite bus. The most pronounced differences in terms of POD occur between GLO-NASS-M and GLONASS-K satellites. First of all, the two types of satellites are of different shapes, i.e., GLONASS-M satellites are elongated and cylindrical, whereas type K satellites are cuboids (see Fig. 1). Moreover, the surface of the solar panels of the M-type satellites is vast when compared 1 3 36 Page 2 of 13 to the K-type satellites, i.e., 30.85 m 2 as compared to 16.96 m 2 . GLONASS-K are significantly lighter than GLONASS-M, i.e., 996 kg vs. 1415 kg. LRA on each type of satellite is different as well. GLONASS-M satellites are equipped with a rectangular LRA consisting of 112 corner cubes, whereas K-type satellites carry ring arrays mounted around the microwave antenna that consist of 123 corner cubes Bury et al. 2019a, b). Apart from geometrical diversities, GLONASS-K satellites can broadcast signals on third, L3 frequency, perform the time transfer, and are characterized by extended lifetime (Zaminpardaz et al. 2017).
The modernization of GLONASS introduces a need for continuous development of the POD strategy. Due to the elongated cylindrical shape of the GLONASS-M satellites, the classical empirical CODE orbit model (ECOM, Springer et al. 1999) was insufficient for modeling GLO-NASS orbits. Therefore, an extended ECOM2  model was introduced to cope with SRP acting on GLONASS spacecraft. However, Prange et al. (2017) indicated a list of satellites for which the new ECOM2 was insufficient, i.e., satellite vehicle number (SVN): 721, 723, 725, 730, 736, 737. Dach et al. (2019) re-estimated the GLONASS satellite antenna offsets, reducing the SLR residuals to orbits of affected satellites. Despite advances in GLONASS POD, their orbits cannot reach the accuracy of GPS or the Galileo system. As a result, the limited accuracy of the GLONASS orbits prevents applications related to a reliable description of the earth system, i.e., determination of high-quality ERPs or geocenter coordinates . One of the greatest limitations for the GLONASS POD is the lack of the accurate optical and geometrical properties of GLONASS satellites. Based on the metadata, it is possible to develop an analytical boxwing model for cuboidal satellites or a tube-wing model for cylindrical satellites to absorb SRP and earth radiation pressure (Rodriguez-Solano et al. 2012a). Such a box-wing model significantly improves the orbits of Galileo satellites for which the metadata are officially published Duan et al. 2019;Li et al. 2019). For GLONASS satellites, approximated values provided by Rodriguez-Solano (2014) are available with the adjusted optical and geometrical parameters derived by Duan et al. (2020). In 2021, IGS published results of the 3rd reprocessing campaign (repro3, Rebischung 2021) of the reanalysis of the GNSS historical observations collected by the IGS network since 1994 for the International Terrestrial Reference Frame (ITRF2020) realization. The recommendation for the SRP modeling was to use the empirical ECOM or JPL GSPM models, or an analytical model based on the properties released by IGS (http:// acc. igs. org/ repro3/ repro3. html). Figure 2 shows that the orbital artifacts are still present for some GLONASS, especially after 2018 when the SLR residuals to the experimental IGS repro3 GLONASS orbits significantly increase. Figure 2 reveals that the GLONASS POD needs further improvements to reach the same quality as currently achieved for the IGS GPS and Galileo combined orbits .
Stations of the International Laser Ranging Service (ILRS, Pearlman et al. 2019) conduct laser ranging to any spacecraft equipped with LRA. To track the satellites, stations need all the necessary data, including orbit predictions. SLR-to-GNSS observations comprise around 10% of all observations collected by the ILRS . SLR data allows for validating the GNSS-based orbit products and may contribute to the GNSS POD (Zajdel et al. 2017;Bury et al. 2021).
The goal of this study is to revise the state-of-the-art strategies for GLONASS POD. We check if the analytical model developed for the GLONASS satellites using unofficial optical properties of the satellite surface elements is beneficial for POD. We also assess the impact of the addition of the SLR-to-GNSS observations on the GLONASS orbit quality. Finally, we try to identify GLONASS satellites for which the orbit solution is inferior, i.e., with the accuracy significantly lower than the IGS declared accuracy of 3 cm. In this study, they are called 'malfunctioning satellites' although the reason for poor orbit quality is usually unknown. The satellites may have the 'healthy' status and can be used for the navigation purposes, but in the precise geodetic applications, the inferior orbit quality prevents sub-centimeter solutions. We

Methodology
We conducted a series of GLONASS POD strategies over the period 2017-2018. We focus on two aspects (1) handling of SRP using a hybrid approach which considers an analytical model with a set of empirical parameters, and (2) the impact of the addition of the SLR observations on the GLO-NASS POD solutions. All the calculations are conducted in the Bernese GNSS Software 5.2 .
We focus on a hybrid SRP modeling, as many studies indicated its superiority against the empirical solutions for Galileo, QZSS, and BeiDou-3 (Bury et al. 2019a, b;Li et al. 2019Li et al. , 2020, especially in the light of a series of GLONASS satellites for which the empirical models are insufficient. We test two approaches H1 and H2, which consider the analytical a priori model and set of ECOM and ECOM2 parameters, respectively. The assumptions of the tube-wing or the box-wing model are consistent with Bury et al. (2019a, b), considering shadows resulting from the earth and moon as a fraction of the eclipsed part of the sun disk comprising antumbra and penumbra periods. The tube-wing model considers the cylindrical shape of the GLONASS-M buses, whereas the box-wing model for GLONASS-K considers the flat surfaces of the cuboidal buses and the optical properties of satellite elements. Both box-wing and tube-wing models treat separately the solar panels, which are described with the equation: and the bus of satellites. Apart from the flat surfaces, which are described by the formula: the GLONASS-M bus also contains the cylindrical surfaces for which the accelerations are described by the equation: In (1-3) A denotes the surface area, m is the mass of the satellite, S c denotes the solar constant (1367 W/m 2 ) rescaled by the AU r s 2 , where AU is the astronomical unit and r s is the instantaneous distance of the satellites from the sun, c is the speed of light. An angle between the unit vector of the surface normal e n and the unit vector of the direction of the (1) illuminating source e ⊙ is described by . All the thermal effects result from an impulse transfer of the absorbed and emitted photons acting on the satellite's surface illuminated either by the sun or the earth's radiation pressure with no lag considered. All the energy is transferred with the assumption of the photon balance, i.e., α + δ + ρ = 1, where α, δ, ρ describe absorbed, diffusely reflected, and specularly reflected photons, respectively. All satellite surface properties are consistent with the recommendation for repro3 processing and were calculated by Rodriguez-Solano (2014). Additionally, to the hybrid test cases, we compute the solution E2 in which we use ECOM2 without applying the a priori model. For the attitude model, we assume the yawsteering mode with the ECOM parameters set to 0 when a satellite enters the earth's shadow. The methodology of the combination of GNSS and SLRto-GNSS observations is consistent with Bury et al. (2021) with the modifications of the weighting strategies for GLO-NASS satellites. The sigma of GNSS phase observations is assumed to be at the level of 1.5 mm. Following Prange et al. (2017), we pre-selected GLONASS satellites whose orbits are of poor quality (GLONASS-BAD, GB), i.e., SVN: 719 (R20), 730 (R01), 733 (R06), 734 (R05-until August 2018), 736 (R16). For GB satellites, the standard deviation (STD) of SLR residuals calculated using observations delivered by the best-performing stations (Zajdel et al. 2019) was greater than 40 mm. For the remaining GLO-NASS (GLONASS-GOOD, GG), the STD of SLR residuals was at the level of 20-26 mm. We used the STD of the SLR residuals as weighting factors in strategy C1A for the two sets of GLONASS spacecraft (see Table 1). We also computed the solution C1C, which treats the SLR-to-GLONASS observations with sigma better than in C1A by a factor of ten, and the solution C1B, which comprises an intermediate case between C1A and C1C. The sigma values for GG and GB satellites in particular solutions are collected in Table 1.
In all the scenarios, the geodetic datum is defined using minimum constraints: no-net-translation (NNT) and no-net-rotation (NNR) imposed on a series of stable stations that are referred to IGb14 (Rebischung et al. 2016). The GNSS processing scheme is based on double difference network processing with the ionospheric-free linear combination. The SLR a priori coordinates come from the SLRF2014 reference frame, which is the ILRS realization of the ITRF2014 (Luceri et al. 2019). The GNSS network is realized consistently with Zajdel et al. (2019), where the set of stable GNSS sites is defined using an iterative Helmert transformation. For the combined GNSS and SLR solution, we treat the GNSS network as in the microwave-only solutions, whereas for the SLR network, we apply NNT and NNR constraints to the set of core SLR stations.
For the POD, we use GNSS observations sampled by 180 s, and for the combined GNSS and SLR solutions, we use SLR normal points with 300 s sampling. Phase center offsets and variations for the GLONASS satellites are consistent with IGS aligned to the ITRF2014 scale. In all solutions, we apply albedo and earth infrared modeling (Rodriguez-Solano et al. 2012b) as recommended for repro3. We assume 100 W for all GLONASS satellites for the navigation antenna thrust to avoid scale inconsistencies with respect to ITRF2014, instead of using the satellite-specific values from ground-based measurements (Steigenberger et al. 2018). We estimated the one-day orbital arcs, with the 5-min orbit integration. Additionally, we estimated stochastic orbit parameters in radial, along-track, and cross-track directions every 12 h. Apart from the orbital parameters, in every test case, we estimated ERPs, geocenter coordinates, GNSS and SLR station coordinates, and implicitly the troposphere.

Accelerations acting on GLONASS satellites
The analytical model allows for the assessment of the perturbing accelerations acting on GNSS satellites. This section analyzes accelerations resulting from SRP and earth radiation pressure.
Due to the greatest impact of the SRP-induced accelerations, the perturbing forces are decomposed into directions: D-pointing from the sun to the satellite, Y-parallel to the solar panel rotation axis, and B-completing the righthanded orthogonal frame . Perturbing accelerations also depend on the elevation of the sun above the orbital plane-the β angle, and the position of the satellite w.r.t earth, i.e., the argument of latitude of the satellite w.r.t. the argument of latitude of the sun-Δu. When describing the satellite reference frame, we use the nomenclature consistent with IGS which is described by Montenbruck et al. (2015). Figure 3 depicts the cumulative impact of the accelerations induced by SRP and earth radiation pressure acting on the GLONASS-M calculated based on the a priori analytical model. The highest accelerations are noted for direction D.
The acceleration differences for positive and negative β result from the earth-sun distance during one full cycle of the β angle which is at the level of 3.2% between aphelion and perihelion. Figure 4 illustrates the accelerations for different values of β angle as a function of Δu. The highest accelerations are visible in direction D and reach -150 nm/s 2 for the highest absolute values of the β angle. The course of the accelerations is almost flat as compared to the lower values of β angles. This results from the yaw-steering mode in which the satellite maneuvers the solar panels to track the sun, simultaneously pointing the navigation antenna toward the earth. When the sun is at the highest elevation, it illuminates the largest area of the GLONASS satellites, i.e., solar panels and the elongated part of the bus.
The flat course of the accelerations for the highest β angles is also noted for direction B. However, these accelerations result mainly from earth radiation pressure. For lower β angles, accelerations in D and B become dependent on the Δu due to the yaw-steering attitude in which satellites maneuver their buses and the area of illuminated surface changes. Note that the accelerations in the direction B are smaller by two orders of magnitude when compared to the accelerations in D. For D direction, the lowest accelerations occur for low β angles when Δu is close to 0°, and the satellite is in the line between the sun and earth when the illuminated area of the bus is the smallest. In the earth shadow for |β| < 14.2° and Δu close to 180°, the SRP accelerations decrease toward the 0 value (Dilssner et al. 2011). The residual perturbing forces result from earth radiation pressure, which is not higher than 1 nm/s 2 . The yaw-steering attitude assumes that the solar panels track the sun direction. If neither lags nor bias are present in the attitude of the solar panels, the Y-accelerations should be close to the 0 value. For GLONASS satellites, the accelerations in Y are by two orders of magnitude lower than those in direction B. Figure 5 depicts the amplitude spectra of the accelerations illustrated in Fig. 4. Based on the Fast Fourier Transform, it is possible to assess the amplitudes of perturbing accelerations and, on their basis, calculate an error of the satellite position when neglecting the particular acceleration (Table 2). We evaluated the position error for typical harmonics of the GLONASS satellite revolution period and checked to what extent the empirical ECOM2 model absorbs them. The satellite position error is calculated using Eq. (3) from Bury et al. (2020). Apart from the constant accelerations in each D, Y, B direction, the ECOM2 model considers even terms in direction D and odd terms in B. The highest amplitudes of accelerations occur twice-perrevolution and cause an error in the satellite position at the level of 55 mm. However, twice-per-revolution terms are considered in ECOM2; thus, accelerations of this type do not affect the orbit solution when using ECOM2. Once-perrevolution accelerations in D directions emerge from differences between + Z and -Z satellite panels, i.e., the front panel with navigation antennas and the back panel pointing toward outer space. The once-per-revolution accelerations in D are neglected in ECOM2 despite that ECOM2 was developed for GLONASS ; however, they were estimated in the classical ECOM by some IGS   analysis centers in the so-called 9-parameter ECOM model. For the neglected in ECOM2 once-per-revolution D term, the amplitude of acceleration is at the level of 1.35 nm/s 2 for low and moderate β angles, which translates into the satellite position error at the level of 56 mm. The resulting accelerations suggest that for GLONASS, the once-per-revolution and the twice-per-revolution accelerations in the D direction must be estimated in pure empirical models without any a priori models. Such an approach is neither applied in the classical ECOM nor ECOM2. For the IGS repro3, some IGS analysis centers employ the empirical ECOM2, which omits the once-per-revolution accelerations in D and results in substantial orbit errors. In terms of direction B, the once-per-revolution accelerations of 2.44 nm/s 2 could cause the periodical position error of 102 mm. The errors emerging from the omission of periodical terms in direction Y do not exceed 0.4 mm, which confirms no need to estimate periodical terms in Y. The analytical model absorbs all the calculated position errors. The residual accelerations can be absorbed by the empirical parameters estimated in the hybrid POD strategy described in the section 'Results.'

Results
In this section, we analyze the results of the GLONASS POD. We evaluate the precision and accuracy of the orbit solution and check the impact of the introduction of the a priori analytical model and SLR-to-GNSS observations for POD.

Orbit misclosures
The orbit misclosures describe the differences between satellite positions from the consecutive orbital arcs. The day boundaries are the least stable part of the arc; thus, they indicate weaknesses of the solution. During the two years of the analysis, only eight out of 24 satellites were tracked intensively by SLR sites, i.e., more than 10,000 observations were collected to these satellites. Therefore, taking into consideration the number of SLR observations, for further analyses, we will focus on solution C1B, which increases the impact of SLR data as compared to C1A. Based on the orbit misclosures for the solution H1, we identify GLONASS satellites whose orbits are deteriorated and confront our list of the GB satellites. Figure 7 illustrates the time series for the GLONASS orbit misclosures in 2017. The preliminarily chosen list of GLONASS satellites based on the SLR residuals to CODE orbits contains satellites of PRNs (SVNs): R01(730), R05(734), R06(733), R16(736), and R20(719). According to the figure, all the satellites from the GB list are characterized by higher orbit misclosures. Moreover, based on Fig. 7, we can identify additional deteriorated orbits for satellites R14(715), R24(735), and R26(801). On average, the median value of the 3D misclosures for the malfunctioning satellites is at the level of 138 mm, as compared to 82 mm for the remaining satellites. In the following analyses, we investigate if the addition of the SLR-to-GNSS observations improves the quality of the orbits, especially those of poor quality.

Empirical orbit parameters
Empirical orbit parameters account for unmodeled perturbing forces acting on GNSS satellites. When no optical and geometrical data describing satellites are available, the empirical models are the most effective tool for the absorption of the direct SRP. When using the analytical-only The box extends from Q1 to Q3 quartile (the interquartile range, IQR), orange lines denote the median, and whiskers extend to the last data point in range covering no more than 1.5 × IQR approach for the SRP handling, the orbit solution is deteriorated due to a lack of D 0 parameter, which absorbs a changeable acceleration induced by direct SRP (Bury et al. 2019a. Therefore, the empirical parameters used together with the analytical model absorb the residual accelerations. Based on the residual accelerations, it is possible to assess the effectiveness of the tube-wing and box-wing models. Figure 8 depicts the empirical parameters estimated for GLONASS-M R07 in solutions E2 and H2, indicating only the impact of applying the tube-wing model. The substantial impact of the tube-wing model is noted for the D 0 parameter that absorbs accelerations acting on the solar panels and the cross section of the satellite bus.  On average, the application of the tube-wing model reduces the D 0 parameter from -147.88 nm/s 2 to 0.26 nm/s 2 for the GLONASS-M satellites. For the GLONASS-K satellites, in the analyzed period, the constellation consists of two spacecraft, R09 (802) and R26 (801), which are treated in the same manner in terms of the box-wing model, although the IGS-MGEX metadata divide GLONASS-K into GLO-K1A (SVN801) which carries an additional antenna for L3 signal and GLO-K1B, i.e., SVN802 and SVN805 launched in late 2020, not included in this study (https:// www. igs. org/ mgex/ metad ata/# metad ata-sinex-format). The D 0 parameter indicates differences between two GLONASS-K spacecraft; in solution E2, the D 0 parameter is at the level of -110.0 and -104.2 nm/s 2 for R26 and R09, respectively. When the analytical box-wing is applied, the D 0 parameter indicates positive values at the level of 2.2 and 8.5 nm/s 2 for R26 and R09, respectively (Fig. 9). The positive values of the estimated D 0 in solution H2 indicate that either optical or geometrical parameters are overestimated or the mass of the satellites is poorly approximated.
The cosine terms of the ECOM2 model, i.e., D 2C and B 1C parameters, are mitigated as they absorb the variable SRP acting on the rotating bus in the yaw-steering attitude mode. An offset (and STD) is diminished by 1.7 (0.5) and 1.6 (1.2) nm/s 2 for D 2C and B 1C , respectively. The parameters Y 0 and B 0 do not significantly change after the addition of the tubewing model. Y 0 is connected with a possible solar panel rotation lag or thermal energy emission. B 0 is not as sensitive to a priori model for GLONASS as for the Galileo satellites due to the specific shape and higher masses of GLONASS-M. Finally, the sine terms D 2S and B 1S are barely different in the two solutions. The sine terms do not have interpretation connected with SRP and absorb thermal re-radiation and spacecraft misorientation effects . Figure 9 shows the residual D 0 parameter estimated using strategy H1 for each GLONASS orbital plane. The more reliable parameters describing the box-wing model, the closer the estimated D 0 parameter to zero. The negative D 0 values mean that the analytical surface parameters are underestimated, and positive values indicate that the parameters are overestimated (or the satellite mass is underestimated). First of all, we notice the differences between the two GLONASS-K spacecraft which are treated in the same manner. The differences between the D 0 parameters describing orbits thereof indicate significant differences between two GLONASS-K satellites. The launches of the two satellites were separated by three years; therefore, the two spacecraft might be different in terms of materials covering particular surfaces or the payload which transfers into different masses. For the GLO-NASS-K, a change of mass at the level of 70 kg refers to a change of the D 0 at the level of 10 nm/s 2 . Thus, an increase in the nominal mass by 60 kg for the GLONASS-K R09 in the a priori model would reduce the D 0 to near-zero values.
In this study, we treat all the GLONASS-M consistently. The residual accelerations absorbed by D 0 for GLONASS-M satellites from planes R-I and R-II oscillate around zero with the increased scatter during the eclipsing period (see Fig. 9). In the plane R-I, we identified the GLONASS-M+ , which started broadcasting a signal at R05 in August 2018. The most pronounced differences are noted for plane R-III for which the D 0 parameter is different for each satellite. However, the largest anomaly occurs for the satellite R23 (SVN732), for which the D 0 parameter is out-of-phase compared to the remaining satellites. Moreover, for the empirical solution E2 (not shown here), the D 0 parameter behaves similarly to other M-type satellites. As a result, there must have been an attitude change of the R23 satellite as the box-wing model for this spacecraft is not suitable as for the remaining M-type satellites. This anomaly has no impact on the POD result, as it is compensated for the D 0 parameter.
The differences between the estimated D 0 parameters result from the equal treatment of all GLONASS-M satellites. The oldest GLONASS-M spacecraft are already 12 years old during the analysis period, whereas the designed GLONASS-M lifetime was 7 years. As a result, the applied optical parameters might become obsolete, and during such a long period, it is inevitable that the construction technology of GLONASS satellites must have been updated; therefore, the mass of the satellites may be different from the assumed value. The differences in the estimated D 0 parameter highlight the necessity for the mission providers to release accurate metadata which consider changes in the optical properties of the satellites and their masses. Nevertheless, Zajdel et al. (2020) indicated that even the a priori model composed using the approximated values of the GLONASS satellite properties helps to reduce spurious signals in the Z-component of the geocenter coordinates and other geodetic parameters.

SLR residuals analysis
SLR residuals are a valuable tool for the quality assessment of microwave-based POD products. However, for the combined SLR and GNSS solution, SLR observations are no longer independent as the same set of the SLR data contributes to the POD and the validation. Nevertheless, Bury et al. (2021) indicated that the laser-ranging data, i.e., the SLR residuals analysis, can still serve as a POD quality assessment indicator when using a reasonable weighting of SLR observations. Figure 10 illustrates SLR residuals for solutions E2, H1, and C1B, as a function of β and Δu angles. For all GLO-NASS, the highest SLR residuals are obtained for solution E2. For GLONASS-M and -M+ , systematic effects are present during orbit midnight; however, for -M+ satellites, the offset of SLR residuals is positive and exceeds 40 mm, whereas for M-type satellites, the effect is negative. For the K-type satellites, the solution E2 is insufficient as the STD of SLR residuals is at the level of 93 mm. The advantage of solution H1 over solution H2 is visible in the SLR residual statistics. Despite reducing the STD of the SLR residuals, solution H2 introduces an offset to the SLR residuals compared to H1 (see Table 3). The introduction of the box(tube)-wing models between solutions E2 and H1 diminished the STD of SLR residuals by 65.0 and 3.0 mm for K-type, and M-type satellites, respectively. The reduction in STD of the SLR residuals for K-type satellites revealed systematic effects for the orbit noon and midnight periods.
The introduction of the tube-wing model for GLONASS-M+ diminished the offset of SLR residuals by 5.2 mm, but the STD remains at the same level as for solution E2. A similar pattern for the M+ satellites was found by Sośnica et al. (2020), who analyzed the quality of the experimental combined orbits of the multi-GNSS Pilot Project (MGEX, Montenbruck et al. 2017). For all GLONASS, the introduction of the tube-and box-wing model shifts the mean values of the SLR residuals toward the zero value (see Table 3 and Fig. 11). The largest impact of the a priori model is noted for the GLONASS-K, but the improvement of the SLR residuals is visible also for M-type and M+ satellites. For example,   Fig. 11 shows that the box-wing model diminished the IQR of the SLR residuals from 45.3 to 33.0 mm and from 36.6 to 33.7 mm for M-type and M+ satellites, respectively. The addition of the SLR observations has no significant impact on the orbit quality compared to the introduction of the box(tube)-wing model. Nevertheless, C1A and C1B cases contribute to mitigating the SLR residuals for all the satellite types, whereas solution C1C indicates a significant advantage over the other combined solution for GLONASS-M+ . Noteworthy, the STD of the SLR residuals decreased below the level of 30 mm for all combined SLR + GNSS strategies for GLONASS-M satellites.
Let us now check the impact of the a priori model and the combined GNSS-SLR solutions on the SLR residuals for each GLONASS satellite. First of all, we focus on the GLO-NASS SVN730, whose SLR residuals are discussed in Fig. 2. When calculating the POD solution using strategy E2, we obtain a similar pattern of the SLR residuals as for the repro3 orbit solution, i.e., the scatter of SLR residuals is increased after 2018. However, when introducing the tube-wing model, the spurious effect diminishes, i.e., the STD of the SLR residuals decreases from 48.7 to 41.2 mm (see Fig. 12). A further improvement is caused by introducing SLR observations-the STD of SLR residuals is mitigated to 37.8 mm. As a result, the reduction in the number of empirical orbit parameters and the a priori tube-wing model absorbed a possible anomalous attitude change of the satellite, which could not be compensated for by the ECOM2 model, whereas introducing SLR observations additionally stabilized the orbit solution.
In Fig. 13, we collect the SLR residuals for all GLONASS satellites considered in this study. We recognize the satellites by their SVNs due to changes of the GLONASS PRNs. Satellites from SVN715-854, excluding spacecraft 801 and 802, are the -M satellites. The SVN801-802 belong to K-type satellites and SVN855-856 are the modified M+ spacecraft. Although all the GLONASS-M and M+ are treated consistently, based on the SLR residual analysis, we distinguish several groups of spacecraft, as the period between the launch of SVN715 and SVN856 reaches 12 years. As a result, the changes in technology used for the construction of the GLONASS-M are highly possible. The optical properties for the oldest satellites became obsolete due to the material aging resulting from a continuous impact of the SRP acting on the satellites. For satellites SVN715-731, we see negative values of the SLR residuals in solution E2. For SVN732-736 in E2, the shift of the residuals is positive, suggesting differences in technology for the two groups of satellites. SVN715-731 and SVN732-736 satellites  The nomenclature of the boxes is consistent with Fig. 6. SVN naming is consistent with that used in the Bernese GNSS Software (dif-ferent GLONASS SVN naming can be found in the literature, e.g., GLONASS-K 801 can also be denoted as 701, GLONASS 851 can be denoted as 751, etc.) were launched between 2006 and 2010. Among the satellites SVN715-736, we identify those of poor orbit quality, i.e., 719 (R20), 730 (R01), 733 (R06), 734 (R05-until August 2018), 736 (R16). After adding the a priori orbit model and the SLR observations, both the IQR and median value of the SLR residuals are shifted toward the zero value, and the differences between the two groups are diminished for satellites with unreliable orbits. Spacecraft launched between 2011 and 2013, and 2014 and 2016, i.e., SVN742-747 and SVN851-854, respectively, are characterized by the orbits of better quality in the solution E2. However, the introduction of the tube-wing model and SLR observations improve the orbits for all the satellites, shifting the offset of SLR residuals toward the zero value and diminishing their STD. As for the GLONASS-M+ , i.e., SVN855-856, and K, i.e., SVN801-802, the addition of the tube-wing and the box-wing model and SLR data mitigates the SLR residuals for both types of satellites.

Quality of the orbit predictions
The stability of the orbit predictions corresponds to the quality of the orbit solution. In this study, we extrapolate a 2-day orbital arc and compare the second day of the orbit prediction with the corresponding day of the calculated orbit. Figure 14 illustrates the STD of differences between calculated and predicted orbit position decomposed into the radial, along-track, and cross-track directions.
We check if the addition of the SLR observations helps to stabilize the orbit prediction, especially for the satellites with poor orbit quality. The impact of the addition of the SLR observations is similar for GG and GB spacecraft. We can clearly distinguish the malfunctioning satellites (R01, R05, R14, R16, R20, and R24), as the quality of the orbit prediction deteriorates comparing to GG satellites. The impact of the SLR observations is marginal as the median STD of the orbit prediction calculated using solution C1B is improved only by 0.5, 14.7, and 1.1 mm in the radial, along-track, and cross-track directions, respectively, when compared to the microwave H1 solution.

Conclusions
GLONASS orbit determination remains challenging in several aspects, as demonstrated for the experimental repro3 orbits. In this study, we checked the impact of applying the analytical tube-wing and box-wing models as well as the SLR-to-GNSS observations on the GLONASS POD. We tested a series of hybrid POD strategies and found that the best hybrid solution contributing to GLONASS POD is based on the analytical box(tube)-wing model with the estimation of the ECOM 5-parameters (H1).
The quality of the developed analytical models was assessed using the estimated empirical parameters, especially the D 0 parameter, which absorbs most of the direct SRP impact. The residual accelerations absorbed by the D 0 parameter indicate the imperfection of the developed model, i.e., its optical properties or mass of the satellite. The analysis of D 0 parameters indicated that we should not treat all the GLONASS-M spacecraft in the same manner because, over 12 years, the optical properties of surface materials or the technology of their payload must have changed. Therefore, the precise information on the optical and geometrical properties and the precise mass of the satellites is crucial for GNSS POD. Nevertheless, using unofficial parameters for the composition of the tube-wing or the box-wing model and the application of the hybrid strategy improves the GLONASS orbit quality, especially for the satellites characterized by deteriorated orbits identified in this study. Orbits of the satellites SVN: 719, 720, 730, 733, 734, 736 have been significantly improved when using the H1 POD strategy.
The further improvement of the GLONASS orbit quality occurred when we introduced SLR-to-GNSS observations, which mitigated the systematic effects in the SLR residuals. For the SVN730, the STD of the SLR residuals has been diminished from 48.7 to 41.2 mm. The introduction of the SLR observations to POD reduced the STD of SLR residuals to the level of 37.8 mm. In general, the introduction of the box-wing model (and SLR observations) for POD diminished the SLR residuals from 93 to 28 (26) mm, and 34 to 31 (28) mm, for -K, and -M spacecraft, respectively. Only eight out of 24 satellites are tracked intensively by SLR sites, i.e., more than 10,000 observations were collected to these satellites during the two years of analysis. As a result, the impact of the combined solutions (C1A-C) on the orbit prediction is marginal when compared to solution H1. As a result, now, the impact of the SLR observations has to be increased by imposing higher weights on the SLR observations. However, an extensive increase in the SLR weights destabilizes the orbit solution, as in the case of C1C. This is why the increase in the number of SLR observations is crucial for the full exploitation of their potential for precise GNSS orbit determination.
Acknowledgements The IGS MGEX is acknowledged for providing multi-GNSS data and satellite metadata, and ILRS for organizing GNSS tracking campaigns. This work was funded by the National Science Centre (NCN), Poland. G. Bury is supported by the grant UMO-2018/29/N/ST10/00289. K. Sośnica, R. Zajdel, and D. Strugarek are supported by the grant UMO-2019/35/B/ST10/00515. G. Bury and R. Zajdel also acknowledge the Foundation for Polish Science (FNP) for the START scholarship.

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 GLONASS POD is obtained from the recommendation for the repro3 reanalysis http:// acc. igs. org/ repro3/ repro3. html.
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/.