Satellite laser ranging to GNSS-based Swarm orbits with handling of systematic errors

Satellite laser ranging (SLR) retroreflectors along with GNSS receivers are installed onboard numerous active low earth orbiters (LEOs) for the independent validation of GNSS-based precise orbit determination (POD) products. SLR validation results still contain many systematic errors that require special handling of various biases. For this purpose, we derive methods of reducing systematic effects affecting the SLR residuals to LEO Swarm satellites. We test solutions incorporating the estimation of range biases, station coordinate corrections, tropospheric biases, and horizontal gradients of the troposphere delays. When estimating range biases once per day, the standard deviation (STD) of Swarm-B SLR residuals is reduced from 10 to 8 mm for the group of high-performing SLR stations. The tropospheric biases estimated once per day, instead of range biases, further reduce the STD of residuals to the level of 6 mm. The systematic errors that manifest as dependencies of SLR residuals under different measurement conditions, e.g., elevation angle, are remarkably diminished. Furthermore, introducing troposphere biases allows for the comparison of the orbit quality between kinematic and reduced-dynamic orbits as the GPS-based orbit errors become more pronounced when SLR observations are freed from elevation-dependent errors. Applying tropospheric biases in SLR allows obtaining the consistency between the POD solution and SLR observations that are two times better than when neglecting to model of systematic effects and by 29% better when compared with solutions considering present methods of range bias handling.


Introduction
Numerous earth observation satellites belong to low earth orbiters (LEOs), which require precise orbit determination (POD) products to fulfill their mission requirements. POD products of LEOs are commonly derived from onboard Global Navigation Satellite System (GNSS; Johnston et al. 2017) receivers, primarily Global Positioning System (GPS), and additionally from Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS, Willis et al. 2010). The independent validation of GNSS-based LEO POD products is often conducted by the Satellite Laser Ranging (SLR) technique. Thus, LEO spacecraft with the demand of the highest orbit accuracy are equipped with GNSS receivers, often laser retroreflector arrays (LRA), and for some missions with DORIS receivers.
Swarm is the European Space Agency's (ESA) mission to provide a global representation of the geomagnetic field and its temporal evolution to improve the understanding of the earth's interior and the geospace environment of the sunearth system (Friis-Christensen et al. 2008). Three identical Swarm spacecraft, A, B, C, were launched on November 22, 2013. The initial altitudes of LEO satellites were about 480 km for Swarm-A and Swarm-C tandem and 510 km for Swarm-B (van den IJssel et al. 2020). All Swarm satellites are equipped with two dual-frequency GPS receivers (Zangerl et al. 2014) and a pyramid-shaped LRA (Neubert et al. 1998) for generating and validating POD products with accuracy at a cm-level.
Swarm POD products are generated using kinematic (KIN, Švehla and Rothacher 2005) and reduced-dynamic (RD, Wu et al. 1991) approaches. In general, KIN solutions represent a series of discrete satellite positions, whereas RD solutions represent a continuous orbit of a satellite integrated and subject to force models. ESA's Swarm Level 2 orbits are generated within the SCARF Consortium (Olsen et al. 2013) of research institutes (van den IJssel et al. 2015). Moreover, other universities and research institutions provide Swarm satellite POD products, e.g., the University of Bern (AIUB, Jäggi et al. 2016), the Graz University of Technology (Suesser-Rechberger et al. 2020), and the German Aerospace Center (Montenbruck et al. 2018).
Independent SLR orbit validation of LEOs equipped with LRA allows for the evaluation of GNSS-based POD products in terms of precision and accuracy. SLR validation of LEOs is based on satellite positions derived from POD products and observations provided by the International Laser Ranging Service (ILRS, Pearlman et al. 2019) stations and results in observation residuals. In the literature, the SLR validation of LEOs considers employing different sets of ILRS stations providing data, e.g., the high-performing stations, the ILRS core station list selected based on LAGEOS observations, or all possible contributing stations (Arnold et al. 2019a;Montenbruck et al. 2018;Strugarek et al. 2019). Moreover, different outlier rejection thresholds of SLR residuals are used, e.g., 20 and 25 cm as shown in Mao et al. (2021) and van den IJssel et al. (2015), respectively, as well as considering or neglecting an elevation angle cut-off, e.g., of 10° (Jäggi et al. 2009).
The systematic errors embedded in SLR residuals are typically mitigated by introducing station coordinate corrections, range biases, and time biases (Exertier et al. 2017;Arnold et al. 2019a;Otsubo et al. 2019). Time and range bias corrections diminish the systematic effects of the SLR observations that originate from the equipment or calibration accuracy limitations at SLR stations. Figure 1 depicts the example of systematic effects found in uncorrected SLR observations residuals to the Swarm-B satellite. SLR residuals show dependencies to different measurement conditions for different SLR stations, i.e., 10-mm offset of residuals, − 0.16 mm/° dependency to the station elevation angle, or variability of residuals with respect to station-satellite azimuth-zenith angles ( Fig. 1), i.e., the horizontal and vertical angles measured at the station to the satellite directions. Analogical dependencies also occur for Swarm-A and Swarm-C and are illustrated in the supplementary material.
Errors affecting the SLR technique can be divided into ranging system errors, timing errors, and modeling errors (Luceri et al. 2019). The origin of systematic effects in SLR residuals to LEOs is not yet explicitly explained and is under investigation by the ILRS (Otsubo et al. 2019). Figure 2 shows how possible sources of systematic effects affect the SLR observations to LEOs. In Fig. 2, top-left, range bias is a constant offset affecting range measurements, which is also seen in residuals for the Wettzell (7827) station (Fig. 1). In Fig. 2, top-right, time bias considers the offset of residuals, which changes the value and sign during one pass of a satellite above the station. An effect similar to the time bias is generated by asymmetric atmospheric conditions above the SLR stations. Tropospheric bias in Fig. 2, bottom-left, varies with respect to station-satellite elevation angle, thus may be increased at low elevation angles and is relatively small in the zenith direction because laser beam intersects different troposphere layers at different incidence angles. This bias may be seen in residuals shown for the Graz (7839) station ( Fig. 1). SLR observations are typically affected by a combination of various effects and biases, which introduce systematic patterns in SLR-based validation results of LEO orbits. The source of range bias is related to a variety of factors, including station calibration error, deficiencies in the determination of station coordinates (station height in particular), errors of the measurement system, e.g., delays in system circuits and cables, or detector performance (Pearlman et al. 1984;Appleby et al. 2016). Tropospheric bias is related to troposphere modeling, errors in the meteorological records on the station, especially the atmospheric pressure, and calibration errors (Drożdżewski and Sośnica 2021). However, the nature of tropospheric biases is similar to elevation-dependent biases, which originate from the limitations of the detectors in the case of different returning signal strength at various elevation angles or other errors, such as the interval counter or event timer bias. A time bias is caused by errors in time measurements, clock stability, event timer resolution, and station clock synchronization with respect to the UTC system (Exertier et al. 2017;Otsubo et al. 2019;Varghese et al. 2019). Barometer biases have been discovered at SLR stations, e.g., Wettzell (7827), Graz (7839), and Borowiec (7811), and can reach the level of 1 to even 5 hPa (Wang et al. 2020;Celka and Schillak 2003), resulting in elevation-dependent systematic errors (Drożdżewski and Sośnica 2021). Therefore, the ILRS investigates methods accounting for tropospheric errors in SLR and organizes a measurement campaign in 2022 using mobile barometers to compare barometer values and discover possible biases at SLR stations.
Range biases are calculated using various strategies, i.e., station-specific range biases from 1-week SLR data to spherical geodetic LAGEOS-1/2 satellites (Appleby et al. 2016), a set of individual station coordinate and range bias corrections based on 1-year data from different LEOs by fixing GPS-derived LEO orbits (Arnold et al. 2019a), or a set of particular station-satellite specific range bias corrections based on 1-year data to spherical geodetic satellites, and fixed GNSS-based orbits of LEO and Galileo satellites (Strugarek et al. 2019(Strugarek et al. , 2021bBury et al. 2021). Time biases are introduced for time synchronization between clocks used at SLR stations (Exertier et al. 2017) by using the time transfer technique (Samain et al. 2008), by calculating timing offsets from a long period of SLR residuals to LEO satellites using GPS-based orbit solutions (Arnold et al. 2019a), or by calculating time biases from each satellite pass (Otsubo et al. 2019).
Modifications in GPS receiver settings and data processing improvements enhance SLR orbit validation results of Swarm to better than 2 cm (van den IJssel et al. 2016). Montenbruck et al. (2018) and Mao et al. (2021) used integer ambiguity fixing and modeled the non-gravitational forces acting on LEO satellites using macro-models, which  Otsubo (2019) improved the consistency between the POD solution and the SLR observations to a level better than 1 cm. Progress in the GPS-based POD of LEOs requires improvement of the orbit validation processing based on SLR.
This study aims to improve the SLR-based verification of LEO POD products by eliminating systematic effects in SLR data. We discuss the modeling of systematic errors affecting SLR residuals (Figs. 1, 2) based on examples of the Swarm satellites. Using various strategies, we estimate different types of biases, i.e., range biases, tropospheric biases, and station corrections. Validation results are studied concerning station groups and measurement characteristics, such as station elevation angle. We propose considering a tropospheric bias as an elevation-dependent correction for diminishing dependencies affecting the SLR validation of LEO orbits. The performance of our methodology is tested on highquality Swarm-ABC GPS-based orbits.

Methodology
We conduct a series of Swarm-ABC orbit validation strategies to mitigate systematic errors affecting SLR residuals. We use Swarm RD POD products provided by AIUB and SLR observations from 32 ILRS stations (Fig. 3) from June 1, 2018, until August 31, 2019. As a reference, we determine a solution without estimating nor introducing additional biases. Then, we compare the validation results with experimental solutions, which consider the estimation or re-substitution of range biases, coordinate corrections, elevation-dependent corrections, and horizontal gradients. We validate all three Swarm satellite orbits, focusing on the Swarm-B satellite due to the largest number of SLR observations. The corresponding Swarm-AC results are shown in the supplementary material. Finally, we verify our method using different POD strategies, i.e., AIUB KIN and ESA RD orbits.

POD of Swarm satellites
Swarm-ABC orbits from AIUB are generated using the Bernese GNSS Software (Dach et al. 2015) in 1-day solutions with 10-s sampling. The orbit processing considers the RD approach by applying GPS antenna phase-center variations (Jäggi et al. 2009(Jäggi et al. , 2016 and is based on the ionosphere-free linear combination of phase measurements as well as pseudo-range measurements used for clock synchronization and fixed wide-lane and narrow-line phase ambiguities (Schaer et al. 2021). The non-gravitational perturbing forces are accounted for by empirical accelerations with a priori satellite macro-models employing 15 flat plates (Montenbruck et al. 2018). Detailed information about POD of Swarm satellites can be found in the description of the corresponding IANG solution in Mao et al. (2021).

SLR validation of Swarm POD products
SLR validation of Swarm orbits is based on the same background modeling for each tested solution except for different handling of biases. For all calculations, we use the Bernese GNSS Software with extended capabilities to model and estimate troposphere biases and horizontal gradients of the troposphere based on SLR. The calculated one-way ranges between stations and satellites are derived using station position from the International Terrestrial Reference Frame ITRF2014 , satellite positions based on GPS-based POD, and various corrections. The corrections consider signal propagation effects or LRA models, such as general relativistic correction; the tropospheric delay for optical wavelengths based on troposphere model and in situ meteorological measurements; the LRA correction, which considers reflections from different LRA corner cubes; satellite center-of-mass correction with respect to the LRA position obtained from the satellite orientation data, attitude, and POD products. The difference between the observed and computed ranges gives the SLR residual-an observable analyzed in this study. Table 1 provides a list of applied models used in the SLR validation of Swarm-ABC orbits. We use earth rotation parameters from the IERS-C04-14 (Bizouard et al. 2019), Swarm satellite positions from POD in IGS14 (Rebischung and Schmid 2016) frame, troposphere delay model from Mendes and Pavlis (2004), and SLR station coordinates in SLRF2014 (Luceri et al. 2019). Station and satellite reference frames are consistent with the ITRF2014 realization as they share the same origin, orientation, and scale. In each tested solution, we eliminate SLR observations with absolute residuals larger than 150 mm or the standard deviation of each station's daily residuals larger than 50 mm.
We calculated six different test solutions for SLR verification of Swarm-ABC orbits, described in Table 2. RES solution serves as a reference, in which we do not estimate any additional parameters. In the solution CRD + RB, in separate processing, we generate the 1-day Normal Equation Systems (NEQ) based on SLR data for each Swarm   Drożdżewski and Sośnica (2018). In the solution RB + TB + G we estimate four parameters per station per day: range bias, tropospheric bias, and horizontal gradients with 0.1, 1.0, and 0.1 m constraining, respectively. Tested solutions vary in terms of different modeling and the number of estimated parameters (Table 2). For example, 1-month period (July 2018), the total number of estimated parameters is 0, 116, 265, 265, 795, and 1060 for the RES, CRD + RB (re-substitution), RB-D, TB, TB + G, and RB + TB + G solutions, respectively.

Results
The validation for all tested solutions consists of residual statistics for each Swarm satellite separated into two groups of stations. We discuss the estimated range bias, tropospheric bias corrections, and horizontal gradients for the group of high-performing stations for Swarm-B. We analyze the dependency of residuals on different measurement conditions, i.e., station elevation angle, time of observations, and station-satellite azimuth-zenith angles for the group of all (32) SLR stations (Fig. 3), providing observations and the group of ten high-performing stations due to good global observing geometry and high performance (Mao et al. 2021), i.e., Yarragadee (7090)

Statistics of SLR residuals
For all tested solutions, we use on average 37,629, 106,974, and 37,063 SLR normal points (NPs) after the residual screening for Swarm-A, Swarm-B, and Swarm-C satellites, respectively. The group of high-performing stations provides 69-72% of all observations to Swarm-ABC. The number of observations to the Swarm-B satellite is almost three times higher than that of Swarm-A and Swarm-C. Close tandem configuration causes difficulties for the ILRS stations to track both Swarm-AC satellites during the same passes. Thus, in the following analysis, we focus mainly on Swarm-B residuals as a representative example. Table 3 shows that when biases are unmodeled (RES), the STDs of residuals for all stations are at the level of 14.3-15.0 mm for all Swarm satellites with a mean offset of 2.3, 0.8, and 2.0 mm for Swarm-A, Swarm-B, and Swarm-C, respectively. In the case of all solutions with modeling of systematic effects, the mean offset of residuals is reduced to near zero. When we introduce a priori range biases and station coordinate corrections (CRD + RB), or daily range biases (RB-D), the STDs of residuals are reduced to the level of 9.8-11.1 mm depending on a satellite. In the case of estimating only distance-dependent tropospheric biases (TB), the STDs of SLR residuals are reduced to the level of 7.5, 8.2, and 8.3 mm for the Swarm-A, Swarm-B, and Swarm-C, respectively. However, solutions with distancedependent corrections, horizontal gradients (TB + G), and additional daily range biases (RB + TB + G) are characterized by the lowest STD of SLR residuals at the level of 5.5-6.3, and 5.0-6.1 mm, respectively. Adding range biases to the TB + G solution insignificantly reduces the STD values, which means that range biases do not absorb any further SLR systematic errors when tropospheric biases are accounted for. Thus, the solution RB + TB + G can be considered over-parameterized because it contains parameters in the functional model that do not account for any significant systematic effects in SLR data.
In the case of ten high-performing stations (Table 3), the offsets and STDs of SLR residuals are also reduced. The CRD + RB and RB-D solutions are characterized by STDs smaller by 2-4 mm than the RES solution. The STDs for TB solution are at the level of 4.6, 6.2, and 5.9 mm for Swarm-A, Swarm-B, and Swarm-C, respectively. The solutions TB + G and RB + TB + G are characterized by the lowest STDs of 3.1 and 2.9 mm, respectively, for Swarm-A, 3.7 and 3.5 mm for Swarm-B, and 4.5 and 4.3 mm for Swarm-C (Table 3). Therefore, the STD of SLR residuals decreases below the level of 3-4 mm, which corresponds to the level of fullrate raw SLR data for high-performing stations when tracking LEO. As a result, the TB + G, RB + TB + G approaches eliminate nearly all systematic errors emerging from SLR data processing and background processing models. Figure 4 illustrates Swarm-B residual statistics in terms of median, 1st, 3rd quartile, maximum, and minimum values for all tested solutions and all stations. In the RES solution, almost all stations are characterized by positive or negative offsets of residual medians and worse than 10 mm interquartile ranges (IQR, 3rd-1st quartile). Modeling of systematic effects reduces the offset of medians to near-zero values; however, in the case of CRD + RB and RB-D for some stations (e.g., 1873, 1889, 7838) the offsets of medians are still visible at a few mm level. The IQRs for the TB solution are reduced for all stations, e.g., from 26 to 18 mm for 1874, from 13 to 8 mm for 7090, or from 32 to 10 mm for 7237 (Fig. 4). Introducing distance-dependent corrections and horizontal gradients (TB + G) and range biases (RB + TB + G) reduces SLR residuals to the largest extent, even further than in the TB and other solutions. For example, the IQRs are smaller than 10 and 2 mm for the 1874-1891 stations and the group of high-performing stations, respectively (Fig. 4).
The currently used methods for reducing systematic effects in SLR residuals to LEOs, i.e., range biases and station coordinate corrections, reduce the offset of residuals but leave their relatively large distribution, which means that station coordinate deficiencies, constant errors of SLR measurement system, and possibly orbit errors are mitigated. Considering only the elevation-dependent corrections (TB) reduces constant systematic errors, station height errors, tropospheric modeling errors, barometer biases and other unknown effects affecting SLR observations to LEOs. The reduction is valid for all stations and leads to 5-6 mm STD of SLR residuals for the high-performing stations. The additional co-estimation of horizontal gradients (TB + G), daily range biases with distance-dependent corrections (RB + TB + G) reduces SLR residuals with respect to the TB solution but increases the number of estimated parameters (Table 2).

SLR residuals with respect to different measurement conditions
We test whether different bias modeling reduces the systematic effects found in SLR residuals (Fig. 1). We focus on an example of SLR stations affected by the dependencies on different measurement conditions, such as time of observation, station elevation angle, and station-satellite azimuthzenith angles. Figure 5 shows residuals with respect to the time of observation for the Wettzell station (7827). In this case, residuals are affected by a 10-mm offset (RES, Fig. 5 topleft). The constant offset in residuals implies systematic modeling deficiencies, calibration errors, or station coordinate errors. Station 7827 joined the ILRS network in 2016 (Riepl et al. 2019) and was not included in the realization of ITRF2014. Thus, the 7827 station coordinates in SLRF2014 were not derived based on ITRF reprocessing, but using the processing of several LAGEOS passes. The offset of coordinates may be related to the deficiencies in station coordinate determination or station systematic errors. All tested modeling of systematic effects successfully reduces the offset of residuals to zero (Fig. 5). However, the distribution of residuals varies for tested solutions, i.e., STDs are at the level of 4.2 and 3.2 mm, for the TB + G and RB + TB + G, respectively, whereas for the CRD + RB, RB-D, and TB solutions, the STDs are at a level of 8.5, 7.6, and 5.8 mm, respectively. Similar offsets and dependencies for Swarm-AC are successfully reduced for tested solutions and shown in the supplementary material. Figure 6 illustrates the dependency of SLR residuals on the station elevation angle for the Graz (7839) station and tested solutions for two data sets of NPs. Wang et al. (2020) reported an offset and drift in barometer measurements installed at the Graz station, which affected troposphere correction, and in consequence, the SLR residuals analyzed in this study. The ILRS corrected the barometer data for Graz and provided the reprocessed NPs, on EUROLAS Data Center (EDC) service (Wang et al. 2020). Thus, we compared solution residuals based on uncorrected ( Fig. 6, left) and corrected (Fig. 6, right) NPs to investigate the performance of different solution strategies for the purpose of diminishing barometer-related biases. In the RES solution, the residuals are shifted toward positive values and show a dependency of residuals of over − 0.16 mm/° with respect to the elevation angle. These characteristics of residuals are caused by pressure measurement errors as well as errors in the station measurement system. The RES-N solution based on the corrected NPs is characterized by reduced residual dependency to the level of − 0.06 mm/°, 1st,3rd quartile,max.,and min. values,in mm) which means that residuals are still affected by some systematic errors. All experimental solutions remove the shift of residuals. However, for the RB-D and RB-D-N, the dependency on elevation angle remains at the level of − 0.13 and − 0.08 mm/°, respectively (Fig. 6). The TB, TB + G, and RB + TB + G successfully reduce the elevation dependency; both TB + G and RB + TB + G diminish the distribution of residuals, i.e., STDs are at the level of 3.4 and 3.0 mm, respectively. Analogical reduction occurs for corresponding solutions based on corrected NPs. Notwithstanding, for the CRD + RB, tropospheric biases are incorrectly absorbed by the station Up component and range bias corrections, reaching 11 and 14 mm, respectively. Thus, considering at least one tropospheric correction reduces the barometer-related offset and drift and better decreases the elevation-dependency of residuals than range bias and coordinate corrections. Estimating tropospheric biases at SLR stations allows for avoiding a situation where the station Up component absorbs barometer-related effects, which equals 14 mm for Graz. Figure 7 shows the dependency of residuals with respect to station-satellite azimuth-zenith angles for the Kunming (7819) station with all satellite passes over the station in the North-South or South-north directions. The dependency of residuals may include coordinate errors because this station was not considered in the ITRF2014 realization. The RES solution shows residuals of − 40 mm in the Northeast direction and 30-40 mm residuals in the Southwest direction. Most of the passes are affected by residual dependency on station-satellite azimuth-zenith angles. The CRD + RB solution successfully reduces some residual dependencies, but some passes and observations have large residual patterns. The estimation of daily biases in the RB-D solution reduces most of the dependencies at high elevation angles (with some exceptions) but increases patterns at low elevation angles (Fig. 7 top-right). The TB, TB + G, and RB + TB + G successfully reduce the residuals to the nearzero values, except for two anomalous passes. The residual values for the same passes for adjacent locations to Kunming station, i.e., Changchun and Yarragadee, and other Swarm satellites (not shown) are not characterized by increased residual values. Thus, the anomalous passes may result from gross SLR measurement errors at Kunming, e.g., due to a wrong calibration rather than from satellite orbit error.

Overview of estimated biases
We compare the medians of estimated daily range biases, troposphere biases from the RB-D, and TB solution for each Swarm satellite, and range bias and station coordinate corrections from the CRD + RB solution to Swarm-ABC, focusing on the high-performing SLR stations (Fig. 8). Values of estimated parameters in each solution are similar within the Swarm constellation. The differences do not exceed 4 mm and are related to a different number of observations for each Swarm satellite. Figure 8 shows that the medians of daily range biases (RB-D) vary within ± 20 mm between different stations. The corrections of the range bias and the station Up component (CRD + RB) are within a similar range of ± 20 mm as RB-D; however, the station horizontal coordinates corrections do not exceed ± 8 mm. The estimated  Arnold et al. (2019b) and Peter et al. (2021), whereas the resulting minor differences are associated with the different number of LEOs used in the solutions. Tropospheric corrections from TB solutions are characterized by medians of ± 6 mm and are two to three times lower than range bias medians from RB-D or CRD + RB (Fig. 8). The medians of estimated parameter formal errors are at the level of less than 2 mm for RB-D and TB solutions. Introducing troposphere biases with horizontal gradients (TB + G), and also range biases (RB + TB + G) increases the median values of estimated biases (see supplementary material) and their errors from a few mm to cm level. The mean observationsto-parameters ratio from TB, RB-D solutions of SWARM-B reaches 13, whereas the ratios for TB + G and RB + TB + G reach only 4. Thus, estimating more than one daily parameter for stations with a relatively low number of observations affects the accuracy of estimated parameters.
In general, range biases and station coordinate corrections are constant values, independent of the observation epoch and, consequently, external measurement conditions, whereas tropospheric biases and horizontal gradients depend on, e.g., station elevation angle, station-satellite distance, azimuth angle, and time. When we estimate additional parameters, i.e., troposphere biases, horizontal gradients (TB, TB + G, RB + TB + G), or coordinate corrections (CRD + RB), the co-estimated range bias values change. Range biases are correlated with troposphere zenith correction and station coordinates, especially the Up component, as shown by Drożdżewski and Sośnica (2021). However, the origin of the systematic effect may be different for each station. For example, the correction of station coordinates, especially for the Up component, changes the co-estimated range bias value (e.g., CRD + RB and 7105, 7839, or 7841 stations, Fig. 8); thus, the error is more likely related to a priori station coordinates than system measurement errors or distance-dependent biases. The expected LAGEOS-based tropospheric corrections for optical measurements are at the level of ± 5 mm, whereas horizontal gradients are within the range of sub-mm to a few mm Sośnica 2018, 2021;Drożdżewski et al. 2019).
The daily tropospheric corrections and horizontal gradients from TB + G and RB + TB + G solutions (not shown) also demonstrate the time variability but exceed dozens of mm and cannot compensate for only troposphere-related effects. Decreased accuracy of parameters originates from time-variable, uneven, and relatively low number of SLR observations. The estimated tropospheric delay possibly absorbs the distance-dependent biases and some errors in station heights because of similar geometrical characteristics in the functional model. In analogy, horizontal gradients may absorb not only the effect of troposphere asymmetry but also time biases. Different satellites at different heights and/ or 7-day solutions, as in Drożdzewski and Sośnica (2021), should be employed to improve the estimates and separate the SLR troposphere delay parameters and horizontal gradients from those that are satellite-and station-specific. Moreover, we tested a 1-day solution with an estimated daily range and tropospheric biases (RGB + TRP, not shown) in analogy to the 7-day, LAGEOS-based solution from Drożdżewski and Sośnica (2021). The estimated parameters from RGB + TRP better correspond to those from LAGEOS, but the reduction of residual statistics is relatively small and similar to the RB-D solution. The estimated parameters in solutions TB, TB + G, or RB + TB + G absorb the inconsistencies in troposphere modeling and other systematic effects affecting SLR residuals, such as intensity-dependent and distance-dependent biases. Consequently, the estimated parameters are more likely elevation-dependent biases and azimuth-dependent gradients than the pure tropospheric corrections.

Verification of method using different orbit solutions
We compare SLR validation results based on AIUB kinematic (KIN) orbits (Jäggi et al. 2016), ESA RD (ERD) orbits (van den IJssel et al. 2015), and AIUB RD orbits analyzed in previous sections for Swarm-B. We calculate the TB, TB + G, and RB + TB + G solutions to investigate if elevation-and azimuth-dependent corrections absorb the orbital errors. The effect of possible orbital errors, which is independent of SLR as associated with GNSS POD, needs to be preserved in SLR residuals to perform the correct validation of satellite orbits. We use the same methodology, period, and SLR observations as in the corresponding solutions from Table 2. Thus, the only difference is the Swarm-B POD product.
We compare the histograms of residuals ( Fig. 9 and supplementary material) and residual statistics (Table 4) from RES, TB, TB + G, and RB + TB + G solutions based on KIN, ERD, and RD orbits. In the case of three different Swarm-B POD products, the TB solution successfully reduces the range of residuals in the histograms. The percentage of the number of residuals within the ± 10 mm range is increased from 49 to 66, 56 to 70, and 68 to 91% for KIN, ERD, and RD orbits, respectively (Fig. 9), despite the same number of SLR observations used for the orbit validation.
The RES solutions with unmodeled SLR errors show the STD of residuals at 18 and 15-16 mm levels for KIN and ESA/AIUB RD orbits, respectively (Table 4). TB solution reduces the STDs of residuals to the level of 12 (14) and 11 (13) mm for KIN and ERD orbits, respectively, and to 6 (8) mm for RD orbits high-performing (all) stations. The differences between statistics for various orbits are minor in the RES solutions, whereas they became more pronounced when correcting for TB. Therefore, the differences in the GNSS LEO orbit quality can be extracted as soon as the tropospheric biases are removed from SLR data. TB + G and RB + TB + G solutions further reduce the STD of residuals to the level of 7-9 and 6-8 mm for KIN and ERD orbits, respectively. The residual offsets are reduced to near-zero in nearly all solutions with bias handling. The reduction of residuals is similar for corresponding solutions from RD orbit (Tables 3, 4), but it does not reach the level of 6-8 mm for TB or 3-5 mm for TB + G and RB + TB + G solutions.
Only the TB approach demonstrates the most PODdependent validation results. Moreover, TB solution shows the best consistency of troposphere bias estimates, with differences below 2 mm for the Swarm satellites and tested POD approaches when compared to other solutions (Fig. 8,  supplementary material). Thus, orbit-related errors are not incorporated in troposphere bias estimates. For tested POD approaches, the differences of SLR validation based on TB estimation become most apparent. Reduction of systematic effects by using troposphere corrections allows us to validate the orbit accuracy and properly assess the quality of KIN, and different RD approaches, characterized by STD of residuals of 6, 11, and 12 mm for high-performing SLR stations, and RD, ERD, and KIN orbits, respectively.

Conclusion
In this study, we search for different methods of reducing systematic effects affecting SLR residuals to LEO Swarm satellites. We investigated a series of solutions with modeling of range biases, station coordinate corrections, tropospheric biases, and horizontal gradients to reduce the residual dependency on different measurement conditions. We evaluated residual statistics of Swarm-ABC satellites for two groups of SLR stations, analyzed residuals for selected SLR stations, and compared the results with additional KIN and RD solutions.
The range biases, tropospheric biases, and horizontal gradients vary in terms of value and range. The daily estimated range bias medians or range bias and station Up coordinate corrections are within ± 20 mm in the analyzed period. The daily estimated tropospheric bias medians are within the range of ± 6 mm. Considering not only tropospheric biases per station per day but also horizontal gradients and/or  range biases decrease the accuracy of estimated parameters due to the relatively small number of SLR observations to Swarm satellites. The estimated tropospheric parameters comprise not only deficiencies in troposphere modeling but also other-station-satellite elevation-dependent effects. The observations of Swarm satellites are insufficient for the separation of troposphere-related effects from distancedependent biases because tropospheric zenith delays have similar geometry to all elevation-dependent biases, whereas horizontal gradients of tropospheric delays absorb not only the effect of atmospheric asymmetry but also errors in a priori horizontal station coordinates. However, elevationdependent biases are the most significant biases embedded in SLR observations. Nevertheless, the estimation of daily tropospheric biases (TB solution) reduces the STDs of SLR residuals of AIUB RD orbits from 10 (15) mm to the level of 6 (8) mm for the high-performing (all) stations. This approach reduced the SLR residuals to Swarm orbits by 3 mm (29%) with respect to the solutions considering daily range biases (RB-D), range biases and station coordinate corrections (CRD + RB), or results by Mao et al. (2021) for a Swarm-C satellite. Range biases reduce only a constant offset of residuals and deficiencies in station coordinate determination. Increasing the number of elevation-and azimuth-dependent corrections further reduces the SLR residuals and dependencies affecting SLR stations. However, the number of estimated parameters should be kept at a minimum level; otherwise, the functional model becomes overparameterized, and the parameters may also absorb orbit-related effects that should be extracted from the SLR validation. On the other hand, the NPs to LEOs that are freed from SLR-based and orbitrelated effects are also required for SLR-based determination of high-quality station coordinates and global geodetic parameters (Strugarek et al. 2019(Strugarek et al. , 2021bLi et al. 2021). The compromise between the number of estimated parameters and the most efficient reduction of SLR errors is considering only elevation-dependent bias (TB solution) with one parameter per station per day.
Estimating tropospheric biases removes completely the recently discovered barometer error for the Graz station. The results for Graz are the same for barometer corrected and uncorrected SLR NPs when estimating tropospheric biases. Thus, adding one tropospheric bias per station successfully removes barometer biases embedded in SLR data. Without correcting the tropospheric bias, the barometerrelated error of 14 mm is wrongly absorbed by the Up station component.
Introducing troposphere biases (TB solution) in SLR solutions allows for proper orbit validation, as well as different GPS POD approaches. The STD of SLR residuals is 6, 11, and 12 mm for AIUB RD, ERD, and AIUB KIN orbits, respectively, for high-performing SLR stations; thus, the differences of orbit quality become unequivocal when correcting for the troposphere biases. The estimated corrections from the TB approach are consistent within the Swarm constellation and tested POD approaches. Therefore, the orbital errors are not absorbed in troposphere bias estimates. Troposphere corrections better reduce most SLR errors than range biases due to the capability of absorbing the constant part of a bias and elevation-dependent biases. Moreover, SLR observations corrected by tropospheric biases allow for an extraction of differences between the orbit qualities derived using different POD methods because the dominating systematic errors in SLR measurements are removed, and thus, the SLR residual differences become more pronounced between different orbits.
Further studies should include an increased set of satellites, longer observation periods and solutions, e.g., 7-day solutions. The increase of SLR observations by using more LEOs or other types of satellites, including GNSS, is crucial for differentiating the particular source of error from the SLR observation model and determining SLR-based station coordinates, geocenter coordinates, or earth rotation parameters. Methods of modeling systematic effects from this study are expected to improve SLR validation to different LEOs and allow for deriving high-quality geodetic parameters and SLR station coordinates that are freed, to the greatest extent, from systematic errors. Data availability The SLR observations used for this study can be downloaded from the Crustal Dynamics Data Information System ftp:// cddis. nasa. gov/ slr/ produ cts and the EDC ftp:// edc. dgfi. tum. de/ pub/ slr/ produ cts. The Swarm POD products are available from the ESA ftp:// swarm-diss. eo. esa. int/ Level 2daily/ Entire_ missi on_ data/ POD/, AIUB ftp.aiub.unibe.ch/LEO_ORBITS/SWARM/RL03/, and authors on a reasonable request.

Supplementary Information
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