Tropospheric and range biases in Satellite Laser Ranging

The Satellite Laser Ranging (SLR) technique provides very accurate distance measurements to artificial Earth satellites. SLR is employed for the realization of the origin and the scale of the terrestrial reference frame. Despite the high precision, SLR observations can be affected by various systematic errors. So far, range biases were used to account for systematic measurement errors and mismodeling effects in SLR. Range biases are constant for all elevation angles and independent of the measured distance to a satellite. Recently, intensity-dependent biases for single-photon SLR detectors and offsets of barometer readings and meteorological devices were reported for some SLR stations. In this paper, we study the possibility of the direct estimation of tropospheric biases from SLR observations to LAGEOS satellites. We discuss the correlations between the station heights, range biases, tropospheric biases, and their impact on the repeatability of station coordinates, geocenter motion, and the global scale of the reference frame. We found that the solution with the estimation of tropospheric biases provides more stable station coordinates than the solution with the estimation of range biases. From the common estimation of range and tropospheric biases, we found that most of the systematic effects at SLR stations are better absorbed by elevation-dependent tropospheric biases than range biases which overestimate the total bias effect. The estimation of tropospheric biases changes the SLR-derived global scale by 0.3 mm and the geocenter coordinates by 1 mm for the Z component, causing thus an offset in the realization of the reference frame origin. Estimation of range biases introduces an offset in some SLR-derived low-degree spherical harmonics of the Earth’s gravity field. Therefore, considering elevation-dependent tropospheric and intensity biases is essential for deriving high-accuracy geodetic parameters.


Introduction
Satellite Laser Ranging (SLR) is a space geodetic technique indispensable for deriving geocenter coordinates, Earth oblateness, standard gravitational constant parameter (GM), the scale of the terrestrial reference frame, validating general relativity effects, determining and validating satellite orbits, and many other applications (Pearlman et al. 2019a). Despite high precision, SLR solutions are still affected by a bunch of systematic errors emerging from limitations in background models or imperfect calibration and meteorological measurements, all of which may bias the final solution and geodetic parameters. All systematic errors require proper identification and handling or estimating additional parameters to compensate for their impact. The International Laser Ranging Service (ILRS, Pearlman et al. 2019b) coordinates collecting, processing, and distributing SLR observations and products. ILRS is currently struggling with fulfilling the Global Geodetic Observations System's (GGOS) requirements which expect the 1 mm of station coordinate stability and the 0.1 mm/year accuracy for station velocities (Plag and Pearlman 2009). Thus, it is overwhelmingly important to investigate the sources of all possible systematic errors.
The scientific discussion about handling systematic errors and biases in SLR is ongoing for over 35 years. Pearlman et al. (1984) proposed to divide the source of systematic errors into the following groups: ranging appliance errors, group cluster calibration as well as synchronization issues, hardware limitations, timing errors (station event timers), and finally, modeling errors (Exertier et al. 2019). Altamimi et al. (2016) indicate that the discrepancy between scale factors derived from Very Long Baseline Interferometry (VLBI) and SLR in the International Terrestrial Reference Frame (ITRF) 2014 at the level of 7-8 mm could be caused by range biases in SLR (Appleby et al. 2016) as well as by antenna gravity deformation in VLBI (Sarti et al. 2009).
Currently, the ILRS community, especially the Analysis Standing Committee (ASC), establishes an approach based on the estimation of one scalar value-a range bias-to absorb the total ranging error with the same value independently from the ranging direction and elevation angle (Appleby et al. 2016;Rodríguez et al. 2019). However, identifying the real source of errors contributing to the determined range biases at SLR stations should also be conducted. Luceri et al. (2019) describe and classify the source of errors affecting on-station measurements with an extended discussion of the impact of range biases on SLR results, including the station coordinates and the scale of the reference frame. However, the errors caused by the effects related to the troposphere delays are not discussed by Luceri et al. (2019).
SLR is the only one space geodetic technique in which the zenith total delay (ZTD) and the slant tropospheric delays are derived based solely on in situ meteorological observations used as fixed values. In SLR, ZTD is based on atmospheric pressure and humidity and the geographical location of the SLR station, whereas the mapping function and thus also the slant delay employ the temperature records collected simultaneously with the laser ranging data. This approach is explained by the relatively low number of observations in SLR in the early years of operating, which prevented the possibility of direct ZTD determination, and a minor sensitivity of optical wavelengths to the wet delays when compared to microwave techniques.
Meteorological sensors installed at SLR sites may provide biased records or may deteriorate over time, resulting in wrong troposphere corrections applied to SLR observations. Recently, the SLR stations equipped with more than one barometer, such as Wettzell and Graz, revealed that both sensors might show a bias of even 0.78 hPa (Wang et al. 2020), which translates to ZTD bias of 2 mm in the zenith and over 8 mm at the elevation angle of 10°. Unexpected occurrence of meteorological sensors issues at SLR stations includes consequences, such as that barometric sensors can drift over time during the aging of devices, whereas a change of the barometric sensor can induce a barometric offset. Additionally, malfunctioning of barometric record updates may result in constant barometric pressure values. When the barometric sensor is not located at the same height as the system reference point, and the height difference is unmodeled in the data processing, the barometric readings may be biased by about 0.1 hPa for every 1 m of the height difference. All of these issues need careful handling to avoid introducing biases in employed SLR tropospheric corrections, which cannot be guaranteed for all SLR stations.
The approach based on the estimation of ZTD corrections, i.e., the corrections to the a priori total delays, has never been considered in SLR. Such an approach can potentially remove the barometric issues in SLR observations by introducing an additional parameter to the functional model. So far, the tropospheric biases were absorbed by the estimation of range biases; however, the nature of both delays is different. Range biases are constant delays acting irrespectively of the elevation angle, whereas tropospheric biases are smallest in the zenith direction and substantially increase for low elevation angles. The estimation of troposphere corrections may require fewer parameters than the solutions with the estimation of range biases because range biases are typically estimated for each satellite separately. In contrast, the ZTD corrections absorbing the biases of barometer readings can be estimated as one parameter for all satellites in the same time interval as typically estimated for range biases.
In this paper, we propose a new approach based on the estimation of troposphere delay correction for SLR solutions, which is a standard approach in microwave techniques. In microwave geodetic techniques, such as the Global Navigation Satellite System (GNSS) or VLBI, the most efficient way to determine the wet part of the delay is the estimation of the tropospheric wet delay due to difficulties in providing high-accuracy a priori delays based on numerical weather models (Balidakis et al. 2018;Teke et al. 2013;Wilgan and Geiger 2019). In GNSS or VLBI, the hydrostatic part of the delay is typically based on numerical weather models or calculated for a station using in situ measurements (Boehm et al. 2007), whereas the wet part is estimated as corrections between the retrieved total delay and the a priori hydrostatic part (Teke et al. 2013). To stabilize the GNSS or VLBI solutions, the troposphere delay is not estimated for each slant observation, but instead, the wet delay is estimated for the zenith direction and projected onto the requested elevation angle using corresponding mapping functions for the hydrostatic a priori part and the estimated wet part. The estimated wet zenith delay from microwave techniques can be used for different purposes, such as assimilation to the numerical weather model predictions (Hanna et al. 2019;Trzcina and Rohm 2019;Le Marshall et al. 2019). However, the estimated tropospheric parameters are typically strongly correlated with receiver clocks and the vertical component of station coordinates in GNSS and VLBI solutions. Boisits et al. (2020) compared the meteorological atmospheric pressure derived from numerical weather models and observations collected at SLR stations and discovered large discrepancies for some stations, e.g., 11.2 hPa for Altay (1879), 15.6 hPa for Zelenchukskaya (1889), and 83.7 hPa for Badary (1890). Wijaya and Brunner (2011) considered using two SLR frequencies for correcting the slant delays; however, the accuracy of the SLR measurements on infrared and green frequencies would have to be substantially improved with respect to the current capabilities of SLR stations. Drożdżewski and Sośnica (2018) calculated the horizontal gradients of the tropospheric delay based on SLR observations to LAGEOS satellites to account for the asymmetry of the atmosphere. Drożdżewski et al. (2019) employed the horizontal gradients from the external Potsdam Mapping Function model to improve the SLR solutions and found that the gradients can improve the consistency of SLR-derived Earth rotation parameters when compared to other space geodetic techniques. However, the zenith troposphere delay corrections have never been directly estimated from SLR observations until now.

Processing strategy
We process nine years of SLR observations to LAGEOS-1 and LAGEOS-2 satellites (Pearlman et al. 2019a) for the period 2010.0-2019.0 using 7-day orbital arcs, in a similar way to the standard ILRS ASC solutions. We follow the a priori models recommended by the Conventions 2010 of the International Earth Rotation and Reference Systems Service (IERS, Petit and Luzum 2010). In all solutions, the meteorological in situ observations collected at SLR stations are used for calculating the total zenith delay based on the model developed by Mendes and Pavlis (2004) and the mapping function developed by Mendes et al. (2002) as recommended by the IERS Conventions 2010.
We determine the 7-day solutions with the estimation of SLR station coordinates, LAGEOS orbits, low-degree gravity field coefficients up to degree and order 2, geocenter coordinates, and Earth rotation parameters (ERPs): pole coordinates and the Length of Day (LOD). Only ERPs are estimated as piece-wise linear parameters with 1-day intervals, whereas all other parameters are parameterized as constant over 7-day solutions. For satellite orbits, we determine 6 Keplerian parameters with constant and onceper-revolution empirical components only in the along-track direction. Other empirical orbit parameters can be correlated with gravity field coefficients, such as once-per-revolution cross-track parameters and C 20 (Sośnica et al. 2014. In the operational ILRS solutions, the gravity field coefficients are not estimated; once-per-revolution orbital parameters in cross-track are estimated, whereas the ERPs are parameterized as piece-wise constant. However, in this study, we check what is the impact of the troposphere modeling in SLR solutions on the low-degree gravity field coefficients, in particular, the oblateness term; thus, we decided to consider these parameters. For the a priori reference frame, we use the ILRS special realization of ITRF2014, i.e., SLRF2014, which shares the origin, scale, and orientation with ITRF2014. The network constraining follows the list of core ILRS stations for the no-net-rotation and nonet-translation constraints with the 25-mm criterion for the coordinate stability for each week following the approach described by Zajdel et al. (2019). Figure 1 illustrates the impact of the geometry of the observation distribution on the estimated values of station heights, range biases, and ZTD in SLR solutions. Having only SLR observation in the zenith direction, it would be impossible to separate these parameters. To decorrelate ZTD from range biases and station heights, observations collected at different zenith angles ( z ) are needed because station height and troposphere parameters depend on the zenith or the elevation angle (e = 90 • − z) , whereas range biases assume the same values irrespectively of the elevation angles. In SLR solutions, the stable station heights are the most desired parameters. From stable station coordinates, the scale of the reference frame can be retrieved, as well as the geocenter coordinates and other geophysical parameters. Fig. 1 Dependency (correlations) between station heights, troposphere zenith delay, and range biases in the SLR solutions. Adopted version based on an analogy to GNSS solutions by Rothacher (2003) 100 Page 4 of 18 Estimated range biases shall absorb all unmodeled effects. However, the nature of systematic errors can be different and not necessarily constant for all elevation angles. Table 1 shows the correlation coefficients between station heights, ZTD, and range biases for different elevation cut-off angles. Most of the SLR stations track the LAGEOS satellites down to 20° or 15°, which results in strong correlations. Despite large correlations between station heights and ZTD, they seem to be better separable than the range biases and station heights. Therefore, we will test all possible solutions: with the estimation of range biases, with the estimation of ZTD, and with the simultaneous estimation of ZTD and range biases to assess their impact on resulting station heights.
To estimate ZTD corrections, proper partial derivatives of slant troposphere delays must be introduced to the functional model described by the first-design matrix. We employ the standard mapping function to estimate ZTD corrections. An assumption of expanding the mapping function with the 1∕ cos (z) term results in the model proposed by Mendes et al (2002) with a 1 , a 2 , a 3 parameters for optical wavelengths, such as that from Eq. 1: In this paper, we propose an approach similar to that commonly used in microwave techniques. The final values of ZTD are calculated in the parameter adjustment process. We add one additional parameter to the first-design matrix with the so-called troposphere bias. This approach introduces an elevation-dependent scaling factor MF(z) for each observation (Hobiger and Jakowski 2017).
We use the partial derivatives of ZTD in the function of slant total delay to establish the troposphere correction in the least-squares adjustment process. The primary difference between the standard GNSS troposphere corrections and the approach in this paper lies in the fact that in GNSS the hydrostatic part of the delay is used as the a priori delay, and the wet delay is calculated. In this paper, we estimate corrections to the total a priori delays; thus, we use the wet and hydrostatic delays as a priori values. Tregoning and Herring (2006) and Hobiger and Jakowski (2017) mentioned that any adjustment required to the a priori value of zenith hydrostatic delay is also estimated as part of the correction to zenith wet delay, so the correction introduced in the SLR technique is the total amount of errors caused at the site related to all troposphere error sources.
In the design matrix, the tropospheric correction is a function of the elevation angle, whereas the range bias contribution is constant (Fig. 2). Observations at lower elevation angles are thus essential to decorrelate the station heights, range biases, and troposphere parameters. ZTD increases substantially with the decreasing elevation angle (Fig. 2). In the standard ILRS solutions, only station height or station heights and range biases are estimated. However, when a troposphere bias exists, the range bias cannot compensate for its value because ZTD has a different elevation dependency in the functional model. Station heights may thus be affected by a missing correction.
The separation of biases to a specific group could improve the investigation of sources of systematic errors in SLR, such as time biases, range biases, detector signature effects, and intensity biases . This paper addresses the issues of systematic errors that can be detected and how their effects on geodetic products could be mitigated empirically. The troposphere biases have a similar nature to intensity biases; thus, troposphere delay parameters may absorb both types of systematic errors. Appleby et al. (2016) showed that the station-related range biases are significant factors limiting the accuracy of SLR-derived coordinates and other crucial parameters, such as the global scale. Moreover, range biases are indistinguishable from satellite center-of-mass corrections (CoM) and can absorb wrong CoM values for spherical satellites (Otsubo and Appleby 2003;Rodríguez et al. 2019). Neglecting the Table 1 Correlation coefficients between the station heights (h), range biases (RGB) and zenith total delays (ZTD) assuming even distribution of SLR observations above the SLR stations, i.e., perfect observational geometry and different values of the elevation cutoff angles. Calculated using the method described by Beutler et al. (1987) Elevation cut-off angle

Fig. 2
Dependency between station heights, range biases, and troposphere biases as a function of the elevation angle of observations systematic effects related to range biases can significantly deteriorate the estimated station heights. We examine different handling of the slant delay (SD) considering the impact of the height difference (h), a priori ZTD, estimated tropospheric correction (TRP), and range bias parameter (RGB). Therefore, four different SD handlings in the functional model are considered (Eq. 2-5), which can be written as: where MF(z) describes a mapping function based on the standard FCULa (Mendes et al. 2002) approach (Eq. 1) employing in situ meteorological observations, ZTD denotes the zenith total delay computed following the IERS 2010 Conventions (Petit and Luzum 2010), also based on in situ parameters (Mendes and Pavlis 2004).
The experiment was conducted in the following way (see Table 2); firstly, we derive the standard reference solution (STD); secondly, we derive a solution with an estimation of troposphere delay (TRP) with 7-day time interval constrained to 0.01 m; thirdly, we derive a solution with the estimation of biases with a priori constraining of 0.01 m; and finally, we derive a solution with the estimation of troposphere delay correction as well as range biases with a priori sigmas at the level of 0.01 m (TRP + RGB). The estimation of ZTD corrections in SLR is realized by a dedicated mapping function (Eq. 1) and strongly depends on the geometry of the observations (Table 1, Fig. 1).
In our approach, the temporal variations of the tropospheric delays are accounted for by in situ meteorological measurements collected together with SLR normal points, whereas the biases of meteorological instruments are captured by constant offsets that are estimated. Estimating TRP every 1 h or 2 h, as typically done in GNSS solutions, would not be possible in SLR solutions because of the low number of SLR data or sometimes the lack of any observations for a few days in 7-day solutions. In our solutions, the highfrequency troposphere delay variations are captured by real meteorological observations without deteriorating the vertical component of station coordinates because of the sparse estimation of TRP parameter-once per 7 days. However, the estimated TRP parameter accounts for the possible offsets of the meteorological readings or improper installation of the device. Therefore, the resulting solution should both be stable because of the minimum number of estimated parameters and should be freed from meteorological-related systematic errors.

Impact of artificial pressure bias
First, we analyze the capability of SLR solutions to reconstruct the tropospheric biases. We introduce an artificial pressure bias of 5 hPa for every observation at all SLR stations. The 5 hPa bias is equal to 11.4 mm of troposphere delay differences in the zenith when using Mendes and Pavlis (2004) model and standard meteorological conditions. The TRP offset of 11.4 mm in the zenith yields 40.6 mm of error for observations provided at 10 degrees of elevation angle. This simulation emphasizes how the pressure biases are absorbed by other estimated parameters, which are correlated with troposphere delays, especially the station heights. Moreover, this experiment can answer the question of whether the estimation of range biases is suitable to handle errors related to the tropospheric delay in SLR measurements.
We derive six solutions and compare them to the standard STD solution with no artificial bias introduced. In the first three solutions, we estimate only one of three analyzed parameters (see Table 3): only RGB, only TRP, and only station coordinates (CRD) to guarantee that only one parameter absorbs the tropospheric bias. Each of the solutions is denoted by the name of estimated parameters. Two solutions consider estimation of two parameter types RGB + CRD or TRP + CRD, and one solution with the estimation of all parameter types: TRP + RGB + CRD. All parameters are correlated (see Fig. 1); thus, we scrutinize what part of the artificial tropospheric bias is absorbed by the troposphere parameter and what part leaks into station coordinates and range biases. In this part, we use a priori station coordinates derived from the previously estimated solution STD.
The results of the experiment are shown on the example of two core stations Yarragadee (Australia) and Graz (Austria). In solution RGB, the total value of the estimated range bias equals − 19.2 and − 21.2 mm for station Yarragadee and Graz, respectively (see Table 3). When we estimate only the troposphere correction, the reconstructed TRP values are between − 11.1 and − 11.8 mm for both stations. When estimating only station coordinates, we obtain offsets at the level for Yarragadee and Graz, respectively. When the meteorological observations and SLR observations are affected by a tropospheric bias, the estimation of range bias does not solve the problem. The up component of SLR station coordinates changes significantly, but the actual position of station heights is somewhere between the solutions with and without estimating the range bias. Adding a range bias as a parameter does not resolve the issue of troposphere delay but instead leads to an overestimation of the total effect with a correction with an opposite sign but similar magnitude in the station heights and thus also in the scale of the reference frame and all station-related parameters, such as geocenter coordinates. Figure 3 shows the time series of estimated RGB, TRP, and station heights from all solutions with a bias of 5 hPa for Yarragadee. The estimated values of range biases are in close range to the tropospheric delays projected onto a low elevation angle rather than the value of the simulated artificial error in zenith. Even if we estimate range biases, station coordinates are still burdened by atmospheric pressure errors. In the solution TRP + CRD, we observe that troposphere delay absorbs more than 90% of the artificial tropospheric bias. The reconstructed TRP value is at the level of a priori artificial error, whereas the effect on station coordinates is absorbed in 75-90% with the residual effect of 2.5 and 3.3 mm for Yarragadee and Graz, respectively. Finally, in the solution TRP + RGB + CRD, we observe a negative offset of station heights for both stations at the level of − 4.9 and − 3.7 mm. The range biases estimated commonly with troposphere corrections are unstable and very noisy because of the correlations (see Fig. 3). Due to the parameter leakage effect, the mean value of the estimated range biases exceeds the level of 10 mm for Yarragadee instead of being fully absorbed by the estimated tropospheric bias. In the solution with the common estimation of tropospheric biases, range biases, and station coordinates, the station heights are affected up to 4.9 mm, which is 40% of the introduced bias or 15% of the station height offset from the CRD solution.
From this experiment, we conclude that the estimation of range biases cannot properly account for the tropospheric biases in SLR observations. Adding the range bias as an additional parameter leads to an opposite offset of estimated station heights with a similar order of magnitude. Only the solutions with the estimation of

Results
Now, we provide results from processing the SLR observation to LAGEOS satellites based on real data, as described in Sect. 2.1. We use the a posteriori sigma of unit weight, that is, the mean global error of the solution, as an indicator to validate the quality of solutions. The sigma decreases when an additional parameter absorbs some effects embedded in observation residuals, and thus the functional model better describes the parameter-observation dependencies. Each sigma value is normalized by the degree of freedom; therefore, the sigma value increases when adding a parameter that does not change anything in the solution. One has to bear in mind that two parameters per station are added in the RGB solution, whereas, in the TRP solution, one parameter per 7-day solution is estimated. Figure 4 shows that the estimation of troposphere delay correction results in improvements of the a posteriori sigma of unit weight by 7% from 11.4 to 10.6 mm. The improvement in the RGB solution with two additional parameters is smaller, although twice more additional parameters are estimated in RGB than in TRP. The solution TRP + RGB also decreases the a posteriori sigma of unit weight, but the improvement is much smaller and equals about 2% when compared to reductions caused by adding TRP parameters (Fig. 4).

Analysis of range biases and TRP correction
for SLR stations Figure 5 shows the total number of observations (top) and the mean value of long-term biases from RGB and RGB + TRP solutions, separately for LAGEOS-1 and LAGEOS-2 (middle). Figure 5 reveals what part of the troposphere biases can be absorbed by the estimates of range biases. A separation between RGB and TRP biases is possible for stations with dense observation sky coverage and regular distribution of measurements. For stations located nearby mountains, where dynamical changes of atmospheric parameters at different layers of the atmosphere are  observed, the troposphere correction constitutes more than half of the average long-term value of range biases, such as in the case of stations Zimmerwald or Graz, which are located at the elevation of 951 m and 495 m, respectively. For the station Graz, the mean RGBs without estimating TRP are at the level of 4.5 mm, whereas when estimating TRP, the RGB decreases to about 0.5 mm for both LAGEOS satellites. The mean TRP correction reaches 2.6 mm for solutions TRP and TRP + RGB, (see Fig. 5 bottom). For station Graz, a recently discovered pressure bias occurred since 2015 (Pavlis and Kuźmicz-Cieślak 2020); thus, this station will be analyzed in detail in subsequent parts of the paper together with the station Wettzell, where also a bias was discovered which could not be attributed to any specific parts of the SLR system (Riepl et al. 2019). Figure 5 shows that the mean TRP value only marginally changes when adding the estimation of RGB, which means that RGB corrections absorb the TRP offsets only to a minor extent. Noteworthy, in the solution with the simultaneous estimation of range biases and troposphere (RGB + TRP), the dominating part of the bias is absorbed by the TRP. The value of the estimated RGB is about four times smaller in the solution where the TRP is estimated compared to the solution without RGB. This suggests that most of the systematic effects at the Graz station should be assigned to the troposphere biases and not to the range bias. For stations Concepcion, Haleakala, Graz, Mt Stromlo, and Zimmerwald, the estimated troposphere offsets significantly change the mean range biases. Figure 5, bottom, illustrates the mean values of troposphere corrections derived from solutions TRP and RGB + TRP. Estimating troposphere biases with range biases leads to an insignificant decrease in troposphere biases (Fig. 5). Only for stations McDonald and Monument Peak with a low number of observations, the mean differences are at the level of 0.5 mm.
For station Graz, the estimation of troposphere delay reduces the mean range biases from 4.5 to 0.5 mm and from 4.5 mm to 0.8 mm for LAGEOS 1 and LAGEOS 2, respectively. The separation of range biases from troposphere biases may help further investigate systematic errors at the stations or be useful for validating new CoM corrections (Otsubo and Appleby 2003;Rodríguez et al. 2019). Figure 6 shows the time series of estimated RGB values for LAGEOS-1 and LAGEOS-2, TRP corrections, as well as the height component for the Graz station. For the troposphere delays from the TRP and TRP + RGB solutions, no significant differences between these two series can be found. The difference of estimated tropospheric corrections never exceeds 0.5 mm between solutions TRP and TRP + RGB. Thus, the estimated troposphere corrections are stable, whereas the RGB values substantially vary between the solutions with and without estimating tropospheric corrections. The vast majority of TRP estimates for the Graz station are positive, suggesting that a barometer bias may exist for this station. In the RGB + TRP solution, the estimated range biases oscillate around zero, whereas in the RGB solution, the estimated range biases are shifted toward positive values. The positive range biases imply that Graz may be affected by a range bias; however, the RGB + TRP solution assigns the bias to the troposphere-related bias. We should keep in mind that the troposphere biases have a similar nature to the intensity-dependent biases because both of them increase with the measured distance and can only be separated when using satellites orbiting at different heights above the Earth's surface. Figure 6, bottom, depicts the time series of station height for each of the solutions. The standard solution (in red) shows a large-scale negative offset and, after 2014, a drift of the vertical coordinate component of the Graz station, which is related to the recently discovered barometer issue. Estimating RGB or TRP or RGB and TRP corrections account for possible biases in on-site temperature or pressure records that substantially affect the station heights. Such biases of meteorological parameters are typically discovered after a few years or even never disclosed at all. Estimating RGB values systematically changes the station vertical component toward positive values. However, the solution with the estimation of TRP and RGB corrections is closer to the solution with the estimation of solely the TRP correction, which suggests that the estimation of only RGB values overestimates the total effect.
In this study, we employ the SLR observations as they were originally published in the ILRS database. It is worth noting that the barometer bias in Graz was discovered in 2020 (Wang et al. 2020) so that five years of Graz observations were replaced on the ILRS servers. In this study, we employ the original (biased) SLR observations to highlight how the proposed method works in terms of identifying and handling undiscovered systematic biases in meteorological data.
For the Wettzell observatory, a discrepancy in meteorological data has been found between the Satellite Observing System Wettzell (SOS-W, 7827) and the Wettzell Laser Ranging System (WLRS, 8834) (Riepl et al. 2019;Pavlis and Kuźmicz-Cieślak 2020). Figure 7 shows a prominent negative secular rate of the TRP corrections, which substantially affect the vertical component of station coordinates in the STD solution in the form of large-scale positive drift. The variations of the station vertical component are strongly mitigated when correcting for TRP or RGB biases.
The spectral analysis of the time series of RGB values estimated with and without TRP corrections (Fig. 8) shows that the noise of the RGB values is much smaller in the RGB + TRP solution. Moreover, the co-estimation of TRP reduces the semi-annual and annual signals from the time series of range biases. The existence of annual signal in the RGB series suggests that RGB values absorb not only the instrumental errors at SLR stations but also some effects of geophysical origin. TRP corrections have an environmental origin, thus may vary with seasons as opposed to RGBs. Therefore, we can conclude that the solution TRP + RGB successfully separates the range biases of entirely instrumental origin from the tropospheric biases of the quasienvironmental origin.

Validation of troposphere delay corrections in SLR measurements
Since 2015, the station Graz has had problems with barometer measurements which reached 0.78 hPa (Pavlis and Kuźmicz-Cieślak 2020). This malfunctioning barometer is best opportunity to validate whether the troposphere correction is sensitive to retrieve those issues. For that purpose, we separate the time series of range biases and troposphere delay correction between 2015.0 and 2019.0 (Fig. 9). We observe that the estimation of troposphere delay correction reduces the mean offset of range bias from 4.1 to 0.2 mm and removes the secular rate from 0.5 mm/year to 0.0 mm/year. The mean offset of the troposphere delay correction is equal to 1.7 and 1.6 mm for solutions TRP and TRP + RGB, respectively (Table 2), which corresponds to approximately 0.7 hPa of pressure bias. This result confirms that estimating troposphere delay corrections properly resolves discrepancies related to ZTD. Riepl et al. (2019) discovered an unexpected differences between range biases from station 7827 SOS-W and 8834 WLRS. Figure 9, bottom, illustrates the time series of range biases from solutions RGB and TRP + RGB with the mean offset of range biases derived from solution RGB equal to − 7.6 mm. We observe also an annual fluctuation with the amplitude equal to 1.6 mm and a secular rate equal to − 0.6 mm/year (see Table 4). The estimation of troposphere delay corrections reduces the mean offset to − 1.8 mm, which agrees well with the value provided by Riepl et al. (2019), who reported a bias of unknown origin of − 1.7 mm. The estimation of troposphere corrections also reduces

Impact of troposphere corrections and range biases corrections on station heights
This section discusses the impact of the estimated RGB and TRP parameters on the station heights. Figure 10 illustrates the distribution of estimated station coordinate north, east, and up components with respect to the a priori SLRF2014 reference frame from the STD, RGB, TRP, and RGB + TRP solutions. No significant differences can be found in the North and East components of the station coordinates, as changes of these are typically related to the horizontal gradients of the troposphere delays and not to the zenith delays ). However, significant differences are found for the height component, which is highly correlated with troposphere correction and range biases. Figure 10 depicts the long-term consistency with respect to the a priori coordinates from SLRF2014. The TRP  solution is closer to the a priori values of station coordinates for stations Tahiti, Concepcion, Grasse, and Zimmerwald than other tested solutions. However, the a priori SLRF2014 does not consider TRP biases, but the scale of the SLRF2014 and ITRF2014 is the mean scale as derived from VLBI and SLR. As shown in Fig. 10, the STD solution underestimates the station heights, and consequently, also the global scale for most of the core stations, such as Zimmerwald, Graz, Mount Stromlo, Yarragadee, Hartebeesthoek, Haleakala. The SLR scale discrepancy has been recently improved by the new center-of-mass correction model , which is employed in this study, but was not used during the preparation of ITRF2014. When RGB or TRP biases are estimated, the up component is shifted toward positive values, and thus the global scale becomes more consistent with the VLBI or ITRF2014 scale. The largest positive shift is for the RGB solution, whereas the simultaneous estimation of RGB and TRP provides the vertical components with median values similar to those from the TRP solution or between the TRP solution and the RGB solution. This suggests that the estimation of only RGB provides to an overestimation of biases' effect at SLR stations. The TRP + RGB solution contains two types of bias parameters in the functional model, thus allows for a proper distinction between the Fig. 10 SLR station coordinate repeatability decomposed into the North (N), East (E), and Up components with respect to SLRF2014. The box describes the 25th and 75th percentiles, the horizontal central red line describes the median value, the whiskers (black dashed line) describes the most extreme data points without outliers bias that is constant, independent from the elevation angle, and the bias that depends on the measured distance.
The up component of station coordinates is correlated with the troposphere corrections and range biases. Now, we need to address the issue of whether any deterioration of SLR station coordinates occurs when adding parameters that are correlated with station heights in terms of the noise and stability of derived parameters. Figure 11 illustrates the repeatability of the SLR station coordinates in the form of interquartile ranges of the estimated station vertical components.
The solution RGB significantly deteriorates the station repeatability with regard to solution STD for more than 80% of analyzed stations. Increased noise of the station vertical component because of co-estimated range bias is expected and discussed by, e.g., Rodríguez et al. (2019) and Couhert et al. (2020). In the case of the TRP solution, we observe an improvement in 50% of analyzed stations and no significant difference, that is below 1 mm, in 19% of cases w.r.t. to the STD solution. From Table 1, we learn that the correlation between station height and troposphere parameters is weaker than the correlations between station heights and range biases. Adding the TRP parameter removes systematic effects related to tropospheric and range-dependent bias and does not substantially deteriorate the estimates of station heights. For some stations, such as Yarragadee or Greenbelt, the TRP solution provides more stable station height estimates than the standard solution. Thus, TRP parameters can be estimated without deteriorating the station heights as opposed to the RGBs, which have a negative impact on the repeatability of station coordinates. Beutler et al. (1987) found that the estimation of absolute troposphere biases significantly impacts the scale factor of the reference frame realized using GNSS observations. Figure 12 shows that one of the crucial ITRF parameters derived from SLR measurements, i.e., the global scale , is also burdened when neglecting ZTD corrections. Thanks to the new center-of-mass corrections, the median offset in solution STD is equal to 0.9 mm with respect to the a priori scale from SLRF2014/ITRF2014. The estimation of troposphere corrections reduces the scale offset to 0.6 mm, whereas in the solution RGB, the mean offset is at the level of − 0.9 mm. Finally, in solution TRP + RGB, the offset equals to 0.0 mm. Before 2014, the discontinuities in the time series of station coordinates were verified and four techniques of space geodesy were connected using local ties in the framework of the ITRF2014 realization . After 2014, we observe the aging of the reference frame realization because new discontinuities and new values of post-seismic deformation in ITRF2014 are not provided. Therefore, a large discrepancy in the scale can be found between the SLR scale and that from ITRF2014, especially in the STD solution.

Scale
We obtain the highest stability of the scale realization from the TRP solution. The standard deviation of scale estimates is equal to 3.4, 3.0, 3.9, and 3.4 mm for solutions STD, TRP, RGB, and TRP + RGB, respectively. Therefore, estimating RGBs causes the largest variations of the scale, whereas the simultaneous estimation of RGB and TRP corrections provides a more stable scale realization, and reduces the offset to 0.0 mm (Table 5). For the realization of ITRF2020, the long-term RGBs are estimated and then re-introduced to the SLR solutions to make the SLR scale closer to VLBI and other techniques of space geodesy. The approach for the ITRF2020 will lead to a similar scale realization as the estimation of RGBs in this study; however, the stability of scale estimates should be better than in the RGB solution because the range bias will be re-introduced. The estimated and re-introduced range biases change, however, the amplitudes of annual signals of some SLR-derived parameters, such as geocenter coordinates (Luceri et al. 2019). A smaller amplitude suggests that estimated biases could absorb some non-instrumental geophysical effects. This study shows that the co-estimation of TRP and RGB in SLR solutions does not harm the realization of the global scale and is entirely plausible as it provides more realistic scale realizations than the standard SLR solutions with neglecting the tropospheric biases. Figure 13 shows the spectral analysis of the reference frame scale factor. The RGB solution increases the semi-annual signal, which becomes larger than in other solutions, and introduces a signal of 149 days, which is not visible in other solutions. TRP and TRP + RGB solutions do not introduce any unwanted signals into the time series of the estimated scale; thus, the SLR-derived scale remains stable despite considerable correlations between the station heights, range biases, troposphere delays, and thus also the ITRF scale when all parameters are estimated together.

Geocenter coordinates
The observations of passive geodetic cannonball satellites, due to the low impact of non-gravitational perturbations on their orbits, are the most prominent targets for recovering the geocenter motion. We expect that the geocenter coordinates are not contaminated by any systematic errors due to biases. The ITRF realizations rely their origin uniquely on the SLR solutions; therefore, wrong or biased geocenter coordinates  would result in the biased coordinates of all techniques contributing to ITRF. Figure 14 shows the time series of the geocenter coordinates derived from SLR solutions with different handling of range and tropospheric biases smoothed by a low-pass Savitzky-Golay filter with 13-epoch windows. Table 6 shows that the STD solution has the smallest offsets for the Z and X geocenter components, which is expected, because the STD solution is most consistent with the modeling used for the ITRF2014 realization. The estimation of troposphere biases increases the long-term mean value by 0.1, − 0.6, and − 1.4 mm, for the X, Y, and Z geocenter components, respectively. The offset of the Z component in all non-standard solutions is different by more than 1 mm when compared to the STD solution with no biases considered.
The amplitudes of the annual signal in the geocenter motion are only marginally changed in solutions with different RGB and TRP handling. However, RGB and TRP biases can significantly change the mean geocenter offset by more than 1 mm for the Y and Z components. Unbiased geocenter has a fundamental meaning for the ITRF realization, whose origin solely relies on SLR. All three solutions, RGB, TRP, and TRP + RGB, seem to provide consistent offsets of the Z geocenter component, which is inconsistent with ITRF2014 and the standard ILRS solution. Also, RGB, TRP, and TRP + RGB solutions provide more internally consistent Y geocenter components than the standard solution. From this part, we conclude that uncalibrated biases substantially affect the geocenter coordinates, which are one of the fundamental products derived from the SLR technique, by more than 1 mm for the Y and Z components.

Impact of troposphere corrections and range biases on the second-degree Stokes coefficients
SLR is the fundamental technique for deriving Earth's oblateness term, C 20 , and other low-degree gravity field coefficients (Pearlman et al. 2019a;Sośnica et al. 2014;Blossfeld et al. 2015Blossfeld et al. , 2018Chen et al. 2021). In this part, we analyze the impact of handling of systematic errors on the estimated degree-2 Stokes coefficients derived from SLR. Figure 15 and Table 7 show the comparison of derived Stokes coefficients w.r.t. the STD solution. For all of analyzed coefficients, the solution TRP + RGB has the highest noise, which can be associated with the largest number of estimated parameters. The simultaneous estimation of many parameters burdens the time series by frequency noise especially in case of the S 21 and C 21 terms. The periodogram also shows that estimation of range biases absorbs some part of the annual signal of the term S 21 which does not occur in solution TRP (see Fig. 15). The term S 21 has a negative offset in RGB and TRP + RGB solutions, whereas the TRP solution agrees well with the standard solution.
The mean formal errors of LAGEOS-based C 20 are 4 ⋅ 10 −12 , whereas for C 21 , S 21 , C 22, S 22 the typical formal errors equal to 1 − 2 ⋅ 10 −11 in 7-day solutions. Table 7 shows that the RMS of difference w.r.t. the standard solution is at a comparable level to the formal errors of gravity field estimates. The mean offsets are smaller than the formal errors from 7-day solutions. However, when using the error propagation law, the mean errors of the offsets derived from 9 years of LAGEOS data are at the level of 2 ⋅ 10 −13 and 1 ⋅ 10 −12 for C 20 and other degree-two coefficients, respectively. Considering offsets exceeding the 2-sigma significance level, we can conclude that the mean offsets of C 20 and C 21 are insignificant; whereas the mean offsets of S 21 , C 22 , S 22 are significantly affected by RGB and TRP handling at two-sigma level. Therefore, the tropospheric and range biases have minor impact on the SLR-derived Earth's oblateness term; however, may significantly bias the estimation of other low-degree tesseral and sectorial gravity field coefficients.

Discussion and conclusions
Retrieving ZTD corrections from SLR observations is pragmatically possible and has a measurable impact on the fundamental products derived from SLR, such as the station coordinates, geocenter coordinates, and the global scale. Estimating troposphere delay corrections absorbs barometric-based pressure biases which occur at the stations due to a failure of meteorological devices, missing height corrections, or incorrect way of barometer installation. Estimating tropospheric delays absorbs also other distance-dependent biases, such as intensity biases, because of their similar nature to the tropospheric refraction and the dependency on the elevation angle. The combination of estimating range and tropospheric biases may absorb some other issues occurring at SLR stations related, e.g., to temperature biases or neglecting the temperature variations in the atmospheric profiles when considering only surface temperature values for the calculation of the tropospheric delays. However, the azimuthal asymmetry of the troposphere cannot be absorbed by RGB nor ZTD, thus, requires the proper handling of horizontal gradients. We conducted four tests with different handling of tropospheric delays and range biases for 9-year LAGEOS-1/2 solutions. Despite high correlation coefficients between the station heights, range biases, and zenith tropospheric delays, the estimation of all three types of parameters together turned out to be possible based even on sparse SLR observations. The a posteriori sigma of unit weight is at the same level in the TRP solution as in the RGB solution, although fewer parameters are estimated in the TRP solution than in the RGB solution, in which individual biases for each satellite are derived.
Estimating tropospheric corrections once per 7-day solution is plausible due to the physical properties of the optical wavelengths, and thus, the low sensitivity to the dynamically changing wet part of the troposphere delay. Such an approach would not be possible in space geodetic microwave-based techniques. In SLR solutions, we used the total delay based on in situ measurements as a priori delays and we estimated small residual tropospheric corrections. The temporal variations of the tropospheric delays are accounted for by in situ measurements, whereas the biases of meteorological instruments are captured by constant offsets that are estimated. Estimating a small number of tropospheric parameters guarantees that the high quality of SLR station heights is conserved. Thus, our approach is similar to that employed for microwave techniques, but not the same. The estimation of the troposphere offsets within 7-day intervals guarantees that the ZTD corrections are recovered from sparse SLR observations, whereas the determination of fundamental SLR parameters does not deteriorate.
The main effort of this paper was to introduce the separation of the systematic errors into two groups which has not been previously described in the literature. The first part of the two error types changes with a function of an elevation angle, whereas the second type of error incorporates biases that are constant irrespective of the elevation angle. Further analysis of systematic errors should consider satellites at different elevation angles and altitudes because, in such a way, the intensity biases can be separated from tropospheric biases and satellite signature effects which result from the detector-specific issues Couhert 2019;Strugarek et al. 2021). These three types of errors are inseparable when using satellites at almost the same altitude, such as LAGEOS-1 and LAGEOS-2. Otsubo et al. (2016) reported a significant dependency between average post-fit residual and the RMS of SLR normal points, especially for single-photon stations. RMS of normal points typically increases for low elevation angles due to atmosphere scintillations and increased refraction, whereas it reaches the minimum for zenith observations. The discovered dependency may thus have an origin in missing atmospheric correction or may be associated with the detector-specific issues causing a signature effect for Compensated Single-Photon Avalanche Diode (CSPAD) detectors. Both the atmospheric and instrumental effects can be mitigated by estimating elevation-dependent TRP bias.
We found that the estimation of RGB, in general, deteriorates the station coordinate repeatability, whereas the solution TRP does not deteriorate or even improves the inter-quartile ranges of station coordinates for most of the high-performing SLR stations. Neglecting the TRP or RGB corrections introduces biases in the estimated geocenter coordinates by about 1 mm in the case of the Z and Y components. The long-term mean value of gravity field coefficients S 21, C 22, and S 22 is significantly biased when using different handling of TRP and RGB. Moreover, the global scale is substantially shifted when adding the biases into the functional model. In the solutions with estimating only range biases, the station coordinates and the global scale parameter deteriorate due to the overestimation of the total effect by the range bias. It is much better when range biases are co-estimated together with tropospheric biases, as this solution provides much more stable station coordinates and scale estimates, whereas the geocenter coordinates are freed from systematic effects. From the comparison of RGB and TRP values derived from different solutions, we can conclude that most of the systematic effects in SLR observations are elevation-dependent, thus are better absorbed by the TRP bias than the RGB bias.