A long-term broadcast ephemeris model for extended operation of GNSS satellites

Correspondence OliverMontenbruck,DLR/GSOC, Münchener Str. 20, 82234Weßling, Germany. Email: oliver.montenbruck@dlr.de Abstract GNSS positioning relies on orbit and clock information, which is predicted on the ground and transmitted by the individual satellites as part of their broadcast navigation message. For an increased autonomy of either the space or user segment, the capability to predict a GNSS satellite orbit over extended periods of up to two weeks is studied. A tailored force model for numerical orbit propagation is proposed that offers high accuracy but can still be used in real-time environments. Using theGalileo constellationwith its high-grade hydrogenmaser clocks as an example, global average signal-in-space range errors of less than 25 m RMS and 3D position errors of less than about 50 m are demonstrated after two-week predictions in 95% of all test cases over a half-year period. The autonomous orbit predictionmodel thus enables adequate quality for a rapid first fix or contingency navigation in case of lacking ground segment updates.


INTRODUCTION
Positioning in global navigation satellite systems (GNSSs) is enabled through pseudo-distance measurements to at least four tracked satellites along with knowledge of the satellites' positions and clock offsets (Langley et al., 2017). The latter information is made available in the form of broadcast ephemerides transmitted by each individual satellite as part of its navigation message. Broadcast ephemerides describe the orbital motion and clock offset variation through a small set of parameters over a limited validity range of less than a few hours. Batches of multiple ephemeris data sets covering a larger period of time are obtained from orbit and clock determination and prediction in the ground segment and routinely uploaded to the GNSS satellites. Depending on the constellation architecture, this information is refreshed on representative time scales of one hour to a day.
Current GNSSs thus largely depend on the availability of ground-based orbit determination and time synchronization (OTDS) as well as regular upload capability to supply orbit and clock information to the users. A greater independence can be achieved through autonomous onboard orbit determination based on measurements with intersatellite links (ISLs; Ananda et al., 1990) but is not generally available in today's GNSSs. GPS has implemented ISLs in Block IIR (Rajan, 2002), Block IIF (Fisher & Ghasemi, 1999), and GPS III (Maine et al., 2003) satellites, but it is not publicly disclosed whether and to what extent they can also be used for an autonomous onboard generation of ephemeris data. BeiDou has implemented ISLs with precise ranging capability on all satellites of the BDS-3 constellation (Wang et al., 2019;Yang et al., 2019;Zhou et al., 2018) and makes use of such measurements in the ground-based ODTS process. Even though autonomous navigation concepts based on these ISLs are addressed in various studies (see, e.g., Tang et al., 2018;Guo et al., 2020), their onboard implementation status is unclear. Galileo considers ISLs to enhance the ODTS process for next generation satellites (Fernandez et al., 2010;Amarillo Fernandez, 2011), but it is likewise open whether ISL measurements will also be used for autonomous on-board orbit determination.
In the absence of full-featured autonomous orbit determination and prediction capability ("autonav"), a routine upload of broadcast ephemeris fitted to a long-term orbit prediction can be used as a backup in case of the extended loss of contact between the control segment and a GNSS satellite. A corresponding "extended mode" is defined in GPS (IS- GPS-200K, 2019), during which the satellites provide a ranging capability with reduced performance for at least 14 days after the last successful upload. The globally averaged signal-in-space range error (SISRE) degrades gracefully over time and is specified to remain below a 95th-percentile value of 388 m throughout this two-week period (DOD, 2020). This value includes orbit and clock contributions and represents an upper limit for the overall constellation performance. In practice, a notably better performance has been demonstrated in a 70-day long-term propagation test conducted with the GPS III-1 satellite in mid-2019. Here, worst-user-location SISRE values remained below 45 m after 20 days and below 600 m after 70 days (Steigenberger et al., 2020).
While the daily upload of broadcast ephemeris batches covering a long period of time can address the needs of extended operations, it implies a high uplink capacity and adequate onboard storage. Even though improved ephemeris compression techniques may contribute to a reduction of the required data volume, they are not expected to change this situation in a fundamental manner and would also require a change of the established ephemeris models and representation.
Within this study, we consider the feasibility of longterm ephemeris prediction through numerical integration. Other than analytical models, numerical orbit models can provide a close-to-reality representation of forces acting on a GNSS satellite and are therefore well suited to describe orbital perturbations over extended periods. In accord with practical needs, propagation errors in numerical orbit prediction are close to zero near the reference epoch and grow gradually over time. This is different from analytical orbit models based on "general perturbation" theories that exhibit a model-dependent best-fit error at all times due to a simplified representation of short-periodic perturbations. In principle, a single state vector and a limited set of auxiliary parameters, such as force model coefficients and Earth orientation parameters, can be used to predict a GNSS orbit over "arbitrary" times with the graceful degradation of accuracy.
Numerical orbit propagation of GPS and GLONASS was studied earlier in Seppänen et al. (2012) as a means to reduce the time-to-first-fix of consumer grade receivers in the absence of network assistance. Here, prediction periods of up to four days were considered and clock stability was identified as a key limitation for the achievable SISRE and the resulting positioning accuracy.
In the present study, we assess necessary force model improvements to bridge periods without broadcast ephemerides of up to 14 days through numerical orbit prediction. The related error budgets in terms of orbit errors, SISRE, and single point positioning (SPP) performance are discussed. As a reference case, we select the Galileo constellation, which does not presently offer ISLs or autonav capabilities, but is particularly promising in terms of long-term orbit and clock modeling. Among others, Galileo employs a highly homogeneous constellation in terms of satellite platforms, which facilitates the use of a common dynamical model in a long-term orbit propagator. In addition, Galileo offers highly stable atomic frequency standards, which helps to confine the respective SISRE contribution. On the other hand, Galileo involves a pair of satellites in moderately eccentric orbits, which provide additional insight into the capabilities of numerical as opposed to analytical long-term orbit models.
We start this study with an assessment of the Galileo orbit prediction accuracy that can be achieved with stateof-the-art models considering a ground-based GNSS orbit determination and subsequent prediction. The respective results are used to identify actual limits in the long-term predictability of a GNSS orbit and provide a reference for the performance assessment of a simplified (analytical or numerical) orbit propagator.
In a second step, a tailored force model is identified that can be used to propagate a GNSS satellite orbit based on an initial state-vector and auxiliary parameters from a ground-based orbit determination process. Aside from a trade-off of gravitational force models, a combined analytical-plus-empirical solar radiation pressure (SRP) model is presented to offer an adequate representation of non-gravitational contributions. Finally, a numerical integration scheme is identified that combines low complexity with proper efficiency and is thus well suited to be used on receiver or on-board processors with limited computational capabilities.
Conceptually, the numerical orbit propagator can be executed inside a GNSS satellite to provide real-time information on its own instantaneous position and velocity or, alternatively, inside a user receiver to compute the orbits of the entire constellation at a time. Given the computational effort associated with long-term orbit integration, GGM01S model (9 × 9), 2 tides (Rizos & Stolz, 1985), no relativity  Beutler et al., 1994) with a priori box-wing-model ; conical Earth and Moon shadow; Earth radiation pressure (Springer, 2009); antenna thrust  3-param. ECOM-1 model (estimated; Beutler et al., 1994) with a priori 2-param. cuboid model (Montenbruck et al., 2015); conical Earth and Moon shadow Numerical integration 8th-order Adams-Bashforth-Moulton multi-step predictor-corrector method with 8th-order Runge-Kutta starting step (Springer, 2009); shadow boundary handling 5th-order Dormand-Prince Runge-Kutta method (Dormand & Prince, 1980) with 4th-order interpolant (Hairer et al., 1987); optional shadow boundary handling we specifically assess the run-time characteristics on a representative micro-processor as a final step of our study.

GALILEO ORBIT PREDICTION
Propagation errors in numerical methods relate to simplifications of the underlying force models and discretization errors of the numerical integration scheme, as well as the knowledge of the actual forces acting on a satellite. Non-gravitational forces, in particular, lack an adequate predictability and pose a technical and practical limit to long-term forecasts of GNSS satellite orbits. Specific aspects include the actual orientation of the satellite body and solar panels, which affect the resulting radiation forces, thermal emissions, and even small thruster firings for attitude control or momentum dumping. In addition, the capability to properly recover the current orbit from measurements of a global network of monitoring stations directly affects the accuracy of orbit forecasts.
To assess the practical limits of orbit predictability, we determined precise Galileo orbits from observations of a global set of reference stations and subsequently predicted them over a two-week time interval. All computations were performed using the NAvigation Package for Earth Orbiting Satellites (NAPEOS v3.3.1;Springer, 2009), which implements state-of-the-art models for GNSS data processing and orbit modeling. Roughly 180 solutions with a oneday shift in the start were obtained over a half-year interval covering the October 2019 to March 2020 time frame. Relevant processing standards are summarized in Table 1.
In comparison with precise orbit products produced by analysis centers of the International GNSS Service (IGS) and its Multi-GNSS Experiment (MGEX; Montenbruck et al., 2017) and based on satellite laser ranging residuals, the resulting orbits exhibit a representative position error of less than 0.3 m (3D RMS) within the orbit determination arc. This is generally negligible compared to orbit errors in the predicted part of the resulting trajectory and allows evaluating the propagation error by comparing a predicted orbit arc with an orbit determination for the same interval. For our analysis, the RMS position errors over all active satellites have thus been determined for the 1st to 14th day of the overall prediction period for each of the daily solutions. Overall, the errors in all axes show a gradual, F I G U R E 1 Variation of Galileo orbit prediction errors in radial (left) and along-track direction (right) using high-fidelity models in the orbit determination and prediction. The graphs show the distribution (median and 5th/95th percentile) of the constellation-wide RMS error at the -th day of prediction over 180 solutions for the six-month test period [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] near-parabolic growth over time but differ substantially in magnitude. While orbit errors Δ R and Δ C in radial and cross-track direction are mostly confined to the few meters level, the errors Δ A in the along-track direction can reach a magnitude of more than 100 m after a two-week prediction.
On the other hand, user positioning accuracy is not affected equally by errors in these three directions, since the modeled range depends only on the line-of-sight component of the satellite orbit error. Considering a geographical average over the surface of the Earth as seen by a Galileo satellite , the orbit-only contribution to the instantaneous SISRE is given by and a corresponding expression applies for the RMS SISRE as a function of the RMS orbit errors. Given the relative magnitude and contribution of the individual components, cross-track errors can actually be neglected, and it is sufficient for the present discussion to consider only the radial ( = RMS(Δ R )) and the along-track ( = RMS(Δ A )) components for the SISRE contribution, of the GNSS orbit propagation errors. Statistics of the RMS orbit prediction errors obtained for the six-month test period at prediction times of 1-14 days are illustrated in Figure 1. In 50% of all cases, the root-mean-square of the along-track error over all satellites in the constellation remains confined to less than 30 m on the 14th day of the prediction arc. Considering both radial and along-track errors, this yields a median value of the geographical-average, orbit-only RMS SISRE of about 4 m and a 95th percentile value of about 11 m. All results refer to the use of the "best" models for orbit determination and prediction (see Table 1) and reflect the practical limitations in forecasting the orbital motion over extended periods of time. Given the high accuracy of models for the gravitational forces acting on a GNSS satellite, non-gravitational forces remain a major source of uncertainty in the trajectory modeling. Even though the adjustment of empirical accelerations can partly alleviate these uncertainties in the orbit determination arc, it cannot compensate for changes in the satellite behavior related to, e.g., attitude variations, solar panel orientation changes, or time-varying thermal radiation in response to the onboard temperature control. While subtle in nature, the resulting changes in acceleration can cause notable position changes when accumulating over sufficiently long periods of time.
Despite these limitations, the results demonstrate that numerical orbit propagation can, in principle, deliver ephemeris information with a degraded, but still reasonable, accuracy to support user positioning even in extended periods without broadcast ephemeris updates. When defining tailored models for autonomous orbit prediction inside a GNSS satellite or receiver, these results will serve as a guideline for the trade-off of model complexity and forecast accuracy.

TAILORED MODEL FOR NUMERICAL ORBIT PREDICTION
In a numerical orbit model, the motion of a GNSS satellite is defined by its position and velocity at an initial time 0 and a differential equation̈= ( , , ), where denotes the modeled acceleration as a function of time, position, and velocity. The resulting initial value problem can then be solved by numerical integration using a diverse set of available methods.
Within this section, we propose a tailored choice of models for reference system transformations, force models, and the numerical integration method in GNSS orbit prediction, which can be accommodated in real-time systems with limited processing capability. Compared to the highfidelity modeling employed in the previous section, the tailored models are intended to offer a notably reduced computational complexity at only a moderate loss in the prediction performance.

Reference systems and transformations
Even though orbital motion is most naturally described in an inertial or celestial reference frame (CRF), broadcast ephemerides of GNSS satellites are traditionally expressed in a terrestrial reference frame (TRF). This choice is particularly convenient for the user segment, since it avoids the need for explicit evaluation of the CRF-to-TRF transformation when computing user positions in an Earth-fixed system.
Rather than performing the orbit prediction in a CRF and transforming the result into the TRF, a formulation of the equation of motion in the rotating TRF may therefore be considered. This approach has long been adopted for the broadcast ephemeris models of GLONASS and SBAS and was likewise found advantageous for real-time onboard navigation systems of low Earth satellites (Montenbruck & Ramos-Bosch, 2008).
For an arbitrary CRF-to-TRF transformation matrix , the equation of motion in the TRF involves additional Coriolis and centrifugal accelerations, which are given by Considering the near constant Earth rotation rate and the close alignment of the instantaneous rotation axis with the -axis of the TRF, Equation (3) can notably be simplified. This yields the familiar expression: (see, e.g., GLO-ICD-5.1, 2008), which adds no relevant complexity to the equation of motion as compared to an inertial formulation and is routinely used in GLONASS/SBAS receivers for evaluation of the respective broadcast ephemerides. However, the same approximation can no longer be applied when integrating the equation of motion over time scales of days or weeks. While apparently small, the neglected contributions to the Coriolis and centrifugal acceleration would result in a notable mismatch of the propagated orbit. A rigorous evaluation of Equation (3), in contrast, implies computation of the time derivativesȧ nd̈. In view of numerical errors, approximation errors, and/or an undue complexity of the resulting expressions, neither difference quotient approximations nor analytical differentiation were found to provide derivatives of adequate accuracy and computational effort. Overall, the effort for a proper evaluation of the Coriolis and centrifugal accelerations rapidly compensates the potential savings of the TRF model in orbit predictions over extended time scales. A CRF formulation of the equation of motion is therefore preferred for use in the long-term orbit propagation.
Common models (Petit & Luzum, 2010;Seidelmann, 2006) describe the transformation as a sequence of elementary rotations with rotation angles defined through polynomials, harmonic series, and observed tabulated quantities. Uncertainties in individual angles do not affect the radial position of the predicted satellite position and also leave the relative geometry of the entire constellation unaffected. However, they result in rotational misalignment and ultimately a shift of the user position on the Earth's surface. By way of example, a 1 ′′ angular error in the CRF-to-TRF transformation translates into a 150 m cross-track or along-track error of a GNSS satellite at 30,000 km altitude and a 30 m North-South or East-West error of the estimated user location. Given the overall error budgets for long-term broadcast ephemerides in extended-mode GNSS operations, a 0.01 ′′ truncation level is considered adequate for individual contributions to the overall transformation.
Since the full accuracy of IERS 2010 transformations is not required for the present purpose, we can safely adopt the legacy IAU 1976/1980 models (Seidelmann, 2006) of precession, nutation, and sidereal time, along with observed values of the Earth orientation parameters. The overall transformation is given by where , , denotes elementary rotation matrices for the respective axes. The precession angles ( , , and ) as well as the obliquity of the ecliptic ( ) are described through third-order polynomials in time, which place no relevant computational burden. The nutation angles Δ and Δ , in contrast, are expressed through harmonic series requiring extensive trigonometric function evaluations.
Compared to the current IAU 2000A nutation model with a total of 1,365 terms (Petit & Luzum, 2010), the IAU 1980 nutation model proposed here for the long-term orbit propagator comprises only 106 terms. Still, evaluation of all terms may induce a heavy computational load for real-time systems. Given the previously defined accuracy threshold, a reduced set of 18 terms with amplitudes above five milli-arcseconds is, in principle, sufficient when converting propagated CRF positions into the TRF. On the other hand, a roughly ten times better accuracy of the nutation model may occasionally be required in the orbit determination, which is used to adjust initial conditions for the longterm propagation. During the eclipse season, in particular, even small errors in the observation model may affect the estimation of empirical SRP parameters with an adverse impact on the drift and along-track motion during the subsequent long-term propagation.
It is therefore preferred to make use of an extended subset of the IAU9180 nutation series with 49 terms >0.5 mas. For a computationally lean implementation, substantial savings may, e.g., be achieved through Chebyshev approximation of the slowly varying nutation angles Δ and Δ . Alternatively, the explicit evaluation of trigonometric functions in the nutation series may be minimized through systematic use of the addition theorems. These allow us to compute the sine and cosine of linear combinations of the mean arguments of luni-solar motion from the sine and cosine of the elementary angles using basic arithmetic operations (Montenbruck & Pfleger, 2000). At the expense of a moderately increased code complexity, substantial reductions in the overall computing time may thus be achieved.
Greenwich mean sidereal time (GMST) in Equation 5 is related to the UT1 time scale through a third-order polynomial, which poses no computational burden, but requires knowledge of the time varying ΔUT1 = UT1 − UTC difference. Along with the pole coordinates p and p , ΔUT1 constitutes a set of three Earth orientation parameters (EOPs) that cannot be described through models but must be retrieved from actual observations of Earth rotation and extrapolated for use in the prediction of GNSS ephemerides.
Since Earth rotation exhibits notable periodic variations due to luni-solar tides, it is common practice to extrapolate a regularized version (UT1R) in which zonal tides with periods of 5 to 35 d have been removed and to add these terms back to the extrapolated UT1R to obtain the prediction of UT1 itself. Conventional expressions of UT1 − UT1R involving 41 harmonic terms are provided in Petit and Luzum (2010) and used within this study. Even though the overall amplitude of these terms is less than 2.5 ms, or, equivalently, 5 m at GNSS altitude, their neglect may adversely affect the adjustment of drift and SRP parameters with notable impact on the long-term prediction. Furthermore, diurnal and semi-diurnal tides need to be added to tabular as well as extrapolated p , p , and ΔUT1 values. To limit the computational effort, the eight-term series of McCarthy (1996) is employed.
For use in specialized, e.g., spaceborne, navigation applications, most GNSSs have incorporated EOP information in their modernized civil navigation messages. EOP predictions including epoch values of p , p , and ΔUT1 along with first-order time derivatives are presently transmitted by the GPS IIR-M/IIF/III and GLONASS-K1 satellites, as well as the BeiDou-3, QZSS, and IRNSS constellations as part of their modernized civil navigation messages. A systematic performance analysis of broadcast EOPs has not been performed so far, and it remains unclear to what extent these data would be suitable for predictions over the two-week time frame considered here. For Galileo, EOPs are not part of the existing INAV and FNAV navigation messages and would in any case need to be included into the ephemeris parameters of a possible long-term broadcast ephemeris model.

Gravitational forces
GNSS satellites in medium Earth orbit operate at altitudes of roughly 4-5 Earth radii, which results in a rapid attenuation of high-degree and -order contributions of the Earth-gravity field. A 12 × 12 expansion is commonly used in precise GNSS orbit determination, but an even smaller 9 × 9 is considered adequate for use in a long-term orbit propagator. Compared to the slightly larger standard, it offers a small advantage in terms of computational load with no relevant impact on accuracy (Table 2). Computationally efficient and stable recurrence relations for computing the acceleration from a given set of Earth gravity model parameters are given in Cunningham (1970). For the present study, the GGM01S gravity model obtained from early measurements of the GRACE mission (Tapley et al., 2004) is adopted. It is intentionally different and older than the models used in Section 2 to show the potential impacts of a reduced quality on the prediction performance. GGM01S involves linear variations of the second-order terms consistent with the definition of the Earth rotation pole but is otherwise static. Given the targeted propagation accuracy, no need is identified to account for seasonal gravity field variations in the present application. However, gravity field variations caused by deformations of the Earth due to luni-solar tidal forces need to be taken into account in at least the lowest TA B L E 2 Performance comparison of model options. Peak errors in 14-day prediction arcs following a two-day orbit determination have been evaluated for eight test dates spread over a three-month period. Median and maximum errors over the tests are provided for various forms of simplifications relative to a comprehensive reference along with the approximate reduction in computational load. Models marked by a × symbol are recommended for use in the long-term propagator degree/order. Compared to a neglect of solid Earth tides, which may cause errors at the 50 m level after the twoweek propagation, the consideration of the 2 contribution (Rizos & Stolz, 1985) already provides a more than ten-fold accuracy improvement (Table 2). Higher degree and order contributions to the solid Earth tides, as well as pole and ocean tides, need not be considered in the present context. Relativistic effects in the motion of Earth orbit satellites are commonly described through post-Newtonian corrections of the equation of motion. The dominating Schwarzschild correction (Petit & Luzum, 2010) results in a perigee rotation of about 1 mas/d (Hugentobler & Montenbruck, 2017) for GNSS satellites. When neglected in the combined orbit determination and prediction, errors at the few-meter level may arise (Table 2), which is well tolerable for the proposed long-term orbit propagator.
In terms of third-body perturbations, a point-mass model for the gravitational acceleration of the Sun and Moon is sufficient, while planetary perturbations can safely be neglected. Since precomputed ephemerides requiring a large amount of tabular data are hardly suitable for use in embedded systems, analytical models of luni-solar motion are preferred for the definition of the longterm orbit propagator.
As shown in Hugentobler and Montenbruck (2017), representative solar perturbations of a GNSS satellite amount to roughly 1 to 2 km after two revolutions for a given initial state, and a 2-3 times higher value applies for perturbations by the Moon. A relative accuracy of typically 4⋅10 −4 , or, equivalently, 50 ′′ in the computation of the lunar and solar position might thus be considered adequate to obtain a propagation error of a few tens of meters over a two-week interval. However, similar to the nutation model discussed above, much larger errors may again arise in a combined orbit determination/prediction due to unfavorable coupling of unmodeled perturbations with the estimation of the semi-major axis and empirical SRP coefficients.
In case of solar perturbations, a Keplerian approximation of the Earth's motion around the Sun is found to be only marginally sufficient for the purpose. Consideration of the dominating planetary perturbations and the offset of the Earth's center from the Earth-Moon barycenter is therefore recommended. Compared to a comprehensive analytical solar ephemeris (Montenbruck & Pfleger, 2000),  (Table 2). In view of larger perturbations, an even better accuracy is required in the formulation of the lunar ephemeris. Here, an analytical model based on the Improved Lunar Ephemeris (USNO, 1954) with mean elements of the Ephemeride Lunaire Parisienne (ELP2000; Chapront-Touzé & Chapront, 1988) has been adopted, and a 2 ′′ truncation is applied in the harmonic series. A computationally efficient implementation that minimizes the evaluation of trigonometric functions in this model is described in Montenbruck and Pfleger (2000).

Non-gravitational forces
Non-gravitational forces acting on a GNSS satellite are dominated by solar radiation pressure. The corresponding acceleration reaches a magnitude at the 100 nm/s level and causes km-level perturbations after a few days of prediction (Hugentobler & Montenbruck, 2017). Other contributions such as Earth albedo, thermal radiation, and antenna thrust are roughly two orders of magnitude smaller and are neglected in the present context in view of insufficient or complex a priori models. Solar radiation pressure acting on the Galileo satellites can largely be described by a three-parameter cuboid boxwing model for a stretched satellite body with solar panels in a nominal yaw-steering attitude (Montenbruck et al., 2015). The model parameters comprise two characteristic accelerations ( C , S ) related to absorption and diffuse reflection of cubic and stretched body contributions, as well as the total solar panel acceleration, sp . Remaining imperfections can be covered by the empirical CODE orbit model (ECOM; Beutler et al., 1994) considering a limited set of three empirical accelerations ( 0 , 0 , * ) that are adjusted as part of the orbit determination process.
Both, the a priori and empirical SRP contributions are described in a DYB frame aligned with the Sun direction (D), the direction perpendicular to the Sun-spacecraft-Earth plane (Y), and the B-axis orthogonal to D and Y. The respective components of the radiation pressure acceleration are given by where denotes the Sun-spacecraft-Earth angle. As a modification of the original ECOM formulation introduced in Beutler et al. (1994), the argument of latitude in the harmonic terms is replaced by the orbit angle since local midnight. Parameterization of the orbit periodic contributions in terms of naturally accounts for the inherent symmetries in the SRP acceleration due to a regularly shaped satellite body (Montenbruck et al., 2015;Arnold et al., 2015) and allows dropping of the sinusoidal term required in the original ECOM formulation. Use of the a priori box-wing model likewise reduces the number of empirical terms that would otherwise be required for a proper SRP modeling and allows us to confine the empirical accelerations to small values. Given the very similar properties of the In-Orbit Validation (IOV) and Full Operational Capability (FOC) satellites , a common parameter set can be used for the entire constellation. Values adopted for the long-term propagator are summarized in Table 3.
The modeling of SRP forces is further complicated by the fact that GNSS satellites are not permanently sunlit. For Galileo, regular eclipse periods with Earth shadow transits arise typically twice per year whenever the Sun elevation above the orbital plane is less than roughly 15 • . In addition, Moon shadow transits may occur on an irregular basis. While comparatively rare, the neglect of such eclipses would result in undue propagation error. By way of example, the Moon shadow transit of Galileo satellite E26 on November 26th, 2019, causes a drift change of roughly 80 m/d. Given their distance from the center of the Earth or Moon, GNSS satellites experience a nonnegligible penumbra phase at the begin and end of an eclipse during which the radiation pressure varies gradually between zero and its full value. A conical shadow model is therefore required for the long-term propagator. This is particularly true for lunar shadow transits, which include at best a very short umbral phase. Within an eclipse, the modeled SRP acceleration as given in Equation 6 is multiplied by an illumination factor 0 ≤ < 1 that accounts for the reduced amount of solar energy received by the satellite.

Integration and interpolation
The choice of a suitable numerical integration method for the long-term orbit propagator is driven by various considerations. Single-step methods are generally preferred for an embedded system to simplify bookkeeping as well as a start and restart of integration at discontinuities. When computing navigation solutions in a receiver at representative update rates of 1 Hz or up, GNSS satellite positions must be available at an equally high sampling rate. In order not to constrain the step size in an undue manner, numerical integration with dense output requires the availability of an interpolant of corresponding order, ideally with no need for additional evaluations of the equation of motion. High-order methods offer improved overall efficiency in terms of the accuracy vs. number of function evaluations over large time intervals, but do not commonly support dense output and become less efficient when discontinuities enforce an integrator restart. Finally, a particularly high efficiency is offered by methods for the direct integration of the secondorder equation of motion̈= ( , ) that does not depend on velocity. Typically, this condition is met for all satellites outside the atmosphere when using an inertial formulation of the equation of motion, but is still violated in the present case. Orbit-periodic accelerations in the ECOM1 radiation pressure model (Section 3.2) are formulated in terms of the orbit angle , which is measured along the orbital plane. Computation of thus necessitates knowledge of the velocity in the acceleration model and therefore inhibits use of, e.g., Runge-Kutta-Nystrøm methods. Two low-order single-step methods for dense output that have earlier been studied for use in real-time navigation systems (Montenbruck & Gill, 2001) include a combination of the classical fourth-order Runge-Kutta with Richardson extrapolation (RK4R; Hairer et al., 1987) and the fifth-order RK5(4)7M Runge-Kutta method of Dormand and Prince (1980). Both methods are effectively of the fifth order and require 5.5-6 stages, i.e., evaluations of the differential equation, per integration step. The Dormand-Prince method, subsequently designated as DP5, offers an embedded fourth-order solution for step size control at the expense of a seventh function evaluation. However, given the near circular nature of GNSS satellite orbits, no step-size control is required here, and only the fifth-order/six-stage formula is therefore considered. Representative step sizes of the order of 100 to 300 s can be used for GNSS orbit propagation with the two methods. State vectors at intermediate epochs are conveniently obtained from a fifth-order Hermite interpolation in case of the RK4R method or a fourth-order interpolant in case of DP5 (Hairer et al., 1987).
All numerical integration methods are based on the inherent assumption that the equation of motion is continuously differentiable up to the order of the method. For satellite orbits, this condition appears readily fulfilled in maneuver-free periods, but is actually violated when considering shadows of, e.g., the Earth or Moon in the SRP modeling. This is obvious for a cylindrical shadow model that causes a discontinuous acceleration, but likewise applies for a conical shadow model, which exhibits a continuous but non-differentiable acceleration. As pointed out in Lundberg et al. (1991), Lundberg (2000), and Woodburn (2000), neglect of shadow boundaries may cause pronounced orbit propagation errors unless appropriately handled by a restart of the integrator at each discontinuity. Shadow boundaries are particularly critical when using high-order integration methods but have a less pronounced impact when working with low-order methods and short step sizes.
To assess the performance of the two methods at different step sizes, orbits of the Galileo constellation have been propagated over a two-week period starting on March 10, 2020 using the tailored force model of Table 1. The test epoch coincides with a phase of deep eclipses for all satellites in plane C and partial eclipses of the Galileo E14 and E18 satellites in eccentric orbits for studying the impact of shadow transits on the overall accuracy. For assessment of propagation errors with different integration options, a solution based on short (50 s) steps and including shadow boundary handling is considered as a reference. Propagated positions were evaluated at a 270 s sampling that is not an integer multiple of the step size. This ensures that a substantial fraction of the output points considered in the performance analysis is obtained from interpolation rather than just the beginning and end of a step. All computations were performed using eight-byte (double-precision) floating point data types, which provide a 10 −16 precision of individual arithmetic operations. With rounding errors of this magnitude, the global integration errors over the 14-day prediction are dominated by the approximation of the employed fifth-order methods rather than numerical precision over the relevant range of step sizes.
The results compiled in Table 4 show a clearly better performance of the DP5 method over RK4R at the identical step size, which cannot be compensated by the slightly smaller number of stages. Even though both methods share the same order, DP5 exhibits a much lower integration error. Given a similar or even reduced complexity of implementation, DP5 is clearly preferable for the present application, even though the associated interpolant exhibits a slightly lower order than the integration method itself.
In all cases, the total position error due to numerical integration is dominated by the along-track component, and peak values are achieved near the end of the prediction arc. Comparison of individual satellites shows that the peak errors are driven by the two Galileo satellites (E14/E18) in slightly eccentric orbits. By way of example, this is shown in Table 4 for the DP5 integrator at 200 s step size but applies equally for the remaining test cases. Satellites in circular orbits show a quadratic growth of propagation errors over time that is almost identical for all non-eclipsing satellites. With shadow boundary handling, the propagation error of eclipsing satellites in circular orbits matches that of the non-eclipsing satellites. Without shadow boundary handling, a quadratic growth can also be observed for the eclipsing satellites, but the overall amplitude differs among the individual satellites ( Figure 2). Despite small advantages, no pronounced benefit of the shadow boundary handling can thus be identified for the given order and step size of the integrator. For completeness, we note that a notably different growth of integration errors is obtained when determining the initial conditions of the orbit propagation in a preceding orbit determination that uses the same dynamical model and integrator as the prediction part (Table 5). Here, errors are mostly lower than in the prediction with fixed initial conditions, and a more pronounced benefit of the shadow boundary handling is obtained. In fact, the orbit adjustment allows us to compensate for numerical integration errors in the subsequent prediction, which reflects itself in solution-specific values of the estimated semi-major axis and empirical non-gravitational force model parameters.
Based on the test results, use of the DP5 integrator and its associated interpolant at a step size of ℎ = 200 s is identified as a suitable compromise between accuracy and computational effort for the long-term orbit propagator. With peak values at the 10 m level, the numerical integration errors remain well below the inherent predictability discussed in the previous section even for the more demanding case of the two Galileo satellites in eccentric orbits. At the given step size, integration across shadow boundaries is tolerable for the orbit prediction, which greatly simplifies the user implementation. On the other hand, a rigorous shadow boundary handling is clearly advisable for the orbit determination and parameters adjustment.

CLOCK PREDICTION
Next to orbit errors, the user positioning accuracy is directly affected by uncertainties in the predicted satellite clock offsets. Other than orbit prediction, which is largely driven by the quality of the dynamical models, the performance of satellite clock prediction is governed by the inherent clock stability. Galileo is renown for the first use of passive hydrogen masers as a time and frequency standard in GNSS satellites (Falcone et al., 2017), and numerous reports have highlighted the exceptional quality and in-orbit performance of these clocks Galluzzo et al., 2018;Huang et al., 2019). During the September 2019 to April 2020 test period TA B L E 5 Peak orbit propagation errors (3D, in m) over two weeks after a two-day orbit determination for different integration methods, step sizes, and shadow boundary handling considered in the present study, all satellites of the Galileo constellation were operating on H-masers, except for the first IOV satellite GSAT101 (E11) that used one of its Rubidium atomic frequency standard (RAFS). Representative Allan deviations of the H-masers amount to 10 −12 ∕ √ ∕s at time scales of 1 to 100, 000 s (Beard & Senior, 2017) corresponding to stochastic clock variations at the sub-nanosecond level around a polynomial clock model over time scales of up to a day. For 24 h clock prediction, peak errors of about 3 ns (1 m) for H-masers and 8 ns (2.5 m) for RAFS are reported in Falcone et al. (2017) for the first Galileo satellites, but only limited information is available for extended forecasts.
For the assessment of clock prediction performance over extended time intervals, Galileo clock offsets obtained in a precise orbit and clock determination process (see Table 1) were approximated by a low-order polynomial: over a two-day interval. Predicted values for up to 14 d beyond the end of the data interval were then compared against precise clock offsets. For all H-masers, best results were obtained with a first-order clock polynomial (i.e., Major changes in clock drift could be observed on satellites E11 (Δ 1 ≈ 12 μs∕d) and E25 (Δ 1 ≈ 1 μs∕d) in November 2019 following a temporary off-period. In the same month, gradual frequency adjustments by 50 to 120 ns/d within 5 to 15 days took place on satellites E07, E30, and E33. Otherwise, all H-masers showed extremely stable clock drifts with typical variations of less than ±5 ns/d over many months (Figure 4).
Other than orbit prediction, which shows a largely homogeneous performance and gradual growth across all satellites, the prediction of clock offset variations in the Galileo constellation reveals a less uniform picture. While most satellites exhibit a surprisingly good performance and allow for a 5 to 10 m RMS error for most of the time even after a two-week prediction, a small subset of satellites may be affected by unforeseeable clock drift changes that can F I G U R E 4 Clock drift variation of Galileo satellites over the sixmonth test period [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] cause errors of 100 to 1,000 m. While exclusion of satellites with Rb-clocks (currently only E11) can certainly alleviate this problem, undesirable frequency changes on the H-masers are also encountered on some of the Galileo satellites in the first months of our test period. On the other hand, the limited number of satellites in the constellation that is concurrently affected by large clock prediction errors and the magnitude of these errors would generally allow us to discard these satellites as part of the receiver autonomous integrity monitoring (RAIM; Parkinson & Axelrad, 1988). Based on this consideration, we exclude satellites with errors larger than 100 m in the predicted clock offset from the statistics but otherwise incorporate all satellites of the constellation irrespective of the employed clock type.
Results in Figure 5 show that the median RMS clock prediction error with 100 m outlier screening increased gradually to about 8 m over the first to seventh day of prediction and essentially stays at this magnitude up to day 14. 95th percentile errors amount to roughly 20 m in our test period. Compared to orbit prediction errors discussed in the previous sections, a balanced contribution of orbit and clock errors to the overall SISRE can thus be achieved for long-term broadcast ephemeris prediction in the Galileo constellation due the use of high performance atomic frequency standards.

PERFORMANCE CHARACTERIZATION
Overall, the long-term orbit propagator concept comprises the provision of an epoch state vector, empirical solar radiation pressure coefficients, and a clock offset polynomial for each satellite, as well as a set of Earth F I G U R E 5 Variation of Galileo clock offset prediction errors. The graph shows the distribution (median and 5th/95th percentile) of the constellation-wide RMS error at the -th day of prediction over 180 solutions for the six-month test period excluding individual clock offset errors larger than 100 m [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] orientation parameters valid for the entire constellation. For a comprehensive assessment, the corresponding parameters were generated over a six-month time frame from October 2019 to March 2020. After characterizing the computational effort on a representative embedded computer system, the achieved accuracy is assessed by comparing the propagated orbit and clock states against a reference solution and by evaluating single point positioning errors.

Parameter adjustment
As discussed in Section 3.1, Earth orientation parameters are required in the long-term orbit propagator concept to transform the propagated state into the terrestrial reference for user positioning. Since availability of EOPs and associated rates enabling a two-week propagation is not ensured in current GNSSs, linear polynomials for pole and UT1R were adjusted to the observed EOPs from the igs96p02 rapid EOP series over a two-day arc prior to the arc used for the orbit estimation. Daily orbit and clock products based on 24 h of observations were used to adjust the orbit and clock parameters for the individual Galileo satellites. For improved accuracy of drift terms and empirical accelerations, precise ephemerides from two consecutive days were concatenated prior to fitting the individual parameters of the long-term propagator. Longer orbit determination arcs did not show relevant changes in the overall statistics. TA B L E 6 Long-term orbit propagator parameters and recommended numerical representation. The number of bits in column 4 refers to a single parameter of the specified type. Depending on the parameter generation process, common epochs may be adopted for EOP and orbit/clock parameters to reduce the total amount of data transmitted to the user Clock offset polynomials of the first or second order were obtained for satellites operating hydrogen masers and the rubidium clock, respectively. To avoid a need for external information, the clock type and required polynomial order were automatically determined based on the clock residuals of a linear fit using a 1 ns threshold. Finally, the epoch state vector as well as the coefficients of the three-parameter ECOM-1 empirical solar radiation pressure model were adjusted using positions from the precise ephemerides as pseudo-observations. Force model and numerical integration options for the orbit adjustment were chosen to match the recommended long-term orbit propagator settings (Table 2) to ensure full consistency of the models on the provider and user side. Table 6 summarizes the set of long-term orbit propagator parameters and provides a recommendation for the required numerical range and precision. The number of bits and the least-significant bit (LSB) suggested for the EOPs match the values currently adopted in the CNAV navigation message of GPS, except for the EOP epoch that is parameterized with an hourly resolution. LSBs for the epoch state vector are based on the consideration that a semi-major axis error Δ gives rise to an along-track error of roughly 10Δ after one revolution and about 250Δ after two weeks. Centimeter-level uncertainties in the position, or, equivalently, micrometer-level velocity errors will thus give rise to meter-level along-track errors at the end of the envisaged propagation time span. LSBs for the clock polynomial coefficients are chosen such as to confine the discretization error to roughly 1 ns over a 14 d interval, and the total number of bits is selected such as to keep the supported range of values compatible with the present Galileo navigation message.
For the performance assessment, all estimation parameters were rounded to the bit resolution specified in Table 6 after the EOP, epoch state, and clock adjustment. Discretization errors are thus considered in the overall error budget.

Orbit and range error
The overall performance of the long-term orbit propagator is summarized in Figure 6. Radial errors, which map directly into the modeled range, reach a median RMS of 4 m after two weeks and a peak RMS value of 7 m over the half-year test period. In the along-track direction, the daily RMS error of the predicted positions exhibits a 40 m median, while the peak RMS error over all days amounts to 90 m. Roughly one eighth of this error, i.e., 5 m (median) and 11 m (peak) would map into the line-of-sight direction on a global average, but worst-user location errors can amount to roughly twice these values. Compared to the results of Figure 1 obtained with "best possible" dynamical models, the tailored models yield a degradation of roughly a factor of two in the radial component but less than 30% in the along-track direction. This illustrates again that the achievable long-term orbit prediction accuracy is driven by unpredictable motion changes that cannot adequately be described by a priori or empirical models.
The combined effect of orbit and clock prediction errors is best described by the associated signal-in-space range errors. Statistics for the globally averaged RMS SISRE over the half-year test period are shown in Figure 7. Similar to the clock error analysis in Figure 5, a 100 m SISRE threshold was employed to eliminate coarse outliers. The median RMS SISRE increases almost monotonically over time and attains a value of roughly 12 m after a 14-day prediction. For 95% of all days in the test period, the SISRE after two weeks remains below 22 m.

Positioning error
As a final test, the single-point positioning accuracy using the predicted ephemeris data has been evaluated for a representative mid-latitude user location (11.6 • W, 48.1 • N).

F I G U R E 6
Variation of Galileo orbit prediction errors in radial (left) and along-track direction (right) using the tailored models of Tables 1 and 2 in the orbit determination and prediction. The graph shows the distribution (median and 5th/95th percentile) of the constellation-wide RMS error at the -th day of prediction over 180 solutions for the six-month test period [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] F I G U R E 7 Evolution of Galileo global average RMS SISRE using the tailored models of Tables 1 and 2 in the orbit determination and prediction [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com and www.ion.org] Results in Figure 8 include the combined effects of line-of-sight orbit errors and clock errors for an individual station taking into account all visible satellites above a 5 • elevation limit. Furthermore, a position dilution of precision (PDOP) threshold of ten and a 100 m outlier rejection limit were applied. Median values of the 3D position error grow roughly linearly with the forecast period and reach a value of about 23 m on day 14. For 95% of all days in the test period, the 3D RMS position errors remain below 50 m. Overall, the position error statistics correspond closely to the SISRE statistics discussed before for a representative PDOP of 2-2.5.

Run-time performance
Aside from accuracy, the feasibility of a long-term orbit propagator also depends on the computational effort. The tailored orbit models presented in Section 3 offer a globally averaged orbit-only SISRE of less than 15 m RMS in 95% of all cases after a two-week prediction period. This is roughly compatible with the achievable clock prediction accuracy for Galileo Hydrogen masers over the same interval and yields a balanced contribution of both error sources to the overall SISRE. On the other hand, the model exhibits a notable algorithmic complexity, and the choice of a numerical orbit propagation concept for the long-term ephemeris model implies a linear growth of the computational effort with the time between the reference epoch of the state vector and the epoch at which a GNSS satellite position is evaluated. This is less critical during continued integration in real time, but may raise concerns after extended off-periods of a receiver.
For an assessment of the run-time performance, the execution time for propagation of a single GNSS satellite orbit over two weeks at a 200 s step size and using the recommended model options of Table 2 was determined for two different 32-bit processors that are widely used in embedded systems (Table 7). They comprise an ARM Cortex A9 processor as used in the Altera Cyclone V System on Chip (SoC) field programmable gate array (FPGA) and a single-core ARM1176JZF-S processor integrated in the Broadcom BCM2835 SoC of a Raspberry Pi 1B+ miniature computer. These processors exhibit a computing power that is roughly two orders of magnitude lower than that of common desktop computers.
In a worst-case scenario, the orbit propagation for the entire GNSS constellation with 30 satellites would thus require 5 to 15 min on these processors when operating a receiver with a two-week old, long-term ephemeris data set. Even though this is comparable to the transmission time of a GPS almanac and probably still tolerable, a more realistic operation scenario would assume operation of a receiver at least once per day. In this case, the propagated state vectors could be stored in nonvolatile memory prior to shut down, and the time to reactivate the receiver a day later would amount to only 20 to 65 s. Both times appear compatible with actual user requirements for operation of a long-term ephemeris model and justify consideration of the approach for reduced-accuracy navigation in the absence of regular broadcast ephemeris updates.
As an alternative use case, one may consider an onboard implementation of the orbit model to continuously propagate its own orbit during extended periods without ground contact. The resulting state vector could then be downloaded to the users in regular intervals using a GLONASS or SBAS-like ephemeris message, or an adjusted set of quasi-Keplerian orbit parameters as used in the common GPS, Galileo, and BeiDou ephemeris models. While still depending on the ground-based genera-tion of the long-term ephemeris data, this architecture would enable a fully autonomous navigation payload operation. Obviously, this concept notably reduces the complexity of user segment modifications, but would require adequate spare resources in the onboard processor of each GNSS satellite.

SUMMARY AND CONCLUSIONS
The feasibility of predicting GNSS satellite orbits and clocks over periods of up to two weeks has been studied using the Galileo system as a test case. A tailored set of force models and reference system transformations is identified, based on which GNSS satellite motion can be propagated using numerical integration. The models are likewise applicable for near-circular and slightly eccentric orbits and properly account for eclipses. Along-track errors of about 40 m median RMS are only moderately worse than those of high-fidelity models and essentially represent the inherent predictability of GNSS orbits in the presence of characterized, satellite-related perturbations.
Given the high performance of the Galileo Hydrogen masers, which allow for clock prediction accuracies of better than 1 m/d, the combined orbit and clock prediction over a two-week time frame enables a 95th-percentile global-average RMS SISRE of less than 25 m. This translates into a 3D RMS user position accuracy of better than 50 m (95th percentile) or 23 m (median). For comparison, a 95th-percentile SISRE of better than 388 m is specified for GPS extended mode operations over a two-week period.
Despite its obvious complexity, the long-term orbit propagator is shown to be compatible with representative processing capabilities of embedded processors and may be considered for either user segment or space segment implementation. Compared to the present GPS extended mode architecture, use of the propagator notably reduces the data volume for the uploading of short-arc broadcast ephemerides. It therefore lends itself as an interesting alternative for the implementation of extended operation modes in GNSSs that do not presently offer autonomous navigation capabilities.

A C K N O W L E D G M E N T S
The International GNSS Service (IGS) is acknowledged for providing GNSS observation data. We would also like to thank Dr. Paolo Zoccarato of the Radio Navigation Systems and Techniques (TEC-ESN) section at ESTEC/ESA for the constructive comments received during manuscript preparation.
Open access funding enabled and organized by Projekt DEAL.