Combination strategy for consistent final, rapid and predicted Earth rotation parameters

The accurate knowledge of the Earth’s orientation and rotation in space is essential for a broad variety of scientific and societal applications. Among others, these include global positioning, near-Earth and deep-space navigation, the realisation of precise reference and time systems as well as studies of geodynamics and global change phenomena. In this paper, we present a refined strategy for processing and combining Very Long Baseline Interferometry (VLBI), Satellite Laser Ranging (SLR), Global Navigation Satellite Systems (GNSS), and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) observations at the normal equation level and formulate recommendations for a consistent processing of the space-geodetic input data. Based on the developed strategy, we determine final and rapid Earth rotation parameter (ERP) solutions with low latency that also serve as the basis for a subsequent prediction of ERPs involving effective angular momentum data. Realising final ERPs on an accuracy level comparable to the final ERP benchmark solutions IERS 14C04 and JPL COMB2018, our strategy allows to enhance the consistency between final, rapid and predicted ERPs in terms of RMS differences by up to 50% compared to existing solutions. The findings of the study thus support the ambitious goals of the Global Geodetic Observing System (GGOS) in providing highly accurate and consistent time series of geodetic parameters for science and applications.


Introduction
The Earth orientation parameters (EOPs) describe variations in the rotation of the Earth in space. The orientation of the Earth's rotation axis, seen from space, is given by the celestial pole coordinates X and Y . Seen from the Earth's surface, it is given by polar motion (PM), i.e. the terrestrial pole coordinates x pole and y pole . In addition, the EOPs describe the spin of the Earth about its rotation axis (rotation angle, UT1). Besides these offset parameters, the EOPs also comprise the temporal variation in the Earth's orientation and rotation in the form of celestial pole and PM rates and the length of day (LOD). Altogether, the EOPs describe the transformation between the terrestrial reference system (TRS) and the celestial reference system (CRS; Petit and Luzum 2010). The sub-group comprising the PM components, their rates, UT1 and LOD, is called the Earth Rotation Parameters (ERPs), which are in the focus of the study presented hereafter. They undergo irregular variations as they are largely influenced by geodynamic processes (Lambeck 1980;Gross et al. 2003Gross et al. , 2004Gross et al. , 2005Seitz and Schmidt 2005). Retrospectively, they are accessible via observations from the four space-geodetic techniques, i.e. Very Long Baseline Interferometry (VLBI), Satellite and Lunar Laser Ranging (SLR, LLR), Global Navigation Satellite Systems (GNSS), and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), whereby it is important to note that VLBI is the only technique giving access to all ERPs including UT1, i.e. the Earth's absolute rotation angle with respect to the inertial space, whereas the satellite techniques are not sensitive to this parameter. Prospectively, the ERPs can be predicted over a limited time span based on geophysical fluid models, time series analysis or machine learning techniques.
Changes in the Earth's rotation are caused by external gravitational forces and geodynamic processes that exchange effective angular momentum (EAM) between the solid Earth, its interior and its fluid envelope (including the atmosphere, the oceans and the continental hydrology). Therefore, accurate observation, modelling and understanding of the ERPs contribute substantially to studying Earth system processes across a wide range of time scales.
The Global Geodetic Observing System (GGOS; Plag and Pearlman 2009) aims at a target accuracy equivalent to 1 mm and a stability of 0.1 mm/year for the parameters of a reference frame used for high-precision navigation on Earth and in space, aviation, monitoring of geophysical processes and many more applications. This GGOS goal corresponds to approximately 30 µas for PM and 2 µs for UT1. Although ERPs may be determined with high accuracy in a combined estimation of all space-geodetic techniques for the past, Belda et al. (2017) and Angermann et al. (2020), amongst others, have pointed out that low-latency ERP estimations, i.e. from several weeks ago until the present day, do not meet the requirements for high-precision real-time applications. Furthermore, both studies found that the treatment of station positions of the underlying reference frame has a direct impact on the EOP estimation. For example, the parameterisation of seasonal station variations (as it is done for ITRF2014; Altamimi et al. 2016) or the consideration of geophysical processes, such as non-tidal loading, via models (as it is done in DTRF2014; Seitz et al. 2022) leads to discrepancies of the resulting EOPs that are about a factor of 3-5 above the GGOS specifications for a terrestrial reference frame realisation with station positions and constant station velocities. Moreover, the transformation between the celestial and terrestrial reference frames cannot be computed with sufficient accuracy yet, stressing the importance of precise and consistent EOPs as the link between them.
Beyond retrospective and (near-)real-time applications for monitoring, positioning and navigation purposes, there is a growing interest in the relation between the Earth's rotation and global change phenomena, e.g. concerning climate change (Chen et al. 2013;Mitrovica et al. 2015;Adhikari and Ivins 2016;Śliwińska et al. 2021). Göttl et al. (2015Göttl et al. ( , 2021 and Chen et al. (2017) show that precise EOPs can contribute to improved modelling of effects related to these phenomena. In this context, the International Association of Geodesy (IAG) Inter-Commission Committee on "Geodesy for Climate Research" (ICCC) Joint Working Group (JWG) C.1 "Climate Signatures in Earth Orientation Parameters" investigates the interdependency between variations in the Earth's rotation and the Earth's climate on short and longer time scales (Poutanen and Rósza 2020a). Reliably determined retrospective and predicted ERPs are thus an important contributor to an improved understanding of the system Earth. Consequently, the GGOS Committee on Essential Geodetic Variables (EGVs) includes the ERPs in the list of observed variables that are crucial to describe key properties of the Earth (cf. Angermann et al. 2022, Sect. 2.2).

Current ERP products
The International Earth Rotation and Reference Systems Service (IERS) is responsible for providing official EOP products on an operational basis (IERS Terms of Reference, see Dick and Thaller 2020). Within the IERS, two Product Centres are responsible for the EOP generation with a focus on providing EOP products at different latencies, namely the IERS Earth Orientation Centre at l'Observatoire de Paris (Bizouard et al. 2019;Dick and Thaller 2020, Sect. 3.5.1) for the final long-term EOP data, provided in the IERS 14C04 EOP time series, and the IERS Rapid Service/Prediction Centre (RS/PC) at the US Naval Observatory (Stamatakos et al. 2012;Dick and Thaller 2020, Sect. 3.5.2) for the rapid data including predictions up to 1 year into the future, made available in the IERS Bulletin A (on a weekly basis) as well as the finals.daily 1 (contains the results from the Bulletin A combination and prediction on a daily basis; Fig. 1).
The IERS 14C04 EOP series results from a combination of different space-geodetic technique solutions which are routinely provided by the IAG Scientific Services dedicated to the organisation and data analysis of the respective technique: the International VLBI Service for Geodesy and Astrometry (IVS; Nothnagel et al. 2017), the International Laser Ranging Service (ILRS; Pearlman et al. 2019), the International GNSS Service (IGS; Johnston et al. 2017) and the International DORIS Service (IDS; Willis et al. 2015).
The combination procedure for the computation of the long-term IERS 14C04 series applied at the IERS Earth Orientation Centre involves several processing steps, which are described in detail by Bizouard et al. (2019). First, the  Centres for the  ITRS, ICRS, long-term and  low-latency EOP products and  the relationships between them. Moreover, the input data and latencies are shown. AAM atmospheric angular momentum; LS + AR least squares + auto-regressive modelling processing-related formal standard deviations of individual input solutions are re-scaled to a realistic level through an intra-technique comparison. Afterwards, a yearly rescaling factor per contributing solution is determined and applied. The technique-specific solutions are aligned to the 'guide series' derived from the latest International Terrestrial Reference Frame (ITRF) and International Celestial Reference Frame (ICRF) solutions by estimating long-term linear trends in conjunction with applying discontinuities at several epochs. However, a full consistency of the EOPs with the ITRF and the ICRF is not achieved yet (Belda et al. 2017;Kwak et al. 2018). The goal to ensure full consistency between the EOPs and the current ITRF and ICRF is, however, difficult to achieve with the combination method applied for generating the IERS 14C04 series, as this combination is done at the solution level (i.e. the level of parameters). If a combination at the normal equation (NEQ) level would be applied that includes all three parameter groups-namely the terrestrial reference frame (TRF), the celestial reference frame (CRF) and the EOPs-much higher consistency could be achieved. Investigations in this direction have already started (e.g. Bachmann and Thaller 2017;Kwak et al. 2018).
The GGOS Bureau of Products and Standards (BPS) compiled an inventory of standards and conventions used to generate the IAG products like the TRF, the CRF and EOPs (Angermann et al. 2020). The inventory summarises the present status concerning the product generation, identifies current deficiencies and gives recommendations for improvements. Section 4.3 of the BPS inventory focusses on the official IERS EOP products, including the procedures for the computation of the final IERS 14C04 EOP series and the realisation of the rapid and predicted ERPs (IERS Bulletin A). The inventory recommends to reduce the inconsistencies among the technique-specific input data, to apply more rigorous combination methods and to increase the update frequency of EAM data for prediction. Besides a redefined algorithm for the IERS C04 EOP product for alignment with new ITRF or ICRF realisations (cf. Fig. 1; Dick and Thaller 2017, Sect. 3.5.1), modifications of the algorithm (e.g. Dick and Thaller 2020, Sect. 3.5.1) have also resulted in significant irregular retrospective changes in the time series (e.g. Pavlis and Kuzmicz-Cieslak 2017).
Besides the IERS Product Centres, also NASA's Jet Propulsion Laboratory (JPL) provides three different ERP products on a routine basis, the most recent release available at the time of this study being the SPACE2018, the COMB2018 and the POLE2018 time series (Ratcliff and Gross 2019). The JPL SPACE2018 time series was combined at the solution level from SLR/LLR, VLBI, and GNSS by applying a Kalman filter approach and provides daily ERP values over a time span between 1976 and 2019. Thereby, the space-geodetic series have been adjusted by considering biases, outliers and realistic uncertainty levels before combination and filtering. The JPL COMB2018 time series additionally contains optical astrometric PM and UT1 measurements determined by the Bureau International de l'Heure (BIH;cf. Li 1985). ERP values are provided with a daily resolution over a time span between 1962 and 2019, therewith covering the full era of space-geodetic observations. The JPL POLE2018 series provides long-term PM values combined from space-geodetic techniques as well as optical measurements. The ERPs have an approximately monthly sampling covering a time span from 1900 to 2019. Dill et al. (2020) transform various ERP time series into the EAM domain and compare them to geophysical fluid models in terms of root mean square (RMS) differences and explained variances. Thereby, JPL COMB2018 provides good results for both the polar and the axial EAM components. For the polar components, the solution performs comparable to IERS 14C04. For the axial component related to UT1, JPL COMB2018 outperforms all other solutions in terms of smallest RMS differences (reduced by more than 50% as compared to IERS 14C04) and highest explained variances (nearly 4 times as large as compared to IERS 14C04) for the short periods between 2 and 8 days. Consequently, JPL COMB2018 is considered one of the benchmark products concerning ERP accuracy. However, the JPL ERP products are final ERP products covering a limited time span, and updates are only provided annually.
In the following, we use both the IERS 14C04 and the JPL COMB2018 time series as references for the quality assessment of our novel approach. Moreover, we investigate the consistency between final, rapid and predicted ERPs by hindcast comparisons between IERS Bulletin A (in the form

Framework and general concept
In 2017, the Navigation Support Office of the European Space Agency (ESA) initiated the project "Independent Generation of Earth Orientation Parameters" (ESA-EOP). It aimed to develop and implement an optimised strategy to combine routinely available VLBI, SLR, GNSS and DORIS data into a single ERP product that contains final, rapid and predicted ERPs with smooth transition, consistent with the latest available terrestrial and celestial reference frame realisations ITRF2014 and ICRF2 (Fey et al. 2015). The final objective of this is the provision of an independent and operational service for EOP products (Schoenemann et al. 2020). Thereby, the approach was required to account for different latencies of the input data. The study was carried out by an international consortium of the Deutsches Geodätisches Forschungsinstitut at the Technical University of Munich (DGFI-TUM; lead), TUM's chair of Satellite Geodesy, the Federal Agency for Cartography and Geodesy (BKG), the Section 1.3 Earth System Modelling of GFZ Potsdam, and the VLBI group of TU Wien. While a combination at the solution level is highly sensitive to (potentially unknown) constraints applied in the processing at the side of the IAG Scientific Services' Analysis Centres (ACs) and would practically always require a correction of the contributing solutions for empirical bias and drift parameters, a combination at the observation level as the most rigorous way of combination would require common processing of all techniques in a single software package. An advantage of this would be that one could apply harmonised standards for data screening (limited by technique-specific requirements and processing standards). However, this approach would exclude the products already provided on a routine basis by the ACs. The most flexible approach for the problem at hand is a combination at the NEQ level where the technique-specific data processing with independent software packages is uncritical. Combination at the NEQ level can be considered equivalent to a rigorous combination at the observation level, provided that one applies common background models and ERP parameterisation (e.g. Seitz et al. 2015;Bloßfeld 2015). We consider it uncritical that a data screening can only be performed for each technique separately before generating the NEQs, assuming that the observations finally used for generating the NEQs do not contain outliers. However, the ERPs in the NEQs provided by the ACs are not harmonised yet across different spacegeodetic techniques, which is equal to the non-parameterised ERPs in a technique-specific NEQ being fixed to their a priori values.
We base our studies on input data provided by individual ACs to the IAG Scientific Services in order to demonstrate the potential improvement for ERP estimation and prediction that is already possible with products available at present. Table 1 shows the input data used for this study and the respective latencies. As outlined in Sect. 4.1, all input data have been provided in the form of normal equations in the SINEX v2.02 format (cf. IERS Message No. 103, 2006). For GNSS Finals, GNSS Rapids and SLR, we rely on the IGS and ILRS contribution, respectively, of ESA-ESOC. The GNSS Rapids have been produced daily and, for this study, also provided in the form of NEQs in the SINEX v2.02 format. For DORIS data, we rely on the IDS contribution of the Space Geodesy Research Group (GRGS), while VLBI data have been provided by DGFI-TUM (VLBI 24-h sessions) and BKG (VLBI-Intensive sessions).
As can be seen, the VLBI 24-h sessions, SLR and GNSS Finals are available only with a latency of several days to weeks, and weekly DORIS SINEX data are provided quarterly only. Consequently, a final ERP combination that routinely includes DORIS data is possible in a reprocessing scenario only but not operationally with short latencies. Rapid ERP estimates must rely on GNSS Rapids and VLBI Intensive (VLBI-INT) sessions. Thereby, GNSS Rapids are processed with < 1 day latency but contain a smaller station network than the GNSS Finals ). Furthermore, one has to take into account the specific VLBI-INT session configurations involving only two or three stations, limiting the sensitivity of the observations to the determi-nation of UT1 with an accuracy that is reduced compared to VLBI 24-h sessions (Hellmers et al. 2019; for the usage of different VLBI session types in EOP determination, also refer to Malkin 2020, Chapt. 5). As a result, the very shortterm ERP prediction suffers from the uncertainty in the latest rapid ERPs, even when introducing EAM forecasts based on geophysical models of the atmosphere, ocean and continental hydrology (Dill and Dobslaw 2010;Dobslaw and Dill 2018). Consequently, our study puts special focus on exploiting the GNSS-derived LOD to bridge gaps between VLBI sessions and to stabilise the UT1 offsets in the rapid combination (see Sect. 4.2 for more details). Belda et al. (2017) have shown that fixing the station coordinates yields the best agreement between TRF-related ERPs and ERPs estimated from VLBI. As our goal is to estimate ERPs consistent with the current ITRF and ICRF solutions, the TRF is consequently realised by fixing a stable network of datum stations to their a priori coordinates. Stations that are not included in the ITRF2014 or stations that are affected by unmodelled discontinuities or drifts are reduced (often called "pre-eliminated") from the NEQ beforehand. Reducing parameters means to remove them from the NEQ system in a way that they remain free but are not explicitly estimated (cf. Bloßfeld 2015). For the present study, this means to reduce all stations that are not included in the ITRF2014, as their a priori coordinates may be inaccurate, while the reference stations from the ITRF2014 have undergone a screening for discontinuities. The same holds for VLBI sources that are not included in the ICRF2. Dill et al. (2020) evaluate different processing options for a combination of ERPs at the NEQ level, comparing the resulting time series to EAM excitation functions derived from geophysical fluid models. Their comparison confirms the importance of introducing precise and consistent a priori TRF coordinates in the ERP determination, which significantly improves the agreement between the combined ERPs and the EAMs. It is important that one has to take care of accurate a priori coordinates also on the long run, as the extrapolated ITRF coordinates lose reliability over time. Thus, we underline the necessity to continuously monitor the reliability of station coordinates and appropriately exclude stations that are not suitable for datum realisation.

Combination strategy
The procedures have been implemented within a software package called ESA-EOP Software. It realises time series of final (FIN), rapid (RAP), and predicted (PRE) ERPs. As shown in Fig. 2, the FIN + RAP + PRE ERP solution thereby contains a final part composed of weekly solved FIN NEQs (latency > 2 weeks), a so-called continuous part composed of the last FIN NEQ plus daily added RAP NEQs solved in a single adjustment (latency < 1 day), and an EAM-based prediction part up to 90 days into the future.
The combination of the space-geodetic techniques is performed at the NEQ level, resulting in FIN (weekly) and RAP (daily) combined NEQs, whereby the ERPs are stacked at the noon epochs of each day. The temporal resolution of the combined multi-technique NEQs matches the limitations given by the input data, as the SLR and DORIS ACs provide weekly NEQs while the GNSS and VLBI ACs provide daily and session-wise NEQs, respectively. A continuous transition from the last FIN ERP to the RAP part of the time series is ensured by two steps: First, the ERPs in the combined NEQs are transformed from offset/drift at the noon epochs to a continuous piece-wise-linear parameterisation at the midnight epochs; for an explanation how this can be done also for the VLBI NEQs, please refer to, e.g., Thaller et al. (2009). Second, the daily RAP combined NEQs are added to the last FIN NEQ via stacking of the ERPs at the day boundaries, resulting in a NEQ we call "continuous". The transition from the last combined RAP ERPs to the prediction is realised seamless by using the last ERP of the "continuous" time series as the initial value for the prediction up to 90 days into the future. The prediction algorithm thereby relies on deterministic signals derived from the FIN + RAP ERP time series (all previous FIN solutions plus the most recent "continuous" solution) as well as EAM analysis and prediction data. All in all, the implemented algorithm realises a seamless FIN + RAP + PRE ERP time series from the past until 90 days into the future. Figure 3 summarises the data flow within the entire procedure.
The technique-specific input data (see Table 1) are provided in the form of NEQs in the SINEX v2.02 format and transformed into a software-internal binary format. Normal equation systems are thereby given in the form of a normal equation matrix N, the right-hand side y and the weighted sum over observed-minus-computed l T P ll l and are solved in  Koch 1999). In the following, we outline the procedure shown in Fig. 3; for detailed graphical visualisation, please refer to "Appendix A". For a detailed description of the mathematical representation of the processing steps as implemented within the developed software, as well as of their (non-)equivalences at the observation or solution levels of the Gauß-Markov model, the reader may refer to Bloßfeld (2015, Sects. 2.1 and 2.2). The most important formulae are summarised in "Appendix B", with Eq. B4 providing the general functional model for parameter transformations. While Eq. B5 is relevant for harmonising the a priori values of parameters before stacking, or to transform given a priori values to a pre-defined reference, Eq. B7 describes the procedure for, e.g., parameter epoch transformations and reparameterisations; Eq. B8 is used to correct the functional model of a parameter for unmodelled biases.
In the first processing steps, we prepare the techniquespecific NEQs for the inter-technique combination. This technique-specific preprocessing comprises reducing ("preeliminating") technique-specific parameters like biases and satellites' initial state vectors and transforming the a priori values of the station coordinates to the ITRF2014. For GNSS, the functional model of the LOD parameter is corrected to take into account the GNSS LOD bias (cf. Sect. 4.2). For VLBI 24-h and VLBI-INT sessions, the ERPs are epoch transformed from the mid-session epochs to noon epochs for stacking with the other techniques ( Fig. 4a; "Appendix B.1"). Thereby, the UT1 and LOD parameters are used together with a regularisation according to Petit and Luzum (2010). In the case of the FIN combination scenario, the station coordinates of all techniques are epoch transformed to the mid Reparameterisation of ERPs from offset/drift representation to two offsets at the day boundaries. Note: the operations apply for regularised UT1 and LOD parameters only epoch of the particular week, whereas in the RAP combination scenario, they are epoch transformed to the mid (noon) epoch of the particular day. Finally, non-datum stations are reduced from the NEQ systems and the datum is realised by fixing the remaining stations to their a priori coordinates. For VLBI 24-h and VLBI-INT sessions, the CRF is realised by fixing the source coordinates to their a priori values in the ICRF2. While the main part of the celestial pole motion can be modelled well, the nutation offsets parameterised in the VLBI 24-h sessions remain free to cover remaining unmodelled effects, most notably the free core nutation (Dehant et al. 2003). The resulting preprocessed technique-specific weekly (FIN scenario) and daily (RAP scenario) NEQs are used as input to the actual inter-technique combination.
The inter-technique combination is performed weekly in the case of the FIN scenario and daily in the case of the RAP scenario. The FIN scenario combines the preprocessed weekly NEQs from GNSS Finals, SLR, VLBI (24-h and VLBI-INT sessions) and DORIS, whereby we apply an equal weighting of all technique-specific contribution. Common ERPs of all techniques are stacked at the noon epochs and subsequently reparameterised to piece-wise linear offsets at the day boundaries, meaning a reparameterisation from PM offsets and drifts, UT1, and LOD at each noon epoch to pairs of PM and UT1 offsets at each of the neighbouring day boundaries ( Fig. 4b; "Appendix B.2"). Also in this step, UT1 and LOD are regularised according to Petit and Luzum (2010). Afterwards, the parameters at the day boundaries between consecutive days are stacked in order to realise a daily ERP polygon covering the entire week. The reason for stacking common ERPs at noon epochs before reparameterising them to the day boundaries is to avoid singularities due to the fact that SLR and DORIS NEQs do not parameterise all ERPs. Finally, the weekly combined FIN NEQ is solved. All available FIN solutions are used to generate a FIN ERP time series by equating the ERPs at week boundaries at the solution level in the form of a weighted mean taking into account the derived variances. By this, we avoid the necessity to set up a normal equation system that increases in size over time and past solutions also do not change any more.
The RAP scenario is limited to combining the preprocessed daily GNSS Rapid and VLBI-INT NEQs. The ERPs of both techniques are stacked at the day's noon epoch and subsequently reparameterised to piece-wise linear offsets at the day boundaries. Afterwards, the last available weekly FIN NEQ and all subsequent daily RAP NEQs are combined into one common "continuous" NEQ, stacking the ERPs at the day boundaries. Missing RAP combined NEQs due to missing VLBI-INT sessions are bridged with GNSS Rapids, stressing the importance of appropriate treatment of the GNSS-derived LOD (Sect. 4.2). Finally, the "continuous" NEQ is solved and all previous FIN solutions together with the "continuous" solution are used to generate the FIN + RAP ERP solution, which is the input for the prediction algorithm. Thereby, the ERPs at week boundaries between weekly FIN solutions and between the FIN and "continuous" part of the time series are again equated at the solution level in the form of a weighted mean.
Finally, the 90-day prediction is calculated by a two-step procedure. First, we derive deterministic signals by transforming the FIN+RAP ERP time series into geodetic angular momentum (GAM) via the Liouville equation and computing the differences between the GAM the EAM analysis data. These deterministic signals are then extrapolated by least squares harmonic analysis (LS) and autoregression (AR) modelling until the end of the 6-day EAM forecasts. After the end of the 6-day EAM forecast period, the full EAM signal (EAM analysis and forecast plus extrapolated GAM residuals) is again extrapolated by LS + AR until the end of the 90-day forecast period, and transformed into the ERP domain via the Liouville equation. A detailed description of the LS + AR prediction procedure and the combination of FIN + RAP ERP with EAM (AAM, OAM, HAM, SLAM analysis + 6-day forecasts) data is described by Dill et al. (2019) and Dobslaw and Dill (2019). The result is a complete FIN + RAP + PRE ERP time series.
The output of the combined and predicted ERP time series is written in a format resembling the IERS C04 format. For the output of the UT1 parameter, the software allows choosing between the official UT1 − UTC and a step-free UT1−TAI representation. To avoid the issue of leap seconds when stacking ERPs at day boundaries, our combination procedure internally uses the UT1 − TAI representation instead of the UT1 − UTC representation reported in the SINEX input files.

Treatment of GNSS and GNSS-derived LOD
The LOD bias describes an offset of the LOD estimates derived by satellite techniques w.r.t. a "physical" LOD. For GNSS, deficiencies in the modelling of orbit perturbations, for example, the solar radiation pressure, systematically impact the motion of the orbital nodes. This leads to biased GNSS LOD estimates which can distort the combined UT1 and LOD estimates (Ray 1996;Senior et al. 2010;Mikschi et al. 2019). Depending on the constellation (number of GPS, GLONASS, and Galileo satellites) and the parameterisation (e.g. daily solutions as used in the present study, or 3-dayssolutions), the behaviour of the GNSS-derived LOD can be characterised by seasonal or abrupt changes. For a homogeneous inter-technique combination as well as for bridging data gaps due to missing VLBI sessions, exploitation of GNSS-derived LOD information with appropriate treatment of the GNSS LOD bias (Kammeyer 2000; Dick and Thaller 2020) is indispensable.
Retrospectively, the GNSS LOD bias can be determined via a comparison of GNSS-derived LOD estimates with VLBI-derived LOD. Based on comparisons between GNSSonly solutions and the official IERS 14C04 time series, we determined a GNSS LOD bias of approximately − 20 µs (cf. Table 2, Fig.5) valid for the GNSS contribution of ESA-ESOC, which is GNSS repro2 until April 2014 and the standard GNSS Final IGS submissions of ESA-ESOC after that date. This time series is based on GPS and GLONASS. The results in Sect. 5 are based on a combination of these GNSS Final data in the FIN combination and the corresponding GNSS Rapid data in the RAP combination, whereby the LOD bias is applied as a correction to the functional model of LOD in the GNSS NEQs ("Appendix B", Eq. B8) before the combination with the other space-geodetic techniques. For the GNSS Final data, it is important to note that the phase centre variation (PCV) models applied are igs08.atx until January 2017 and igs14.atx after that date. While this switch in the PCV models would have a significant impact on the realisation of the scale of a terrestrial reference frame, it is considered uncritical for the realisation of the ERPs (Rebischung and Schmid 2016).

Validation
We assessed the quality of the ERP time series computed with our homogeneous combination and prediction methods in two steps. In the first step, our FIN ERP solution was validated against the external final ERP products IERS 14C04 and JPL COMB2018 (Sect. 5.1). In the second step, the internal consistency of ERP solutions was investigated by so-called hindcast experiments. Thereby, the RAP and PRE

Combination strategy
This subsection summarises the results of a direct comparison between final ERP solutions, taking into account the following time series: • IERS 14C04 (Bizouard et al. 2019).
• FIN ERP solution (this study). Figure 6 shows the differences between the ERP products IERS 14C04 and JPL COMB2018. As a measure for the consistency between these two solutions, the RMS deviations are given, which are about 19 µs for UT1 and 66 µas and 50 µas for x pole and y pole , respectively. It has to be noted that neither IERS 14C04 nor JPL COMB2018 can be considered as an error-free reference for the external accuracy assessment of any ERP time series obtained with our new approach. Figure 7 shows the comparison between our FIN ERP solution and the IERS 14C04 series. For UT1 (upper panel of Fig. 7), the RMS deviation is 22 µs, which is in the range of the IERS 14C04 versus JPL COMB2018 comparison. Similar results (lower panel of Fig. 7) are also found for polar motion with RMS deviations of 58 µas for x pole and 50 µas for y pole . Taking into account the reported formal uncertainties of different series (see Table 3), the obtained differences between those series are not significant.  Figure 8 shows the comparison between our FIN ERP solution and the JPL COMB2018 series. With 12 µs, the RMS deviation of UT1 is the smallest of all three pair-wise comparisons carried out. For PM, we find RMS deviations of 85 µas and 74 µas for x pole and y pole , respectively, which are larger than the differences between our FIN ERP solution and IERS 14C04.
Given the three pair-wise differences shown above, the results indicate that our FIN ERP solution and IERS 14C04 are closest to each other for both PM components, whereas for UT1, we observe the best agreement between our FIN ERP solution and JPL COMB2018. These fits suggest that the quality of our FIN ERP solution is comparable to the two other series that are typically considered as state of the art.
For the year 2018 (the last year contained in the JPL COMB2018 solution), Table 3 summarises the average formal errors of the FIN ERP solution and the reported uncertainties of the external ERP solutions compared. Since the formal errors strongly depend on the assumed stochas- tic behaviour of the input data, they are only of limited significance to assess the performance of the underlying combination procedure but have to be viewed in combination with the input data used. In our case, the variance propagation may lead to a too optimistic level of the formal errors, as GNSS yields systematically small formal errors due to known mismodellings (Schön and Kutterer 2007;Schön and Brunner 2008).

Internal validation by hindcast experiments and availability of VLBI 24-h sessions
The quality of rapid and predicted ERPs was assessed by determining their consistency with the corresponding final ERP solution. For each comparison, we generated ensembles of approximately 500 daily hindcast experiments for the years 2018 and 2019. Thereby, each hindcast experiment realises a final, rapid and predicted ERP time series as can be assumed to be available at the respective epoch, taking into account the input data latencies. The time series from our experiments are FIN + RAP + PRE ERP time series calculated for each day of the hindcast period, the prediction based on the EAM data available at the respective epoch. We calculated the two following ensembles of hindcast experiments (Table 4): -The first series of hindcast experiments (H1) has been performed under realistic assumptions, meaning that the RAP ERPs are combined from GNSS Rapids and VLBI-INT data. -In the second series of hindcast experiments (H2), the RAP combination has been performed under the assumption that, in addition to VLBI-INT sessions, also VLBI 24-h sessions were provided with < 1 day latency.
The comparison contrasts the FIN + RAP + PRE ERP time series from the hindcast experiments against the FIN ERP solution covering the full hindcast period. As a reference for the comparison, we performed similar hindcast evaluations contrasting the IERS Bulletin A rapid/predicted ERP product (in the form of finals.daily) available at the epoch against the IERS 14C04 final ERP product. Figure 9 shows the forecast-horizon-related mean RMS values resulting from the comparison of each hindcast ensemble with the respective final ERPs. Thereby, the plots distinguish between the combined part (shown is the combination for a forecast horizon from − 30 to 0 days; note: in the following, non-positive forecast horizons refer to the combined ERPs while positive forecast horizons refer to the actual prediction) and the predicted part (shown is the forecast horizon from 0 days to + 30 days of the 90-day prediction; note: day 0 refers to the last set of combined ERPs, i.e. the set of initial values for prediction) of the time series.
The validation reveals that our approach performs well for the PM parameters (Fig. 9, uppermost and middle panels); the consistency between RAP and FIN ERPs outperforms that between the IERS products until a forecast horizon of −1 day. The improved agreement between the PM components from IERS Bulletin A and IERS 14C04 for day zero was not anticipated. We assume that it might be related to a high weight of GNSS-derived PM in both the final days of IERS Bulletin A and the IERS 14C04 product. The predictions from our study again outperform IERS Bulletin A for a forecast horizon between 2 and 30 days for x pole and between 2 and 23 days for y pole . This improvement demonstrates the advantage of a prediction based on consistently combined FIN and RAP ERPs with smooth transition. The additional availability of VLBI 24-h sessions with low latency, however, would only have a minor effect on the determined PM as GNSS observations are highly sensitive to PM and usually also dominate the combined PM estimates.
For UT1 (Fig. 9, lowest panels), both scenarios H1 and H2 perform significantly better than IERS Bulletin A, demonstrating the importance of a homogeneous combination of final and rapid ERPs by appropriately exploiting GNSS-derived LOD. The availability of VLBI information for the determination of UT1, however, turns out to be crucial for improving the rapid part of the time series: As scenario H2 demonstrates, currently achievable accuracies of rapid ERP combinations could be improved by over 50 % Fig. 9 Mean RMS values of the ERP time series for the forecast horizon − 30 days to 0 days ("continuous" vs. FIN ERPs, left panels) and 0 to + 30 days (PRE vs. FIN ERPs, right panels) and mean RMS values of IERS Bulletin A (rapid and predicted ERPs) against IERS 14C04 (final ERPs). For the hindcast scenarios, please refer to Table 4 if low-latency processing of VLBI 24-h sessions could be performed, yielding a much more accurate start value for the prediction.
In all cases, it is visible that our FIN + RAP + PRE ERP solutions involve a direct transition between the FIN and RAP part of the time series, as both parts are fully consistent at a forecast horizon of − 30 days. The comparison of Bulletin A against IERS 14C04, on the other hand, reveals the inconsistencies between two products generated by two different institutions with their particular processing choices and inhomogeneous input ERP data. Table 5 lists the obtained mean RMS values of predicted ERPs from our study and IERS Bulletin A for different forecast horizons up to 90 days into the future. At forecast horizons in the range between 5 and 10 days, the achieved accuracy of our prediction is almost 50% better in comparison with the IERS product. For longer-term predictions over 30 days, the IERS Bulletin A provides better results for PM, whereas our prediction performs better for UT1 at all forecast horizons. The reason for this is most probably a different tuning of the prediction algorithms: our algorithm being especially optimised for short-to medium-term predictions up to

Conclusion and outlook
The developed strategy combines space-geodetic observations at the normal equation level and delivers a consistent set of ERPs. It offers seamless transition, starting from final processing using archived observations extending back into the past for decades, continuing with rapid processing of the most recent data, and finishing with predictions based on deterministic signals derived from the combined ERPs, EAM analysis and EAM forecast data for up to 90 days into the future. We could demonstrate that our combination strategy performs at a level comparable to IERS 14C04 and JPL COMB2018, the ERP products commonly seen as the stateof-the-art benchmarks. ERP predictions over several days into the future outperform the IERS Bulletin A by almost 50 % in terms of consistency with the corresponding final ERPs, proving the benefit of consistent processing of the final and rapid ERPs following a common approach. The results demonstrate that a refined rigorous ERP combination leads to more reliable and accurate combined retroand prospective ERPs. Reducing inconsistencies inherent to a combination at the solution level, our combination approach at the NEQ level, however, still fully relies on the information content provided in the technique-specific NEQs. Consequently, it would be valuable to further homogenise the technique-specific processing strategies (i.e. parameterisation and background models; cf. Angermann et al. 2020) of the IAG Scientific Services in combination with thorough monitoring of the GNSS LOD bias.
To further improve the consistency of the combination, we consider it highly beneficial if all technique-specific NEQs parameterised either one full set of ERPs per day, or piece-wise linear PM offsets and UT1, independent from the sensitivity of the technique-specific observations to the parameters. Examples for this are the NEQs provided for both VLBI 24-h and VLBI-INT sessions that parameterise PM offsets and rates, although the usually 1-h long VLBI-INT observations are not sensitive to the terrestrial pole and LOD.
As VLBI is the only space-geodetic technique giving direct access to UT1, the technique is essential to determine the full set of ERPs. To date, VLBI 24-h sessions are usually covering two consecutive days starting at around 17:30 h UTC, while VLBI-INT sessions usually cover a period of 1 h starting around 07:30 h or 18:30 h UTC. 2 As most VLBI Analysis Centres parameterise one set of ERPs at the session mid epoch, this means that the parameters must be extrapolated to either mid-day or the day boundaries. To eliminate extrapolation errors, we consider it valuable to schedule the VLBI 24-h sessions from 0:00 h to 24:00 h UTC. Another measure to optimise the VLBI contribution to a regularly-spaced series of UT1 would be to organise the VLBI data processing of the standard VLBI sessions (i.e. scheduled from 17:30 h UTC onwards) in a way that all sessions are directly parameterised as piece-wise linear offsets at midnight epochs. For the current session scheduling across 2 days, this would lead to two intervals of piece-wise linear ERP estimates, i.e. the first interval covering 17:30 until 24:00 h UTC, and the second interval covering 24:00/00:00 until 17:30 h UTC.
To conclude, we formulate several recommendations that may contribute to improving accuracy and consistency of next-generation ERP products: (a) The standards for the processing of technique-specific observations should be further harmonised across all IAG Scientific Services. Besides using common background models, this should include parameterising either one set of PM offsets and rates, UT1, and LOD per day, or piecewise linear offsets at the day boundaries. Both final and, if available, rapid data should be provided routinely in the form of NEQs in the commonly adopted SINEX format. (b) An automated low-latency processing should be envisaged not only for the VLBI-INT sessions but also for the VLBI 24-h sessions. This would significantly improve the achievable accuracy of UT1 in rapid ERP solutions. (c) For an improved estimation of UT1, it is essential to apply constellation-dependent LOD biases for the GNSS input series. The applied LOD bias can be kept constant as long as the constellation-dependent mean values (without seasonal signals) vary within a range of about ±5 µs (the observation accuracy of VLBI; cf. Dick and Thaller 2020, Sect. 3.5.2). (d) To date, the ITRF and the ICRF reference frames are not fully consistent, thus introducing inconsistencies into the combination via VLBI. A homogeneous realisation of both CRF and TRF consistently with retrospective final ERPs is thus crucial for realising highly accurate rapid and predicted ERPs (Kwak et al. 2018). In this context, the International Union of Geodesy and Geophysics (IUGG, Resolution No. 3, 2011;see Drewes et al. 2012) and the IAG (Resolution No. 2, 2019; see Poutanen and Rósza 2020b) both urge highest consistency between ITRF, ICRF and EOPs for all future realisations.
So far, our combination approach does not imply a filtering of the data, as it is applied, e.g., for the JPL COMB2018 series. Consequently, one focus of future studies might be to investigate the benefits of complementing a combination at the NEQ level by an appropriate filtering strategy, combining the advantageous properties of the approach presented here with, e.g., the benefits of the JPL COMB2018 approach.
Concerning the prediction, we emphasise that the algorithm developed within the ESA study currently includes EAM forecasts only for up to 6 days into the future. Deterministic forecasts from ECMWF and other numerical weather prediction centres typically extend to 10 or 14 days into the future, so that more information is potentially available. Moreover, seasonal forecasts extending over several months ahead of time are now being regularly produced in an ensemble setting that could also be valuable for ERP prediction. It would even be possible to extend the forecast range beyond 90 days as soon as sufficient urgent user needs for longer predictions will be made known.   In the case of UT1 and LOD, the relationship between non-regularised parameters p,ṗ and regularised parameters p,p is with r andṙ being the regularisation of p andṗ, respectively, according to Chapter 8.1 of Petit and Luzum (2010). The regularisation is a nonlinear part of the functional model, resulting in the functional model for the regularised transformation to be inserted into Eqs. B4 and B7.

B.3 ERP re-parameterisation to piece-wise linear representation
The transformation is described by and leads to the functional model for the linear transformation In the case of UT1 and LOD, the relationship between non-regularised and regularised parameters is according to Eq. B11 and the regularisation is calculated according to Chapter 8.1 of Petit and Luzum (2010).
The regularisation is a nonlinear part of the functional model, resulting in the functional model for the regularised transformation to be inserted into Eqs. B4 and B7.