On the relation of GNSS phase center offsets and the terrestrial reference frame scale: a semi-analytical analysis

Phase center offsets (PCOs) of global navigation satellites systems (GNSS) transmit antennas along the boresight axis introduce line-of-sight-dependent range changes in the modeling of GNSS observations that are strongly correlated with the estimated station heights. As a consequence, changes in the adopted PCOs impact the scale of GNSS-based realizations of the terrestrial reference frame (TRF). Vice versa, changes in the adopted TRF scale require corrections to the GNSS transmit antenna PCOs for consistent observation modeling. Early studies have determined an approximate value of α=-0.050\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha =-0.050$$\end{document} for the ratio of station height changes and satellite PCO changes in GPS orbit determination and phase center adjustment. However, this is mainly an empirical value and limited information is available on the actual PCO-scale relation and how it is influenced by other factors. In view of the recurring need to adjust the IGS antenna models to new ITRF scales, a semi-analytical model is developed to determine values of α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} for the four current GNSSs from first principles without a need for actual network data processing. Given the close coupling of satellite boresight angle and station zenith angle, satellite PCO changes are essentially compensated by a combination of station height, zenith troposphere delay, and receiver clock offset. As such, the value of α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} depends not only on the orbital altitude of the considered GNSS but also on the elevation-dependent distribution of GNSS observations and their weighting, as well as the elevation mask angle and the tropospheric mapping function. Based on the model, representative values of αGPS=-0.051\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _\text {GPS}=-0.051$$\end{document}, αGLO=-0.055\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _\text {GLO}=-0.055$$\end{document}, αGAL=-0.041\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _\text {GAL}=-0.041$$\end{document}, and αBDS-3=-0.046\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _\text {BDS-3}=-0.046$$\end{document} are derived for GPS, GLONASS, Galileo, and BeiDou-3 at a 10∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^\circ $$\end{document} elevation cutoff angle. These values may vary by Δα≈0.003\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta \alpha \approx 0.003$$\end{document} depending on the specific model assumptions and data processing parameters in a precise orbit determination or precise point positioning. Likewise changes of about ±0.003\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm 0.003$$\end{document} can be observed when varying the cutoff angle between 5∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5^\circ $$\end{document} and 15∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$15^\circ $$\end{document}.


Introduction
Compared to other space geodetic techniques, terrestrial reference frames (TRFs) based on global navigation satellites systems (GNSSs) benefit from a particularly dense network, which facilitates TRF access to common users. On the other hand, the technique still suffers from limitations in the independent realization of the TRF scale. As a one-way measurement system, GNSS station heights are strongly correlated with the receiver clock and other parameters. Most notably, the long-standing lack of absolute phase center calibrations for GNSS satellite transmit antennas has, until recently (Villiger et al. 2020;Zajdel et al. 2022), inhibited the use of GNSS observations for TRF scale determination. Instead, phase center offsets (PCOs) of GNSS satellite antennas had to be estimated from observations of a global network to align GNSS station heights with the TRF scale derived from Very Long Baseline Interferometry, VLBI, and Satellite Laser Ranging, SLR (Schmid et al. 2007. Among the first studies on this subject, Springer (1999) describes selected experiments to assess the relation of GPS PCO and scale changes in precise orbit determination (POD). For a 1 m PCO change, a −53 mm shift in station height, corresponding to a −8.3 ppb change in the TRF scale, was obtained along with a 5 mm change in tropospheric zenith delay. In a related effort, the impact of GNSS PCOs on scale changes in POD and in network solutions based on precise point positioning (PPP) is assessed in Zhu et al. (2003). Using a network of 42 stations, and observations to 24 GPS satellites with an elevation mask of 15 • , the authors demonstrate a −51 mm change in station height h, when increasing the z PCO component of all satellites by a common value of 1 m. The same ratio α = Δh/Δz PCO was obtained when adjusting the GPS orbits as well as working with fixed orbits. This reflects the fact that the orbital radius in GNSS POD solutions is tightly constrained by the dynamical models and largely insensitive to the assumed phase center location.
The coupling between GNSS satellite PCOs, tropospheric zenith delays, station clocks offsets, and the terrestrial reference scale was further analyzed by Rebischung (2014), who studied the correlation of these parameters in quantitative terms based on dedicated simulations. Among others, the author highlights the adverse impact of high elevation cutoff angles and a down-weighting of low-elevation observations in the POD process, which may result in an almost full correlation of those parameters. With respect to multi-GNSS solutions, a decreased correlation is identified when jointly processing two constellations with different orbital radius due the use of a common station height and troposphere zenith delay parameter.
Building on the results of Zhu et al. (2003), Ge et al. (2005) assessed the impact of satellite block-specific PCO errors and showed that the overall TRF scale varies with the fraction of individual block types within the constellation. Complementary to these studies that focused on the global scale impact of satellite PCOs, Cardellach et al. (2007) examined geographically dependent errors in the position and velocity fields of global networks caused by PCO errors. Using a simulation approach based on a simplified observation model, the authors demonstrate the presence of vertical network distortions at the level of 1% of the PCO uncertainty for the GPS constellation.
As a consequence of the aforementioned limitations, satellite-specific PCOs have repeatedly been estimated by the International GNSS Service (IGS; Johnston et al. 2017) to achieve consistency of IGS products with specific TRF versions (Ray et al. 2013). For GPS, Schmid et al. (2007) present PCO estimates adjusted to match the scale of the IGb00 reference frame, which is itself aligned to the year 2000 release of the International Terrestrial Reference Frame (ITRF2000; Altamimi et al. 2002). Later, these values were complemented with ITRF2005 (Altamimi et al. 2007) compatible GLONASS-M PCO estimates derived by Dilssner et al. (2010) and Dach et al. (2011). Following the release of ITRF2008, consistent PCOs for GPS and GLONASS were adjusted in Schmid et al. (2016) and a new IGS antenna model, igs08.atx, compatible with the respective scale was introduced for IGS processing. Latest updates for igs14.atx and IGb14 are discussed in Rebischung and Schmid (2016) and Rebischung (2020).
With the public release of manufacturer calibrations for the PCOs of QZSS (Cabinet Office 2022), Galileo (GSC 2022), BeiDou (CSNO 2019), and GPS III satellites (USCG 2022), a new interest has been raised to study the potential contributions of GNSS to the TRF scale and to understand the relation between PCO and station height or scale changes. Villiger et al. (2020) study the impact of manufacturer calibrated Galileo transmit antenna PCOs and demonstrate a 6.4 mm offset in station height, i.e., a 1 ppb scale difference, compared to the scale of ITRF2014 (Altamimi et al. 2016), when working with chamber-based receiver antenna calibrations for the Galileo frequencies. In a similar context, Xia et al. (2020) and Zajdel et al. (2022) analyze the consistency of manufacturer calibrated BeiDou-3 transmit antennas with the ITRF2014 reference frame for different signal combinations. More specifically, Zajdel et al. (2022) obtain a ratio of −0.052 between estimated station heights and PCO shifts of the BeiDou-3 MEO satellites, which is close to and even slightly larger by magnitude than the value derived by Zhu et al. (2003) for GPS.
The most recent ITRF2020 reference frame release offers improved consistency of the SLR and VLBI scales with a difference of 0.15 ppb as compared to a 1.38 ppb offset between the SLR and VLBI contributions in ITRF2014 (Altamimi et al. 2022;Pavlis et al. 2022). However, the scale changes with respect to ITR2014 necessitate a renewed adjustment of GNSS transmit antenna PCOs for the IGS antenna model. For a consistent translation of scale changes to PCO changes for all current GNSSs, it is important to know the respective ratios α = Δh/Δz PCO for each individual constellation and to understand the respective uncertainties.
With this background, the present study aims to investigate the mutual relation of station heights and satellite PCOs in GNSS POD and PPP as well as the various parameters that affect their ratio for each GNSS. A tailored observation model is used for this purpose that enables a semi-analytical computation of α without a need of performing full-scale network solutions with real observations. Constellation-specific parameters of interest in this framework are summarized in Table 1 for the medium Earth orbit (MEO) satellites of the four current GNSSs. After introducing a tailored observation model in Sect. 2, the semi-analytical adjustment of station height, receiver clock offset, and zenith The table provides the mean orbital radius a, inclination i, orbital period T , and the ground repeat period ratio characterized by the ratio n/m of the integer number of revolutions n completed in m sidereal days. Furthermore, the number of operational satellites in May 2022 is given path delay is described in Sect. 3 along with an analysis of the elevation-dependent distribution of GNSS observations in global receiver networks. Results for the four GNSS constellations are presented and discussed in Sect. 4, considering different choices of POD/PPP processing models and standards. Finally, an independent validation of the semianalytical results using selected network solutions is given in Sect. 5, before presenting a summary and conclusions.

Observation model
The transmit antennas of the four global navigation satellite systems (GPS, GLONASS, Galileo, and BeiDou) are permanently oriented in nadir direction. From the law of sines in the Earth-station-satellite triangle, the zenith angle z of a satellite observed at elevation E = π/2 − z and the associated off-boresight (or off-nadir) angle ζ are thus related by the expression where R ≈ 6378 km is the radius of the Earth and a is the orbital radius (or semi-major axis) of the circular satellite orbit ( Fig. 1). A change of the GNSS station height by Δh results in a variation of the modeled satellite-to-station pseudorange by while a change of the phase center location in boresight direction, i.e., in the z-direction of the spacecraft and antenna reference system, by Δz PCO changes the pseudorange by Despite a somewhat different zenith angle dependence, both functions exhibit a monotonic decrease with z. In combination with the receiver clock offset, PCO changes can thus In practice, zenith-angle-dependent changes of the modeled pseudorange are also caused by variations in the tropospheric delays, which are commonly expressed as the product of a zenith delay change ΔT and an elevation-or zenithangle-dependent mapping function M(z). For the present study, we compare three mapping functions of different complexity including two fully analytical models and a tabular model representative of more elaborate tropospheric mapping functions used in state-of-the-art GNSS processing software. Neglecting Earth curvature and assuming a plane atmosphere, the mapping function can be approximated as A more realistic formulation for the wet troposphere established by Chao (1972) is given by Other than more refined mapping functions used in actual GNSS data processing, the Chao mapping function does not depend on site-and time-dependent meteorological data. It is thus convenient to use for assessing the impact of the tropospheric delay adjustment on the relation of PCO and station height changes in a GNSS PPP and POD. For further analysis, we also evaluated the mean values and the scatter of the global mapping function (GMF; Böhm et al. 2006) for the wet troposphere over a one-year period with meteorological data from the global pressure and temperature model (GPT; Böhm et al. 2007) for a network of 85 globally distributed stations over the elevations of interest. Differences between the mean GMF value and the Chao mapping function are less than 2.5% of the mapping function value for elevations down to 5 • and the difference between both mapping functions overbounds the 1-σ value of the regional and temporal GMF variations down to elevations of 1 • . Despite the lack of a closed-form analytical formulation, a tabular representation of the mean-GMF can be used as an independent model to assess the influence of tropospheric models on the ratio of PCO and station height variations and complements the two analytical mapping functions in our semi-analytical study.
Since the mean effect of station height changes, PCO changes, and tropospheric delay on the modeled pseudorange over all observed satellites will usually be absorbed in the estimated receiver clock offset dt rcv , we also consider an associated contribution in our model, where c denotes the speed of light. Overall, the pseudorange change related to the above impact factors is then described by the relation or, with Δτ = Δ(cdt rcv ) − Δh − Δz PCO , when normalizing the height and PCO contribution at zenith (z = 0) to zero.

Parameter adjustment
Equation (9) represents a reduced observation model that focuses on contributions with an explicit zenith-angle dependence of the modeled pseudorange. While ignoring numerous other contributions that are commonly considered in the modeling of GNSS observations, it is particularly suited to assess the influence of changes in station height, zenith troposphere delay, receiver clock offset, and PCO that are common to all observations of a given GNSS constellation. In particular, it allows us to assess the impact of common PCO changes on estimated station coordinates when solving for site coordinates of a global receiver network in a PPP or POD. This is accomplished by considering that range changes due to Δz PCO will be compensated by changes in Δh, ΔT , and Δτ in such a way as to minimize the square sum of the residuals Δρ, when performing a parameter adjustment in a PPP or POD. Based on Eqn. 9, the ratios between a given PCO change and the corresponding changes in height, clock, and troposphere can thus be determined by solving the relation with in a least-squares sense over the set of observations covering all stations and all satellites of the constellation. Given a discrete set of observations i = 1, . . . , n with elevation dependent weighting w i = (σ 0 /σ (z i )) 2 , the solution of the respective normal equations is given by with A i j = f i (z j ), y i = y(z i ), and W ik = w i δ ik for i, k = 1, . . . , n and j = 1, . . . , 3. In the limiting case of an infinite number of observations, the elements of the normal matrix N = (N i j ) = A T W A and the right-hand side B = (B i j ) = A T W y can be replaced by continuous integrals and covering the visible hemisphere from zenith down to the limiting elevation E min = π/2 − z max considered in the PPP or POD. The distribution function ν(z) in the integrands describes the observation density, i.e., the fractional number of observations per interval dz at zenith angle z. For a random distribution with a constant number of satellites per solid angle, it would be given by ν(z) = sin(z), while a constant function ν(z) = 2/π would apply for a uniform distribution in elevation as assumed by Cardellach et al. (2007). However, neither of these is fully realistic, and the actual distribution for all current GNSSs is in fact more closely described by a linear function ν(z) = 8z/π 2 ≈ 0.8z .
This is illustrated in Fig. 2, which shows the observation density functions for the four GNSS constellations as obtained for a representative global station network and the actual satellite orbits over one or more full ground repeat cycles (Table 1). While notable elevation dependencies may be observed for individual stations, the linear relation (16) provides a very good approximation of the globally averaged distribution and is thus well suited for the current analysis. Similar to the observation density, the elements of the normal equations depend on the weighting function w(z) = (σ (0)/σ (z)) 2 that accounts for the increase of the GNSS measurement error standard deviation σ with zenith angle. For the present study, we make use of five different weighting functions w 4 (z) = (a + (1 − a) cos(z)) 2 (a = 0.15) , and w 5 (z) = a 2 + b 2 a 2 + b 2 / cos 2 (z) (a = 5.5, b = 3.5)

Results
The relations given in the previous section provide a mathematical framework for a semi-analytical computation of PCO-related station height changes in PPP or POD solutions for a global network. Overall, the solution of the normal equations and thus the ratio of station height change and PCO change depend on the orbital radius a of the GNSS, the weighting function w(z), the elevation threshold E min , the tropospheric mapping function M(z), and the observation density ν(z). Even though some of the integrals in Eqs. (14) and (15) allow for a fully analytical treatment, we prefer numerical quadrature when forming the normal equations for selected test cases, to benefit from greater flexibility in the choice of w(z), M(z), and ν(z). Figure 3 illustrates the compensation of satellite antenna PCO changes by joint adjustment of the station height, the receiver clock offset, and the tropospheric zenith delay using the GPS constellation as an example. The test case is based on a 15 • elevation mask and the w 2 weighting function that aim to replicate the conditions of the experiment described in Zhu et al. (2003). Except for a small difference of Δα = 0.002, which possibly reflects the limited traceability of the actual processing models in that study, a good agreement of the semi-analytical model (α = −0.053) with the result of α = −0.051 given in Zhu et al. (2003) is obtained. The relative contributions of clock and troposphere changes are described by the corresponding ratios β = −0.006 and γ = +0.005 that are roughly ten times smaller than α. Corresponding results for the three other constellations are given in Case 1 of Table 2. As a general trend, the absolute value of α decreases with the altitude of the constellation and is typically about 30% smaller for Galileo than for GLONASS. This relates to the fact that the maximum offboresight angle amounts to 14.5 • for GLONASS, but only 12.4 • for Galileo. As such, PCO-induced range variations are smaller for that constellation and require a reduced change in station height for compensation. As already noted in Zhu et al. (2003), pronounced variations may also be recognized when considering different elevation thresholds (Cases 2-4). Here, lower α values are obtained, when using a lower cutoff angle E min . Even though lower elevation limits imply a larger range of off-boresight angles and therefore more pronounced pseudorange changes for a given PCO shift, their impact on station heights is reduced in this case due to an improved condition of the normal equations (Eq. (13)).
Use of different weighting functions (Cases 4-8) causes peak variations of ±0.002, which corresponds to ±2 mm station height change per 1 m of PCO shift. On the other hand, virtually no changes (Δα ≤ 0.001) were found when selecting different tropospheric mapping functions (Cases 4 and 9-10). For completeness, we note that test Case 3 largely reflects the conditions adopted in the experiment of Springer (1999). Our result of α = −0.051 shows an error of Δα = 0.002 with respect to that work (α = −0.053), which is of the same magnitude but opposite sign as obtained in comparison with the result of Zhu et al. (2003).

Validation
For validation of the semi-analytical model, dedicated PPP solutions were computed with two independent software packages, namely NAPEOS (Springer 2009) and the Bernese GNSS Software (Dach et al. 2015), using code and carrier phase observations from a global network of 60 stations in 1day batches over a 10-day period. Tropospheric delays were consistently modeled in both cases with the GMF (Böhm et al. 2006) and global pressure and temperature (GPT; Böhm et al. 2007) models, adjusting one zenith delay per 2h interval along with daily gradients. Next to the nominal PCOs as given in the IGS antenna model, solutions with few-decimeter-level offsets on the PCOs of one constellation at a time were computed. Based on these, the corresponding ratios of station height and PCO changes were obtained (Table 3). A 5 • elevation mask and the w 1 = cos 2 (z) observation weighting were applied in these PPP tests in correspondence to Case 4 of Table 2. The PPP simulation results obtained with the two software packages show a good agreement at the ±0.001level among each other and with Table 2 for all constellations, thus indicating the overall validity of the semi-analytical model.
Complementary to the PPP solutions, dedicated POD solutions with the Bernese S/W were computed to study the impact of PCO changes in the joint adjustment of GNSS satellite orbits and monitoring stations. Considering ten consecutive 1-day data sets, changes of less than 3 mm in orbital radius per meter of PCO change were observed for the individual GNSSs. This confirms the finding of Zhu et al. (2003), that the orbital radius is largely unaffected by the adopted PCO due to the dynamical constraints of the employed orbit models and provides additional justification that orbit changes can be neglected in the simplified semi-analytical model. As shown in Table 3, the ratios of station height and PCO changes obtained in POD and PPP solutions show indeed a good consistency among each other and with the analytical model. Despite inevitable limitations in the num- For comparison, the semi-analytical results from Case 4 of Table 2 are also given ber and coverage of the presented test cases, the reference PPP and POD solutions thus provide additional confidence in the validity of the semi-analytical model and support its use for parametric studies covering various impact factors such as data weighting, elevation cutoff angles, and tropospheric models.

Summary and conclusions
A semi-analytical framework is presented that allows us to assess the mutual relation of global station height changes and constellation-wide satellite phase center changes in GNSS point positioning and orbit determination without a need for comprehensive network solutions. It builds on the consideration that elevation-dependent range variations introduced by PCO changes can, to a large extent, be reproduced by a linear combination of station height changes, receiver clock offset changes, and tropospheric zenith path delay changes in the pseudorange and carrier phase observation model. For global navigation satellite systems in medium Earth orbit, representative ratios α = Δh/Δz PCO of station height and z-PCO shift between −0.04 and −0.06 are obtained. A 10 cm PCO reduction common to all satellites in a constellation will thus translate into a global increase in station heights of 4-6 mm and thus a TRF scale change of 0.6-0.9 ppb. Aside from the dominant impact of the elevation cutoff angle already identified in earlier studies for GPS, the impact of other factors such as orbit altitude, elevation dependent weighting, and the tropospheric mapping function can conveniently be studied in the semi-analytical formulation. Considering the currently operational four GNSSs, α varies by roughly 0.015 between GLONASS at about 25 000 km and Galileo at about 30 000 km altitude. It may be recognized that α is not proportional to the orbital radius, but exhibits a clearly nonlinear dependence. In the relevant altitude range of current MEO GNSSs, a 20% change in semi-major axis translates into a 25% change of α.
Variations in the elevation cutoff angle between 5 • and 15 • result in changes of α by up to 0.008, while different weighting functions in current use by IGS analysis centers result in a 0.004 spread of the associated α values. Vice versa, almost no influence is obtained for changes of the tropospheric mapping function. In comparison with network solutions for selected test cases, the semi-analytical model offers an RMS accuracy of about 0.001 in the predicted station-height-to-PCO ratios.
The results of this study are considered of particular interest in the context of GNSS PCO adjustments that are routinely required within the IGS to adjust GNSS satellite PCOs to GNSS station heights upon the release of new ITRF versions. They provide relevant background information for the understanding of differences between analysis-centerspecific PCO solutions as well as the combination of results from different GNSS constellations in a multi-GNSS processing and enable a consistent translation of existing antenna models in case of TRF scale changes.
By its very nature, the semi-analytical model presented here builds on the assumption of a globally distributed station network and uniform visibility conditions. These prerequisites are adequately met by the GLONASS, GPS, and Galileo constellation as well as the MEO component of the BeiDou-3 system, which all represent Walker constellations with multiple orbital planes and several tens of satellites. On the other hand, they do not apply equally well for the purely regional systems or regional supplements made up of geosynchronous satellite in either inclined or equatorial orbits. As a general rule, an overall decrease of α values down to roughly 0.02 is obtained for these satellites in view of their higher altitude. This means that PCO changes will only have a low impact on station height, but small station height changes will map into much larger PCO shifts than for MEO satellites. However, the exact relation cannot be predicted with the semi-analytical model and will depend more critically on the specific distribution of monitoring stations and the employed processing standards in a specific network solution.
Author Contributions OM proposed the study concept and developed the relevant algorithms building on initial ideas and contributions of PR. OM, PS, and AV performed the semi-analytical simulations and complementary PPP/POD network solutions. OM formulated the draft manuscript and prepared the individual figures. All authors contributed to the discussion and interpretation of results and critically reviewed the final manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL.
Data availability Orbital parameters of GNSS satellites as used for simulating the distribution of GNSS satellites are freely provided at https:// www.space-track.org/. GNSS observations used for validation of the semi-analytical results through PPP and POD network solutions are publicly available at https://cddis.nasa.gov/archive/gnss/data/ and other IGS data centers.

Conflict of interest
The authors have no conflicts of interest that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.