Troposphere delay modeling in SLR based on PMF, VMF3o, and meteorological data

Satellite laser ranging (SLR) requires accurate troposphere delay models to properly correct the observed distances to satellites and derive fundamental geodetic and geodynamic parameters. The currently used models for the tropo-spheric delay employ in situ meteorological data collected simultaneously with laser measurements. However, the standard models assume full symmetry of the atmosphere above the SLR stations because all meteorological data come from one sensor. In this study, we evaluate various methods of troposphere delay modeling based on numerical weather models, such as the Potsdam Mapping Function (PMF) and Vienna Mapping Function for optical frequencies (VMF3o), in situ measurements, the Wrocław gradient model (WGM) and the combination of different models. We found large discrepancies between pressure, temperature, and humidity records between in situ measurements and numerical models. The best results for the zenith delays are obtained when using in situ meteorological data with the estimation of tropospheric biases. For stations with some deficiencies in proper humidity measurements, e.g., Zimmerwald in Switzerland, the best results are obtained when using hydrostatic zenith delays based on in situ data and wet delays based on numerical weather models. Finally, we found that using horizontal gradients of the tropo-spheric delay is indispensable to avoid biases in the SLR-based Earth rotation parameters of approximately 20 µ as for polar motion. The horizontal gradients successfully account for the asymmetry of the troposphere above SLR stations and can be derived from PMF, VMF3o, or a parameterized WGM model with similar accuracies.


Introduction
The integrated space geodetic techniques contribute to the International Terrestrial Reference Frame (ITRF) (Altamimi et al. 2023) realizations by deriving crucial information about the origin, orientation, and geometric scale datum.The influence of the lowest layer of the atmosphere, the troposphere, is one of the major sources of errors in data processing of space geodetic technique observations.In recent years, we observed a vast number of studies related to tropospheric modeling for microwave-based geodetic techniques, such as Global Navigation Satellite Systems (GNSS) (Tregoning and Herring 2006) and Very Long Baseline Interferometry (VLBI) (Boehm et al. 2006;Zus et al. 2014;Landskron and Böhm 2018;Wang et al. 2022).Moreover, the tropospheric delay is also a subject of scientific interest in other techniques, such as satellite altimetry or Synthetic Aperture Radar Interferometry (InSAR), where local atmospheric variations significantly deteriorate interferograms (Wilgan et al. 2019).One of the main factors limiting the accuracy of standard Satellite Laser Ranging (SLR) solutions is neglecting some atmospheric effects, such as atmospheric asymmetry above the stations.The impact of horizontal refractivity on SLR was introduced by Abshire and Gardner (1985), who proposed implementing two-color laser systems at several stations located on different continents and tracking radiosondes equipped with retroreflectors and comparing the results to atmospheric models.Degnan (1993) reported that although measurements are conducted in a clear sky in many SLR observations, the emitted photons pass through the invisible parts of cirrus clouds, introducing unconsidered atmospheric delays.
The currently used tropospheric delay model in SLR consists of a common mapping function for the sum of hydrostatic and non-hydrostatic delays.The model is based on one year of radiosonde data (Mendes et al. 2002).The Mendes-Pavlis model for the mapping function (Mendes et al. 2002) requires a temperature record derived from the meteorological station at the epoch of observation.Hulley and Pavlis (2007) proposed the raytracing method to compute the slant total delay for each observation at the site.The simulation proposed by Hulley and Pavlis (2007) improved the residuals of observations, especially for low elevation angles at 10 degrees; additionally, the proposed model considers the atmospheric asymmetry above the station.Hulley and Pavlis (2007) demonstrated that the variance of SLR observation residuals can be reduced from 10-43% depending on the source of the numerical weather model (NWM) used.Wijaya and Brunner (2011) derived a formula that can be used to determine the total troposphere delay for stations equipped with two laser frequencies (2F-SLR) by considering the impact of the non-hydrostatic delay.In that case, the water vapor pressure could be delivered from the GNSS station co-located with the SLR station.However, the accuracy of current photodetectors is beyond the range to render this formula useful in the operational mode in two-color ranging.The further proposals of troposphere delay models, mapping functions, and horizontal gradient models for SLR are based on NWM.These models are called the Potsdam Mapping Function (PMF) (Zus et al. 2015) and the Vienna Mapping Function 3 for optical frequencies (VMF3o) (Boisits et al. 2020) and employ ERA5 data and operational products of the European Center of Medium Weather Forecast (ECMWF) for PMF and VMF3o, respectively.However, none of the above-mentioned atmospheric models can fully handle the atmospheric turbulence that occurs above the stations.In PMF, due to the low sensitivity of SLR observations to the wet part of the troposphere delay, the separation of the mapping function for laser frequencies is neglected.This approach is consistent with the statement from the International Earth Rotation and Reference Systems Service (IERS) Conventions 2010 (Petit and Luzum 2010) concerning the delays of laser techniques in the troposphere.In contrast, VMF3o proposed by Boisits (2018) introduced the concept of the separated mapping function to the hydrostatic and nonhydrostatic model of tropospheric delay, as well as separated zenith total delay (ZTD) and horizontal gradients.The VMF3o approach allows additional experiments to be conducted with a mixed approach of zenith wet delay (ZWD).In almost all of the analyzed cases in this paper, we use the zenith hydrostatic delay (ZHD) and ZWD based on meteorological measurements derived from the co-located meteorological stations at the SLR site.Drożdżewski et al. (2019) noticed that the Mendes-Pavlis model for the ZTD and ZWD extended by PMFbased horizontal gradients improves the consistency between SLR and microwave-based techniques (GNSS, VLBI) by reducing the mean long-term offset at the level of 20 µ as for the X and Y components of the polar motion.Drożdżewski and Sośnica (2021) introduced the troposphere delay correction in SLR in the form of estimated tropospheric biases.The authors also mentioned that other errors, which have similar nature to the troposphere delay correction, such as intensity biases, are also absorbed by estimated tropospheric biases.The authors described the impact of tropospheric biases, such as station coordinates, geocenter coordinates, Earth rotation parameters, and low-degree spherical harmonics of Earth's gravity field based on SLR observations to LAGEOS satellites.Strugarek et al. (2022) confirmed that the estimation of troposphere correction removes the skewness from the observation residuals when processing SLR data to low Earth orbiters.Strugarek et al. (2022) confirmed that using the estimation of troposphere delay removes other elevation-dependent errors.However, one should be careful with estimating tropospheric corrections due to strong correlations with other parameters, such as horizontal components of station coordinates and horizontal gradients.Moreover, for a limited number of SLR observations, constraining tropospheric parameters is indispensable.In this paper, we provide a comparison of the recently developed tropospheric models for SLR, i.e., PMF and VMF3o, and the impact of the mapping function differences, horizontal gradients, and ZTD.These models are shown first as a comparison of input data, and second, as the impact on the crucial geodetic parameters, such as station coordinates and geocenter coordinates, on an example of SLR observations to passive geodetic satellites LAGEOS 1 and LAGEOS 2. PMF and VMF3o are compared to the standard Mendes-Pavlis model that employs meteorological records for tropospheric delays.The corrections to the tropospheric delays are estimated for each model to verify the level of tropospheric biases.Finally, the mixed approach is tested with selected data from meteorological sensors and some from NWM.

Methods/experimental
The currently recommended IERS mapping function FCULa (Mendes et al. 2002) is based on one year of radiosonde data from 1999 and in situ temperature, pressure, and humidity records derived from meteorological stations.The new approaches to determining mapping functions for SLR, such as PMF and VMF3o, are based only on NWM ERA5 and ECMWF operational forecasts, respectively.Table 1 provides the selected differences between the ERA5 and ECMWF operational models (Hersbach et al. 2020).Drożdżewski et al. (2019) showed the benefits of applying PMF based on the ERAinterim reanalysis.The new release of ECMWF reanalysis the so-called ERA5 distinguishes better spatial and temporal resolution with respect to its predecessor, ERA-interim.Both NWMs are also based on in situ measurements of temperature, pressure, and humidity.Therefore, the mixed approach of mapping function and gradients based on NWM and ZTD based on in situ data, as tested in this paper, should still be consistent.Huang et al. (2023) emphasize that in ERA 5, the annual pressure biases are reaching 0.07 hPa on the surface of 31 km squared.McGarry et al. (2019) defined that the in situ atmospheric pressure measured in NASA's Next-Generation Satellite Laser Ranging (NGSLR) system should measure the pressure with accuracy at the level of 0.08 hPa.However, in tests conducted by Drożdżewski et al. (2019) the ZTD derived from NWM resulted in worse SLR results than when using in situ pressure data.
Even though both models, PMF and VMF3o, are based on the NWM meteorological analysis, the reader should be aware of some crucial differences.VMF3o proposes a mapping function dedicated to laser observations based on the VMF3 concept (Landskron and Böhm 2018) which splits the hydrostatic and non-hydrostatic part for: (1) ZTD, (2) the mapping function, and (3) horizontal gradients.Splitting the ZTD and mapping function is commonly used in GNSS.Due to differences in the physical properties of the microwave and optical signals, dedicated mapping functions for SLR are expected (see Table 2).The PMF model assumes that due to low sensitivity to the wet part of the troposphere delay, the separation of the mapping function for laser frequencies could be neglected for observations provided at an elevation angle above 10 degrees.This concept is consistent with the IERS 2010 Conventions that recommend the Mendes et al. (2002) model for the mapping function which combines the hydrostatic and wet parts of the delay.In PMF, all coefficients are tabularized, whereas VMF3o uses the expansion into spherical harmonics for selected coefficients.Finally, PMF applies nonlinear gradients using 7 coefficients to describe the atmospheric asymmetry, whereas VMF3o provides 2 linear gradient coefficients for each hydrostatic and wet delay.

Standard solution based on a common mapping function for ZTD
The standard approach for SLR tropospheric delay modeling (Mendes and Pavlis 2004)

PMF solution
The PMF in the original version included 2 linear coefficients for the gradients and 3 additional coefficients for nonlinear dependencies (Drożdżewski et al. 2019).The updated formula of PMF is consistent with that of Zus et al. (2021), who extended the gradient model by highorder terms.The PMF dedicated to SLR is adopted for SLR frequencies equal to 532 nm and reads as: where m PMF (e) is the mapping function, A describes the azimuth of each observation, G N and G E represent the horizontal gradients for the north and east components, respectively, and GZ1-Z4 represent the higher-order horizontal gradient terms.The mapping function for the gradients considers both the hydrostatic and wet part of the delay with k = 0.0032 (Chen and Herring 1997) being consistent with the dominating hydrostatic part:
(5) d atm =d atmh + d atmw . (6) where d atm h describes the hydrostatic delay consider- ing horizontal refractivity projected on elevation angle e, and G Nh and G Eh describe the hydrostatic part of hori- zontal gradients for the north and east components.The a (h,w) , b (h,w) c (h,w) describe the mapping function coeffi- cient for the hydrostatic part (h) and for the wet part (w).G N w, G E w represent the wet part of the horizontal gra- dients for the north and east components.The m gh , m gw are functions for projecting horizontal gradients, where k h is equal to 0.0031 and k w is equal to 0.0007 (Chen and Herring 1997).In the other approaches without separation into hydrostatic and non-hydrostatic parts, the same function m g is typically used with a k coefficient equal to 0.0032.VMF3o provides tabular values for the mapping function: a wet and hydrostatic coefficients and zenith delays d h and d w , whereas b and c are provided as an expansion into spherical harmonics separately for the mean values and seasonal signals and with a separation between the hydrostatic and wet parts.Moreover, VMF3o provides tropospheric data from NWM, such as temperature, pressure, and water vapor, interpolated for the location of the SLR station, allowing for direct comparison with in situ data.The latency of the operational VMF3o products is 24 h; the new data are available at approximately 18:00 UTC every day for the previous day.(11) m gw (e) = 1 sin(e) • tan(e) + k w .
Predictions for the next day are generated at 9:00 UTC; however, they are not publicly available on the VMF website.Therefore, VMF3o can be applied for operational SLR products, as the latency of the product is approximately one week.

WGM solution
The Wrocław Gradient Model (WGM), derived in this study, is based on time series of linear gradients derived from NWM fitted into an empirical function.The WGM employs the standard FCULa mapping function derived by Mendes et al. (2002).For gradients, WGM is based on a fit of the first-degree polynomial with the annual and semi-annual signals into 16 years of data for the east and north components of the total delay, i.e., the sum of the hydrostatic and wet parts.This solution allows for the reconstruction of most of the major signals included in the horizontal gradients.Figure 1 shows the stacked amplitudes for the north component of 43 SLR stations.
For most stations, characteristic annual and semi-annual signals are observed, as well as daily and monthly fluctuations, which were described by Boisits et al. (2020).Nevertheless, the daily and monthly fluctuations have variable and small amplitudes; therefore, annual and semi-annual signals with stable amplitudes were selected for the simplified models of gradients.
where a 0 describes the offset, a 1 is the drift, a s1 and a c1 are the annual components, a c2 , a s2 describe the semi- annual components, t describes the elapsed days from the beginning of a year and T is the length of the year in days.Finally, the slant delay can be estimated as a sum of the ZTD multiplied by a mapping function and gradients multiplied by the gradient mapping function developed by Chen and Herring (1997): To determine this model, the discrete gradients dedicated to SLR optical observations from both PMF and VMF3o sources can be used with no significant differences. (12) (13) However, for WGM derivation in this study, PMF was employed.The WGM can also be applied for both operational and postprocessing products.The polynomial consists of the annual and semiannual signals that retrieve only the systematic, deterministic part of the gradients (see Fig. 2).The mean value a 0 for the north component of gradients has a negative value in the Northern Hemisphere, with the maximum value at the level of − 0.5 mm for Japanese stations (see Fig. 3).The mean value of the north gradient in the Southern Hemisphere is positive.

Results and discussion
In this section, we analyze the impact of different approaches to troposphere delay modeling in mixed configurations with respect to the standard approach recommended by the 2010 IERS Conventions (Petit and Luzum 2010).Figure 4  Figure 5 shows the horizontal gradients derived from PMF and VMF3o projected at 10 degrees elevation angle.We observe the highest differences for station Mt Stromlo at the level of 3 mm.In general, we can perceive very good agreements between PMF and VMF3o solutions, which confirm that the separation of mapping function in SLR measurements can be neglected.

Meteorological data records comparison
The major difference between PMF, VMF3o, and FCULa is the source of input data.PMF and VMF3o are based on the pressure, temperature and humidity derived from NWM, but only VMF3o products make publicly available the discrete values of temperature, pressure, and water vapor pressure with a 6-hour time span (Boisits et al. 2020).Hence, we decided to compare the in situ meteorological data with discrete data derived by VMF3o interpolated to the epochs of actual SLR observations.Figure 6 shows differences in the time series of pressure records for selected stations.Drożdżewski and Sośnica (2021), introduced the estimation of troposphere delay corrections as a method to detect and handle barometer malfunctioning.The estimated ZTD values can detect pressure jumps that occurred recently, e.g., at the Wettzell and Graz stations.Boisits et al. (2020) showed the differences between in situ pressure and temperature measurements and NWM.However, the question of these inconsistencies and the application of this information by the SLR community still remains open.Although the pressure derived from NWM is still beyond the accuracy of the in situ pressure measurements, the analysis of pressure differences with an assessment of estimated station coordinates may help expose the undetected and unmodeled pressure jumps in the Inernational Laser Ranging Service (ILRS) data time series.The SLR products are very sensitive to pressure biases; as a result, undetected pressure issues may deteriorate crucial geodetic parameters, such as station coordinates, the scale of the reference frame, or geocenter coordinates (Drożdżewski and Sośnica 2021).
Figures 7 and 8 show the comparison of in situ temperature and humidity measurements and VMF3o values for selected stations with an approximate separation of the seasons of the year.For station Haleakala (7119), we observe a correlation between in situ measurements and VMF3o of 0.58 and an interquartile range (IQR) of differences between records at the level of 3.6 • C (see Table 3).This value indicates a rather poor agreement and a substantial discrepancy between in situ measurements and VMF3o, which will be investigated in the further part of this article.Station San Juan (7406) is located in Argentina at a height of 727 m above sea level.For this SLR site, the lowest temperature is noticed in the middle of the year (red dots in Fig. 7); for that period, the highest discrepancies can be noticed.Figure 8 illustrates the comparison of water vapor pressure (WVP) between in situ and VMF3o records.The WVP impact in SLR Fig. 3 The mean offset of the WGM gradient model (mm).Topnorth component, bottom-east component measurements contributes to approximately 1% of the ZTD values, which is approximately equal to subcentimeter values in zenith.However, even several millimeters of the discrepancies in the zenith direction may result in a few centimeters of errors for observations projected at low elevation angles.For stations Matera (7941) and Zimmerwald (7810), the in situ measurements agree quite well with VMF3o in the winter, spring, and autumn, whereas large discrepancies occur in the summer, where VMF3o is remarkably lower than the in situ records, even by one third.

Troposphere delay bias as a tool for mapping function validation
We test different approaches to zenith delay modeling with and without separation between wet and hydrostatic delays for zenith delays and mapping functions, considering horizontal gradients from VMF3o, PMF, and simplified model WGM, as well as mixed approaches, in which some information is taken from NWM and some is based on in situ measurements.The tested solutions are shown in Table 4.Moreover, ZTD corrections to a priori values are also estimated in selected solutions to identify tropospheric biases (TRP).The different approaches to modeling the troposphere delay in SLR may be assessed by the estimated troposphere delay corrections for each station (Landskron and Böhm 2018).The same constraint for ZTD correction equal to 0.01 m with respect to the a priori value was applied in all solutions.The first two solutions from Table 4, Mendes-Pavlis (M-P), correspond to the approach currently used by the ILRS in operational SLR products, except for the ZTD bias estimation in the second solution (M-P TRP).In the ILRS approach, all tropospheric delay components are derived from in situ measurements.In other solutions from Table 4, some components were replaced by corresponding values based on NWM, e.g., ZWD or ZHD.
Figure 9 shows the box plots of estimated ZTD corrections from solutions using different a priori delays and modeling.For some stations, the estimated corrections agree between solutions using in situ and NWM data; for some other stations, the corrections are different, especially for those where atmospheric pressure records are employed from NWM.For station Simosato (7838), the models based on atmospheric pressure derived from NWM, i.e., VMF3o ZTD and PMF ZTD, lead to a positive estimated ZTD bias, as opposed to all other approaches, in which ZHD is based on in situ atmospheric pressure records (M-P, VMF3o ZWD, VMF3o, M-P model).The atmospheric pressure records from NWM can thus be biased because solutions using NWM pressure lead to a much larger bias for Simosato than when based on in situ pressure data.The pressure differences together with estimated troposphere delay correction are reliable methods to identify and handle such kinds of pressure errors.For Greenbelt and Concepcion, the solutions based on NWM lead to large negative tropospheric biases of several millimeters.For Wettzell, most of the solutions provide similar results, suggesting that the in situ atmospheric pressure measurements are incorrect.In Concepcion, using in situ atmospheric pressure measurements Fig. 5 Comparison of PMF and VMF3o Horizontal gradients projected onto 10 • elevation angle for selected stations: Herstmonceux, Mt Stromlo, Graz and Yarragadee.O1 denotes the linear gradients from PMF; O2 denotes the higher-order gradient terms in PMF leads to a positive bias of approximately 2.5 mm, whereas the NWM pressure values lead to estimates of large negative biases out of range of Fig. 9.This implies that both NWM and in situ data are inconsistent with actual tropospheric states in Concepcion.For Potsdam (7841) and Yarragadee (7090), differences between estimated tropospheric biases are below 1 mm for different approaches.Although PMF and VMF3o are based on the ECMWF products, the results slightly differ from each other for the same stations.For Potsdam (7841), Grasse (7845), and Tahiti (7124), the median value of tropospheric bias is the smallest in VMF3o ZTD, which uses NWM for the a priori atmospheric pressure.The complete PMF model (PMF ZTD) provides the smallest biases for Hartebesthoek (7501), Wettzell (8834), and Yarragadee (7090); however, for the mean biases, PMF-based biases do not change much for these stations.In the case of utilizing the wet part of the troposphere delay from NWM (VMF3o ZWD), a ZTD bias reduction at the level of 0.2 mm is noticed for station Matera (7941).In the case of solutions PMF and VMF3o, the small estimated ZTD bias is related to improvement emerging from using the mapping function based on NWM data instead of in situ measurements.Using PMF leads to reduced biases for Tahiti (7124), Greenbelt (7105), Matera (7941), and Zimmerwald (7810).The mean bias reduction introduced by VMF3o is at the levels of 0.1, 0.2, 0.1, and 0.2 mm in the zenith direction for stations Graz (7839), Potsdam (7841), and Matera (7941), respectively.For stations Zimmerwald (7810) and Mt Stromlo (7825), the reduction of the interquartile range (IQR) in solution VMF3o ZWD has been noticed as well as a reduction of the range of the whiskers.For Zimmerwald, discrepancies between in situ measurements and NWM were detected in the summertime (see Fig. 8).The smallest IQR values from the VMF3o ZWD TRP solutions mean that the best solution can be obtained when the hydrostatic zenith delay is based on in situ data, whereas the wet delay is derived on the basis of NWM for Zimmerwald.Therefore, the unique feature of VMF3o with the separation between the wet and hydrostatic delay may be useful to account for the errors of WVP recording.For station Matera (7941), where the inconsistency between WVP derived from in situ measurements and VMF3o has been noticed, we do not observe any improvements in the solutions with ZHD based on NWM.

Station coordinate repeatability
Figure 10 illustrates the 3d RMS of station coordinates' repeatability based on observations from the period 2010.0-2019.0 for selected representative solutions.Only the M-P solution does not include the estimation of tropospheric corrections, whereas all other solutions use ZTD values from different sources with zenith troposphere corrections estimated.Figure 10 shows that the difference in the station coordinate repeatability is minor irrespective of the a priori ZTD values used when estimating tropospheric corrections.Figure 11 shows the difference in station coordinate repeatability with respect to the M-P solution with no troposphere corrections estimated for the north, east, and up components.The horizontal gradients from NWM or a simple WGM model especially improve the horizontal components of station coordinates, which is visible for the north component.The estimated ZTD correction mainly impacts the up component of the station coordinates.Improvements in station coordinate repeatability above 1 mm were observed for stations Shanghai (7821), Simosato (7838), and Changchun (7237).Improvement at the level of half mm was noticed for stations Haleakala (7119), Grasse (7845), and Zimmerwald (7810).The deterioration of average values of station coordinate repeatability has been noticed for stations Monument Peak (7110) and Yarragadee (7090).For station Shanghai (7821), we observe an improvement exceeding 1 mm.For stations Changchun (7237) and Simosato (7838), the improvement is at the level of 2 mm even for solutions with ZTD derived from NWM w.r.t. the standard M-P solution.The improvement of horizontal components is mainly related to considering horizontal gradients.Drożdżewski et al. (2019) noticed that neglecting horizontal gradients leads to the deterioration of the consistency between SLR and IERS-C04-14 Earth Rotation Parameters (ERPs).Indeed, Table 6 shows that the horizontal gradients reduce the mean offsets for the X and Y pole coordinate components at the level of 20 µ as between the SLR solutions with and without gradients and the IERS-C04 series.The horizontal gradients can be taken from a simplified WGM parametrized model of horizontal gradients (M-P model) or NWM (PMF or VMF3o).In the case in which the ZTD model is used from NWM without estimating ZTD corrections, the X and Y pole coordinates deteriorate in terms of offset between the IERS-C04-14 combined solution and mean sigma (see Table 5).When the ZTD from NWM is used, the troposphere delay correction has to be estimated; otherwise, the ZTD derived only from NWM introduces deterioration of the X and Y ERP components at the levels of 22 and 8 µas, respectively (see Table 5).In solutions, VMF3o ZTD TRP and PMF ZTD TRP, the length-of-day (LOD) parameter deteriorates at the level of 18 and 20 µ s/day for solutions PMF ZTD no est.and VMF3o ZTD no est., respectively.Thus, the estimated troposphere delay correction is not able to fully account for the large sigma and offsets observed in pole coordinates when the a priori ZTD values are biased.The aim of testing the solution VMF3o ZWD TRP is to show the impact of the separation mapping function on the crucial geodetic parameters as well as to validate the SLR in situ humidity records.However, the solution VMF3o ZWD TRP is consistent with the solution PMF TRP, where the separation of the mapping function is neglected.The most consistent solution with the IERS-C04-14 product is the solution PMF TRP.Table 5 does not indicate one solution with superior performance with respect to the other solutions.However, the best solutions are obtained when (1) horizontal gradients are considered, (2) ZHD is based on in situ meteorological data, and (3) tropospheric correction is estimated.The selection of the mapping function and the separation between wet and hydrostatic delays play a minor role.The set of best-performing models includes the VMF3o TRP, VMF3o ZWD TRP, PMF TRP, and M-P model TRP.Therefore, it is sufficient to consider tropospheric biases and gradients from the WGM to reduce the mean pole offsets by approximately 20 µas.

Geocenter coordinates
We compare the solutions with horizontal gradients, the mapping function derived from NWM (PMF or VMF3o) and the ZTD correction with the M-P standard ILRS solution for the geocenter coordinates.Next, the solutions with ZTD derived from NWM with the estimation of the troposphere delay correction are compared to the standard M-P solution and to the mixed solutions: PMF TRP and VMF3o TRP. Figure 12 shows the differences in the estimated geocenter coordinates from the standard M-P solutions.Figure 12 (top) reveals large differences and fluctuations of 1 mm at the end or beginning of the year, which is probably related to the lower number of observations and poor observation geometry that becomes sensitive to individual SLR stations contributing to the solution.However, the Z component of the geocenter coordinates reflects substantial differences at the end of the year at the level of 5 mm. Figure 12 (bottom) shows the Fast Fourier Transformation (FFT) results of the geocenter differences, confirming the moderate consistency of the annual signal for all of the analyzed solutions at the submillimeter level.
Figure 13 shows the comparisons of differences in geocenter coordinates for solutions in which the ZTDs are derived from NWM with respect to the standard M-P solution.Although the troposphere delay correction is estimated in solutions PMF TRP and VMF3o TRP, highfrequency noise can be noticed (Fig. 13, bottom), and artificial signals occur.The annual signal of geocenter coordinates is consistent at the level of 0.3 mm with other analyzed solutions (see Table 6).SLR observations provide stable information about the origin of the reference frame that coincides with the long-term mean value of the geocenter.Table 6 shows, however, that the geocenter coordinates are mostly sensitive to the estimation of tropospheric corrections.The offset of the Z component changes between − 1.7 and − 2.4 mm in the solutions with different handling of tropospheric biases, i.e., with and without estimating TRP biases.Thus, the estimation of the tropospheric biases plays a main role in the recovery of the geocenter motion, whereas a priori ZTD values, horizontal gradients and mapping functions play a minor role without any substantial offsets.

Conclusions
The comparison of meteorological parameters (pressure, temperature, humidity) derived from in situ measurements and VMF3o pressure records together with troposphere delay corrections as an indicator of biases demonstrates that ZTD corrections should be estimated to properly handle biases in meteorological data.The estimation of troposphere delay corrections slightly increases the offsets between the ERPs derived from SLR and the combined IERS-C04 solution.In the case  of solutions with the troposphere delay model based on NWM, the mean offsets are absorbed by troposphere delay correction, and the ZTD does not deteriorate the consistency between SLR results and the IERS-C04 series.The ZTD model derived from the NWM is not the most suitable approach for SLR data processing when used as an a priori model without estimating corrections.However, the mixed approach of mapping function based on NWM and zenith delays based on in situ data may mitigate some malfunctioning of in situ records.
The troposphere delay correction reduces the deterioration of station coordinate estimates and the mean offsets in pole coordinates to some extent in the case of solutions in which the ZTD is based only on the NWM.The m gh (e) = 1 sin(e) • tan(e) + k h .
For the east component, the values of the mean offset are in the range from − 0.15 to 0.15 mm; only for a few stations located in South America the values are at the level of ±0.3 mm.In central Asia and South America, larger fluctuations in the mean east component occur due to local atmospheric conditions and orography.

Fig. 1 Fig. 2
Fig. 1 Stacked amplitudes of horizontal gradients for the north component of 43 SLR stattions (mm).The x-axis is scaled in days shows the differences in delays projected at low elevation angles of 7, 15, and 20 degrees for two high-performing SLR stations located in the Northern and Southern hemispheres, Wettzell (8834) and Yarragadee (7090).The top figure illustrates the differences between the PMF and Mendes-Pavlis (MP) mapping function (FCULa), the middle figure shows the differences between VMF3o and Mendes-Pavlis, and finally, the bottom figure shows the differences between PMF and VMF3o.The largest differences are obtained for observations provided at an elevation angle of 7 degrees.The PMF and VMF3o, which are both based on NWM, are consistent even for 7 degrees of elevation angle at the level of 1-2 mm.In the case of the comparison with the Mendes-Pavlis model, the mapping functions generate differences at the level of ± 8 mm.

Fig. 4
Fig. 4 Differences between mapping function values in 2008 for stations Wettzell and Yarragadee located in the Northern and Southern Hemispheres, respectively.The differences are projected at the levels of 7, 15, and 20 degrees and shown in red, green, and blue, respectively

Fig. 6
Fig. 6 Differences between pressure records derived from in situ records at SLR stations and VMF3o based on NWM with discontinuities detected in the period from 2008 to 2019 for Grasse, France (7845), Simosato, Japan (7838), Wettzell, Germany (8834) and Matera, Italy (7941).The green and red solid line represent the annual signal fitted by least squares adjustment

Fig. 7
Fig. 7 Comparison of temperature from in situ measurements at SLR stations and VMF3o interpolated for the SLR sites in years 2008-2012 in Celsius degrees as a function of the day of the year (color scale: winter-blue, spring-green, summer-red, fall-yellow in the Northern hemisphere)

Fig. 8
Fig. 8 Comparison of water vapor pressure (WVP) from in situ measurements at SLR stations and VMF3o interpolated for the SLR sites in 2008-2012 in hPa as a function of the day of the year (color scale: winter-blue, spring-green, summer-red, fall-yellow in the Northern Hemisphere)

Fig. 9
Fig. 9 Box plots of troposphere delay correction derived from solutions using different approaches to troposphere delay modeling based on in situ and NWM data, according to Table 4.The black box describes the 25th and 75th percentiles, the magenta central line describes the median value, and the whiskers describe the most extreme data points without outliers.Stations are sorted by the number of observations

Fig. 11
Fig. 11 Differences in station coordinate repeatability for the north, east, and up components of station coordinates in the 2010.0-2019.0period.The positive values denote the deterioration of station coordinate repeatability with respect to the M-P solution, whereas the negative values denote improvements

Fig. 12
Fig. 12 Time series of differences in geocenter coordinates for solutions with ZTD based only on the standard ILRS solution (top) w.r.t. the reference solution (M-P).and periodogram of geocenter coordinates (bottom) for the X, Y and Z components, respectively.The solid line represents the Savitzky-Golay filter with 11 weeks of windowing

Fig. 13
Fig. 13 Time series of differences in geocenter coordinates for solutions with different approaches of ZTD modeling (Top) w.r.t. the M-P solution.The solid line represents the Savitzky-Golay filter with 11 weeks of windowing.The bottom figure describes the periodogram of geocenter coordinates for the X, Y, and Z components

Table 2
Characteristics of signal propagation in the atmosphere for GNSS and SLR solutions with respect to the IERS 2010 Conventions Based on Mendes et al. (2002)-radiosonde data Troposphere asymmetry modeling Estimation of horizontal gradients, possibly with a priori values derived from NWM Not considered Geometry of observations Good geometry and dense distribution of observations A limited number of observations Ionosphere Dispersive medium No ionosphere delay pressure, temperature, and humidity resulting in the slant delay ( d atm ):

Table 3
Statistics of temperature and humidity differences between in situ measurements and VMF3o for selected SLR stations

Table 4
Description of analyzed solutionsTRP denotes that the trosphere bias is estimated.M-P model denotes using the a prori WGM for horizontal gradients.Temp.denotes temperature

Table 5
Difference in pole coordinates with respect to the IERS-C04-14 series

Table 6
Mean values of geocenter coordinates with respect to the SLRF2014 and the amplitude of the annual signal