Evaluating Processing Choices for the Geodetic Estimation of Earth Orientation Parameters With Numerical Models of Global Geophysical Fluids

Different Earth orientation parameter (EOP) time series are publicly available that typically arise from the combination of individual space geodetic technique solutions. The applied processing strategies and choices lead to systematically differing signal and noise characteristics particularly at the shortest periods between 2 and 8 days. We investigate the consequences of typical choices by introducing new experimental EOP solutions obtained from combinations at either normal equation level processed by Deutsches Geodätisches Forschungsinstitut at the Technical University of Munich (DGFI‐TUM) and Federal Agency for Cartography and Geodesy (BKG), or observation level processed by European Space Agency (ESA). All those experiments contribute to an effort initiated by ESA to develop an independent capacity for routine EOP processing and prediction in Europe. Results are benchmarked against geophysical model‐based effective angular momentum functions processed by Helmholtz Centre Potsdam GFZ German Research Centre for Geosciences (ESMGFZ). We find, that a multitechnique combination at normal equation level that explicitly aligns a priori station coordinates to the ITRF2014 frequently outperforms the current International Earth Rotation and Reference Systems Service (IERS) standard solution 14C04. A multi‐Global Navigation Satellite System (GNSS)‐only solution already provides very competitive accuracies for the equatorial components. Quite similar results are also obtained from a short combination at observation level experiment using multi‐GNSS solutions and SLR from Sentinel‐3A and Sentinel‐3B to realize space links. For ΔUT1, however, very long baseline interferometry (VLBI) information is known to be critically important so that experiments combining only GNSS and possibly SLR at observation level perform worse than combinations of all techniques at normal equation level. The low noise floor and smooth spectra obtained from the multi‐GNSS solution nevertheless illustrates the potential of this most rigorous combination approach so that further efforts to include in particular VLBI are strongly recommended.


Introduction
The orientation of the solid Earth with respect to the celestial reference frame needs to be precisely known for a number of applications including ground-based astrometric observations, communication with satellites, including probes in deep space, and also Global Navigation Satellite System (GNSS) nowadays used for the positioning of sometimes rapidly and even autonomously moving objects on the ground or in the air. Space geodetic techniques such as GNSS at permanent stations, very long baseline interferometry (VLBI), satellite laser ranging (SLR), or Doppler Orbitography and Radio-Positioning Integrated by Satellite (DORIS) provide information about time variations in the position of the terrestrial pole (polar motion), the phase angle of Earth's rotation ΔUT1 expressed as the difference between universal time (UT1) and the coordinated universal time (UTC), and the celestial pole offsets (nutation). Those five (time-variable) parameters are conventionally referred to as Earth orientation parameters (EOPs). The drift parameters related to each of these EOP can be also determined by the space geodetic techniques. Therein, ΔLOD plays an important role related to the spin rate of the Earth.
Due to the advent of more precise sensors, denser measurement networks, and the availability of (at least partly) redundant observation techniques, the precision of space geodesy has improved over the most recent decades. Commonly, the available sensor data are combined into intratechnique EOP solutions in a least squares sense to arrive at best possible solutions with minimal errors. A number of intratechnique EOP solutions is subsequently combined by various approaches to arrive at one single EOP time series. However, in view of the high internal precision of the individual techniques, it becomes increasingly important to enforce consistency among the different techniques to avoid the introduction of spurious artifacts. This includes in particular all aspects of the realization of the terrestrial reference system. Similar attention should be devoted to geophysical background models required to process individual observations like, for example, solar radiation pressure effects on individual satellites, or ocean tide models including ocean tidal loading that affect space geodetic observations in numerous and typically highly systematic ways. A more rigorous way for the combination of the individual space-geodetic technique solutions would be the combination at the normal equation (NEQ) level of the Gauss-Markov model before solving for EOP. Ideal from a theoretical perspective would be the combination at observation level using one single software with identical parametrizations and background models to invert the observations from all techniques at once. So far, no publicly available EOP time series is applying any of the latter two approaches.
Polar motion and ΔLOD are governed mainly by terrestrial processes associated with the redistribution of masses of the near-surface geophysical fluids. Variations in ΔLOD are largely dominated by zonal tropospheric winds (Salstein, 1993), whereas atmospheric surface pressure and ocean dynamics are rather equally important for the excitation of high-frequency polar motion variations (Ponte & Ali, 2002). On seasonal time scales, large-scale variations in terrestrial water storage are dominant (Chen et al., 2012). On decadal-to-centennial periods, prominent contributors to polar motion are the low-frequency changes in the continental ice masses (Adhikari & Ivins, 2016), whereas ΔLOD is also affected by core-mantle coupling effects (Holme & De Viron, 2013).
The quality of available models of global geophysical fluids relevant for the excitation of Earth orientation changes has increased tremendously in the more recent past. Atmospheric reanalyses produced by Meteorological Services like the European Centre for Medium-Range Weather Forecasts (ECMWF) are now routinely available (Dee et al., 2011). Particularly, the mass component estimates of ocean and land hydrosphere models have benefited from the availability of time-variable gravity field obtained with the GRACE mission (Göttl et al., 2019;Śliwińska et al., 2020). We therefore consider it nowadays as a viable option to use a geophysical model data set as the reference against which different geodetic combination time series are compared. Although geophysical models cannot be considered as error-free, typical error sources of geodetic space techniques-like dependencies of the solar radiation pressure modeling on the satellite's beta angle (elevation of the Sun above the orbital plane) or spacecraft geometry-are not inherent in geophysical models, and therefore should become visible in such a comparison.
The paper is structured as follows: We collect three of the most commonly used EOP series that were calculated from a combination of different geodetic space techniques, and additionally introduce four experimental EOP combination series processed specifically for this study within a project of the European Space Agency to improve EOP (section 2). Subsequently, we derive so-called geodetic excitation functions (GAM) out of the EOP that can be readily contrasted against geophysical effective angular momentum (EAM) functions (section 3). Time series comparisons are provided in terms of root-mean-square differences, Taylor plots, and explained variances for different frequency bands (section 4). Since largest differences among the geodetic solutions are found for periods shorter than 8 days, we specifically discuss spectra for those highest frequencies (section 5). The paper closes with a discussion of the differences found in the geodetic solutions and some recommendations for future improvements in the processing of combined geodetic EOP solutions.
For completeness, we note that the celestial pole offsets are largely governed by gravitational attraction of different bodies of the solar system. Only a very tiny fraction of the nutation is caused by (seasonally modulated) diurnal tides in oceans and atmosphere that additionally deform the solid Earth by means of surface loading (Nastula &Śliwińska, 2020). Albeit formally a part of the set of Earth Orientation Parameters, we entirely disregard celestial pole offsets in this study.

Selected EOP Time Series
The Earth Orientation Center of the International Earth Rotation and Reference Systems Service (IERS) at Paris Observatory is the official provider (Bizouard, 2020) of daily estimates of polar motion and ΔUT1. The excess length of day ΔLOD that is related to the Earth's rotation spin rate equals the difference of consecutive UT1-UTC estimates.

C04-08: IERS 08C04
The combination solution IERS 08C04 aligned to the ITRF2008 (called C04-08 in the reminder of this paper) results from a combination of intratechnique EOP series obtained from GNSS, VLBI, SLR, and DORIS (Gambis & Bizouard, 2009). One or two representative series from each technique are considered for the pole coordinates. For ΔUT1, the whole set of VLBI series available from the International VLBI Service for Geodesy and Astrometry (IVS) is taken into account, because no space geodetic techniques other than VLBI is able to determine ΔUT1 in an absolute sense.
The intratechnique EOP series entering into the combination are made compatible by rescaling the formal uncertainties and by realigning to both the International Celestial Reference Frame (ICRF) and the International Terrestrial Reference Frame (ITRF). Pole coordinates are smoothed by an epoch-dependent Vondŕak-Filter (Vondrak, 1977) and are interpolated linearly to equidistant daily epochs. The trend of the ΔUT1 series derived from GNSS and SLR is aligned to the trend of ΔUT1 obtained from VLBI. The final series are again smoothed by Vondrák-filtering to remove spurious variations likely introduced by the applied numerical procedures. Vondrák smoothing coefficients can be found at page 4 of the C04 description document (ftp://hpiers.obspm.fr/iers/eop/eopc04_08/C04.guide.pdf). Since C04-08 refers to the axis of the nowadays outdated ITRF2008, a slow degradation of the overall accuracy can be expected in particular for epochs in the year 2009 and later.

C04-14: IERS 14C04
The EOP combination procedure applied at Paris Observatory has been thoroughly upgraded to calculate a new series IERS 14C04 (Bizouard et al., 2017), called here C04-14. This EOP solution is realigned to the most recent ITRF, thereby also improving the numerical combination procedure by the introduction of more realistic weights for the intratechnique solutions. Updated Vondrák smoothing coefficients are reported in Table 3 in (Bizouard et al., 2019). Pole coordinates of C04-14 are now consistent with ITRF2014, whereas nutation offsets and ΔUT1 are aligned to the ICRF2 and ICRF3 before and after the year 2019, respectively. The series C04-14 has been reprocessed back until 1962 and is officially recommended by the IERS. It is updated 2 times per week, with an average latency of about 30 days. Differences to the previous solution C04-08 are as large as 50 μas in polar motion and 5 μs in ΔUT1, and are primarily related to the selected terrestrial reference frame.

JPL-Comb2018
Earth Orientation Parameters are also processed at the Jet Propulsion Laboratory (JPL) of the National Aeronautics and Space Administration (NASA) in a manner that is fully independent from IERS. The so-called JPL-Comb2018 solution utilizes tracking data from Lunar Laser Ranging (LLR), the Global Positioning System (GPS), VLBI, SLR, and historic optical astrometric observations by means of a Kalman Filter approach (Ratcliff & Gross, 2019). Rotational variations caused by solid Earth (Yoder et al., 1981) and ocean tides (Kantha et al., 1998) were removed from the ΔUT1 values prior to the combination and added back afterward.
As the individual space geodetic techniques might use their own realizations of the terrestrial reference system when solving for EOP; for example, EOP(IGS) 00 P 03 for the GNSS solutions provided by the International GNSS Service (IGS), both bias-rate corrections and uncertainty scale factors were determined for each single-technique EOP time series. Each individual series was compared to a combination of all other remaining series to estimate those parameters individually for each technique. The procedure was repeated iteratively until convergence among all considered single-technique solutions had been reached.
It should be noted that updates to this series are only published annually. For routine applications JPL provides daily updates including short-term predictions by additionally incorporating rapidly available observations such as the GPS and AAM forecasts from NCEP (https://keof.jpl.nasa.gov).

Experimental Solutions by DGFI-TUM and BKG
The European Space Agency (ESA) is currently working toward establishing an independent capacity for calculating EOP out of multiple space geodetic data products processed within its Navigation Support Office (OPS-GN) at the European Space Operations Center (ESOC). An external team is currently being tasked with the development of a new combination software suitable for routine EOP estimation and prediction. This group consists of scientists from Deutsches Geodätisches Forschungsinstitut (DGFI-TUM) at the Technical University of Munich, Federal Agency for Cartography and Geodesy (BKG), Chair of Satellite Geodesy at the Technical University of Munich, Research Group Advanced Geodesy at the Technical University of Vienna, and the Earth System Modelling group at the Helmholtz Centre Potsdam GFZ German Research Centre for Geosciences (ESMGFZ). The work is in particular based on previous experience gained at DGFI-TUM as an IERS ITRS Combination Center (Seitz et al., 2012), and at BKG which is operating the IVS Combination Center jointly together with DGFI-TUM (Bachmann et al., 2016).
All input data to the combination software are provided in terms of technique-specific NEQs given in the solution-independent exchange format (SINEX) by ESA with the exception of the VLBI solutions (BKG). Before combination, the technique-specific NEQs undergo a set of preprocessing steps. Whereas GNSS, SLR, and DORIS already contain EOP parameterized at noon epochs, the VLBI-based EOP need to be transformed from session midpoints to the nearest noon epochs. The functional model of the ΔLOD parameter in the GNSS NEQs is expanded in order to account for a potential ΔLOD bias. In this study, we apply a fixed correction value of −20 μs which is based on an analysis (with respect to C04-14) of the ESA ESOC GPS+GALILEO LOD time series between 2016 and 2019. Daily GNSS NEQs and session-wise VLBI NEQs are then accumulated to weekly technique-specific NEQs in order to match the weekly resolution of SLR and DORIS. The TRF datum for all techniques is kept by fixing all station coordinates to their a priori values, which ensures consistency between the estimated EOP and the a priori reference frame (Belda et al., 2017).
The combination of the weekly technique-specific NEQs to a common weekly NEQ is performed by summing up all NEQs with equal weights. Thereby, all technique-specific EOP at noon epochs are stacked to combined EOP at noon epochs. Parameterized are pole offsets, pole rates, ΔUT1, and ΔLOD. Each daily set of EOP at noon is transformed to the respective day boundaries as a pair of midnight offsets at 0 and 24 hr UTC, taking into account the effect of tidal deformation on the Earth's rotation in the transformation of ΔUT1 and ΔLOD according to the IERS Conventions (Luzum & Petit, 2012). As described in Chapter 8.1 of the conventions, all periods from 5 days to 18.6 years are considered for regularization. Afterward, the pole offsets and ΔUT1 at the day boundaries between consecutive days are stacked. As VLBI is the only space geodetic technique that allows for the direct observation of ΔUT1, this procedure ensures that gaps between VLBI sessions are bridged with ΔLOD information from the satellite techniques. Thus, the combined NEQ system is invertible without any further EOP constraints. After inversion, weekly solutions with full sets of EOPs at the day boundaries (eight different epochs) are obtained. A time series of consecutive daily EOP estimates is subsequently generated by stacking the EOP values at the week boundaries at solution level, that is, by calculating a weighted mean of the estimates. With that software and general processing strategy, the following two experiments E1 and E2 were performed.

Experiment E1: Combination at NEQ Level
For Experiment E1, we use NEQs of GNSS and SLR solutions processed at the Analysis Center (AC) ESOC as regular contribution to the IGS, and to the International Laser Ranging Service (ILRS), respectively. In addition, 24-hr VLBI solutions are used from the IVS AC at DGFI-TUM, whereas VLBI Intensive solutions are taken from the IVS AC at BKG. Station coordinates as given in each intratechnique NEQ are not changed in this experiment. The main problem arising from this treatment of the routine products as is is that the ITRF realization to which the coordinates are referred changes over time, so the results have to be taken with care. Moreover, the NEQs provided by the IAG services do not necessarily contain station coordinates that are fully consistent with the ITRF2014, as technique-specific realizations of this TRF are used.

Experiment E2: Combination at NEQ Level Aligned to ITRF2014
In order to improve the consistency of the datum definition across all techniques, in the second experiment (E2) the station coordinates from ITRF2014 stations have been transformed to the ITRF2014 datum in advance. This procedure reduces datum inconsistencies for all stations given in the ITRF2014, but introduces some inconsistencies within the networks between ITRF2014 and non-ITRF2014 stations. However, these inconsistencies remain neglectable in the beginning of the processed period as the vast majority of sites processed is contained in ITRF2014. Later on, the inconsistencies become more relevant, as more stations not considered in the ITRF2014 are added especially to the GNSS network. Apart from the transformation of the a priori values before fixing the station coordinates, the combination approaches of experiments E1 and E2 are identical.

Experimental Solutions by ESA
We hypothesize that consistency of the contributions from the different geodetic space techniques is a key element to achieve a best possible EOP accuracy. To achieve that goal, ESOC reprocessed archived observation data from the International Doris Service (IDS), IGS, and ILRS in a single homogenized solution (Otten et al., 2012) by using their own software NAvigation Package for Earth Orbiting Satellites (NAPEOS). This approach allows for the most rigorous combination of IDS, ILRS, and IGS reference frames by using space ties. ESA is aiming for combining all space geodetic techniques on observation level (GNSS, SLR, DORIS, and VLBI). However, to understand the impact of the different observation types, the solution is carefully extended by adding only one new observation type at a time. We use in this article two intermediate solutions.

Experiment E3: Multi-GNSS Solution as Contribution to the Third IGS Reprocessing Campaign
The experiment E3 used in this study is identical with the ESA contribution to the third reprocessing campaign of the International GNSS Service (IGS). The EOP solution is based on the daily analysis of undifferenced pseudorange and carrier phase observations of 150 globally distributed multi-GNSS IGS tracking stations considering on average 31 GPS and 24 GLONASS satellites as well as, starting from January 2014 also Galileo satellites. Initially only four Galileo satellites were included, but the number increased to 24 until December 2018. As the data from the three constellations are jointly processed, with common receiver clocks estimates allowing for corresponding intersystem biases, the solutions can be considered as combined at the observation level with highest consistency. In view of a full set of EOPs, it is important to emphasize that especially VLBI is missing in experiment E3. Thus, ΔUT1 cannot fully be determined.

Experiment E4: Combination of GNSS and SLR at Observation Level
We introduce also a very recent solution that combines GNSS observations with tracking data of Sentinel-3A and Sentinel-3B as low Earth orbiters for this space link. Both GNSS and SLR observations to those satellites are considered. The data are rigorously combined at observation level. So far just 12 months of data from experiment E4 were completed so that a full evaluation of this series by means of model-based EAM is not possible. Therefore, we will discuss E4 in section 5 only. Please note that ΔUT1 can be expected to be determined similarly poorly as in experiment E3.

Effective Angular Momentum Functions
Changes in the orientation of the solid Earth are conveniently studied by applying the principle of conservation of angular momentum in the whole Earth system including its surrounding fluid layers. Relevant are both the instantaneous mass distribution (matter terms) and the relative angular momentum changes associated to winds and currents (motion terms). Globally integrated angular momentum changes are multiplied with empirically derived parameters to account for the actual rheology of the Earth including, for example, the anelasticity of the mantle, the partly decoupled rotation of the core, and the associated equilibrium response of the oceans (Brzeziński, 1992;Gross, 2007). It is important to note that in contrast to EOP time series, EAMs are free of the dominating Chandler wobble in the equatorial components.
Globally integrated changes in angular momentum of each of the subsystems can be described by effective angular momentum functions (EAM) derived from individual global numerical models. Customarily, those contributions are named as atmospheric angular momentum (AAM), oceanic angular momentum (OAM), and hydrological angular momentum (HAM). The additional effect of a time-variable barystatic sea level in response to a net transfer of water mass from the land into the ocean is sometimes assumed to be part of the OAM, but sometimes also kept separated and labeled as sea-level angular momentum (SLAM).

ESMGFZ: Geophysical Model-Based EAM
The various geodetic solutions will be evaluated against a model-based EAM time series provided by the Earth System Modelling group at Deutsches GeoForschungsZentrum (ESMGFZ). The daily updated nontidal EAM data are given in terms of dimensionless effective angular momentum functions of the matter and motion terms individually for the Earth's subsystems atmosphere, ocean, and terrestrial water storage. The underlying mass redistribution for atmospheric surface pressure is taken from the European Centre for Medium-Range Weather Forecasts (ECMWF), ocean bottom pressure and vertically integrated ocean currents are simulated with the Max-Planck Institute for Meteorology Ocean Model (MPIOM) (Jungclaus et al., 2013), and terrestrial water storage is simulated with the Land Surface and Discharge Model (LSDM) (Dill, 2008). Please note that contributions of the 12 most prominent tidal frequencies have been removed from atmosphere and ocean. Additional contributions arising from major earthquakes (Chao & Gross, 1995;Yun, 2019), electromagnetic jerks at the core-mantle boundary (Ron et al., 2019), or glacial processes in the continental ice sheets (Mitrovica & Wahr, 2011) present in the geodetic observations are not covered by this model-based data set. As the geophysical models do only represent mass variations and mass exchange but provide no access to the absolute atmospheric, oceanic, and terrestrial water masses, a long-term mean (2003-2014) has been already removed from the EAM products. Further information on the product is provided via the web page https://esmdata.gfz-potsdam.de:8080/repository and in the product description document (Dobslaw & Dill, 2018).

Geodetic Angular Momentum
To obtain excitation functions out of observed EOP, the Liouville equation with pole coordinates p = p 1 + ip 2 and complex Chandler frequency c = 2 (1 + i∕2Q)∕T c is deconvoluted (Wilson & Vicente, 1990) to transform pole coordinates into so-called geodetic angular momentum functions (GAM) = 1 + i 2 . We use a Chandler period of T c = 434.2 days with a damping of Q = 100, which is consistent with the parameterization of the rotational deformation applied in the model-based EAM calculations. The axial component 3 follows from GAM are available for every day since 1962. Those GAM should be therefore understood as the excitation required to change Earth orientation in a way as it is observed by space geodesy. Effects of long-period tides 10.1029/2020JB020025 were removed from ΔLOD as recommended in the IERS conventions (Luzum & Petit, 2012) to make 3 comparable to the nontidal EAM from ESMGFZ.
As an introductory example, we show time series of GAM derived from JPL-Comb2018 together with the sum of model-based EAM functions from ESMGFZ ( Figure 1). We note that model-based EAM explain almost 90% of the intra-annual signal in 3 , which is related to the dominance of seasonal variations in zonal tropospheric winds that are very well captured by present-day atmospheric reanalyses. For the equatorial components, residuals are much larger (≈50%) with both strong high-frequency variability and a distinct long-term trend. The equatorial components are rather sensitive to mass distributions in polar regions with both strong variability in the wind-driven ocean dynamics and slow mass loss of ice sheets and glaciers. Nevertheless, a considerable fraction of the signal seen by JPL-Comb2018 is explained by the model-based EAM so that it is sensible to use the geophysical model as a reference to evaluate the different geodetic solutions.

Time Series Analysis
GAM series are calculated according to section 3.2 from all EOP series introduced in section 2. Except for experiment E4, all series are available to us with daily sampling from January 2009 to April 2019. EAM are taken as sum of AAM, OAM (both sampled from 3 hr sampling to the daily epochs of GAM), HAM, and SLAM. A third-order Butterworth filter with varying cut-off periods is applied to split all time series into three frequency bands: (1) 2-8 days, (2) 8-20 days, and (3) 20-100 days. In addition, also the (4) combined band of 2-100 days, and the (5) unfiltered series that includes all periods above 2 days are considered. We calculate various metrics commonly applied in time series analysis to quantify the correspondence of GAM and EAM. In particular, we use root-mean-square differences (RMSD), standard deviations (STD), the Pearson correlation coefficient (CORR), and explained variances (EXVAR).
RMSDs quantify the residual variability after subtracting ESMGFZ EAM from any of the GAM series, reduced by their mean offset over the analyzed period ( Figure 2). For the periods above 8 days, we find very consistent results across the six GAM series considered. The only exception is the experiment E1, which has 5% higher RMSD in 1 . Differences among the geodetic series are more pronounced at the highest frequencies: For the pole, E1 fits rather poorly to ESMGFZ when compared to the other solutions. For ΔLOD, both E1 and C04-08 have the largest misfit, whereas both experiments E2 and E3 are even slightly better than C04-14. In all components, JPL-Comb2018 provides the best fit to the model, and the largest margin with respect to the competing geodetic series is found in the third component.
To properly interpret the RMSD, it should be viewed in relation to the standard deviations of the two time series involved. It should be noted that the RMSD can be readily calculated from STDs and CORR according to where indices t and ref denote the time series to be tested and the reference time series, respectively. That relation equals the law of cosines where STD ref and STD t are the length of the sides of a triangle, and CORR t, ref referring to the cosine of the angle between those sides. Hence, RMSD t, ref is the length of the third side of the triangle vis-a-vis to the correlation angle. Equation 3 therefore provides a geometrical relationship between the different metrics that can be used to display all three metrics jointly within a so-called Taylor diagram (Taylor, 2001). The Taylor diagram shows the agreement of any time series with a reference series not only by means of the RMSD, but provides a synopsis with the corresponding STD and CORR values. In principle, we are looking for results with a low RMSD, a STD similar to the reference series, and a high CORR coefficient.
In the following, we present Taylor diagrams that display results not only for the different GAM series (each by a separate color) but also for the different filters applied (each by a separate marker). For every category, the STD of the geophysical model-based time series ESMGFZ is given at the axis of abscissa as the reference point. The Euclidean distance from the reference point to the marker (STD t , CORR t ) of an individual series gives the RMSD t that is equal to the values given in the bar plots of Figure 2.
For both equatorial components (Figure 3, top row), we generally find a good correspondence of all GAM series with the modeled EAM. Results for 20-100 days (stars) are very close to each other, and also the results  for 8-20 days are quite similar for all six geodetic series considered. For the shortest periods below 8 days (squares), we find a substantially larger spread: C04-08 and C04-14 are still very close to each other, with slightly smaller RMSD and slightly higher correlation for the more recent series from IERS. JPL-Comb2018 has a notable smaller STD than C04, which nevertheless does not always lead to a smaller RMSD misfit. We also find a huge reduction in STD for E2 when compared to E1: Since both experiments only differ in the treatment of the station coordinates (as given in the SINEX files for E1; taken from ITRF2014 where possible for E2), this result clearly underlines the importance of precise a priori coordinates for the determination of EOP.
We further note that experiment E3 always has the smallest STD from all geodetic time series considered. We recall that this is a multi-GNSS solution only and VLBI, SLR, and DORIS observations are not included in this experiment. We nevertheless note that correlation and also RMSD are already quite competitive with respect to the other geodetic series. This indicates that pole coordinates are indeed very well determined from GNSS information alone. It is important to recall the (relatively) good performance of E3 might arise from the fact that all geodetic solutions except E3 have to deal with different parametrizations for the station positions adopted by the various Analysis and Technique Center which have a direct impact on the EOP solutions (Bloßfeld et al., 2014). For completeness, we also present the results for the band 2-100 days (pluses) and the unfiltered series (dots). The results basically reflect the findings of the weekly band and do not need to be reiterated here.
For the axial component ( Figure 3, bottom row), we find again very consistent results across all geodetic series for the lower frequencies and significant scatter only for the shortest periods of 2-8 days. For this component, C04-14 is a substantial improvement over the older series C04-08 with much reduced STD of the series, leading to both a smaller RMSD and a higher CORR with the geophysical EAM. This improvement is mirrored by the difference between E1 and E2, highlighting again the importance of a consistent terrestrial reference frame for EOP estimation. E3 has again the smallest STD of all series considered, but CORR and RMSD are much worse than experiment E2, thereby strongly underlining the well-known importance of VLBI for the determination of ΔUT1 and consequently ΔLOD. The best results in this comparison are obtained with JPL-Comb2018, where a similarly small STD is connected with CORR and small RMSD, indicating that a good compromise has been found in this series to suppress high-frequency noise while retaining the relevant short-period signals. As for the equatorial components, the results for the other frequency bands are also included in the plots for completeness, but do not provide additional insights.
As an additional evaluation metric not captured by Taylor plots, we define the explained variance (EXVAR) as with STD 2 err as the variance of the unexplained signal, that is the difference between the time series and its reference. Note that this quantity is also sometimes called coefficient of determination in the statistical literature. For identical time series, EXVAR equals 100%, and for time series not fitting at all it might even become negative.
For the polar motion excitations 1 and 2 , EXVAR reaches values between 30% and 75% depending on the period band considered (Figure 4). Differences among the six geodetic solutions are very small apart from the shortest periods between 2 and 8 days. Here, four series have a similar level of EXVAR for both 1 and 2 , whereas experiment E1 has very small and barely positive values only. As the a priori station coordinates were kept as given in the intratechnique NEQs and it is not mandatory that the technique-specific realizations of the terrestrial reference system are aligned to each other, station coordinates in E1 might differ among the techniques. Those differences in the station coordinates were eliminated in E2, which consequently does not contain anymore such spurious high-frequency signals that almost entirely mask the real geophysical signal contained in the geodetic observations. Best results in this comparison are again obtained by JPL-Comb2018.
In the axial component 3 , the largest spread between the geodetic solutions is also found at the highest frequencies. C04-08 and E1 have largely negative explained variances. C04-14 and E2 reveal significant improvements, with E2 outperforming C04-14 by a substantial amount. It is interesting to note that the experiment E3-the multi-GNSS solution-is also already outperforming C04-14 and lags only slightly behind E2. The best performance, however, is found again with JPL-Comb2018.

Spectral Analysis
We calculate amplitude spectra for all GAM time series and their residuals against the model-based EAM from ESMGFZ. For the longer periods of the equatorial components 1 and 2 , the residuals are dominated Figure 5. Amplitude spectrum of geodetic angular momentum time series GAM derived from different EOP solutions and model-based EAM from ESMGFZ, for 1 , (top), 2 (middle), 3 (bottom). For better readability the individual spectra were smoothed (five-point boxcar) and shifted by 0.5 · 10 −8 for 1 and 2 and 0.2 · 10 −10 for 3 . Excitation functions are unitless.
by a peak at 13.66 days not present in the EAM and possibly related to errors in the fortnightly tides (Ray et al., 2017). For the highest frequencies between 2 and 8 days, the spectra of the residuals against EAM differ substantially (Figures 5,top and 5,middle). We note very high variability and several significant peaks in both C04-08 and also E1. Those peaks somewhat reduce for C04-14 and E2, but remain much larger than in JPL-Comb2018, where the energy found at the highest frequencies is even lower than in the geophysical model. The experiment E3 instead has very little energy at the highest frequencies, which is between 2 and 3 days even smaller than in JPL-Comb2018. This is indeed interesting, since GNSS information with high temporal resolution has been ingested by the solution.
Results are quite similar also for the axial component 3 (Figure 5, bottom). Prominent peaks are found in E1 and E2 at 7 days, which corresponds conspiciously to the chosen weekly NEQ accumulation interval. Less prominent peaks are also visible at the associated overtones of 3.5 and 2.3 days. A similar characteristic is also seen in C04-08, but disappeared almost entirely in C04-14, which is known to suppress high-frequency variations by a strong smoothing algorithm. JPL-Comb2018 and also E3 instead do not contain such prominent peaks. For the highest frequencies, JPL-Comb2018 and E2 are approximately at the same level as ESMGFZ. It should be noted, however, that VLBI 24-hr sessions are performed regularly twice a week (Mondays and Thursdays), which might contribute to the identified systematic. Moreover, no smoothing is applied in experiments E1 and E2. In contrast, the amplitude spectra of E3 calculated only from GNSS information reveals much smaller variability at those subweekly periods than predicted by the geophysical model, thereby clearly suggesting that important variability is not captured by the selected observing system configuration.
We also present here results from a preliminary combination of GNSS and SLR at observation level (experiment E4), which is only available to us over 12 months from July 2018 to June 2019 so that it could not be readily included into the analysis presented above. From the comparison of the residuals against experiment E3 ( Figure 6) it becomes obvious that the combination at observation level closely follows the multi-GNSS solution with no obvious systematic differences. Differences between E4 and E3 are more than one magnitude smaller than the RMS of E3 to our reference ESMGFZ. Deviations of E4 from E3 are also smaller than the deviations to other EOP series, for example, JPL-Comb2018. However, because of the limited time span, we cannot conclude how far the addition of SLR improves the multi-GNSS EOP solution E3. Nevertheless, the results are generally encouraging and should further motivate ESA to extend the combination to a longer time span and include other geodetic techniques in order to allow for an in-depth analysis of EOP obtained from this most rigorous combination approach.

Summary and Conclusions
Three publicly available time series of terrestrial pole coordinates and ΔUT1 estimates are augmented for this study by four experimental EOP series processed by DGFI-TUM, BKG, and ESA that are all transformed into time series of geodetic angular momentum for contrasting against global geophysical fluid models. All geodetic series reveal very similar variations for periods longer than a week, but show systematic differences among each other at periods between 2 and 8 days. We therefore conclude that individual processing choices during the geodetic data analysis significantly affect the resulting EOP, in particular in the shortest periods.
A comparison against geophysical model-based excitation functions from ESMGFZ by means of various metrices (standard deviations, correlations, RMSDs, explained variances) documents the relative improvements achieved by the IERS with the transition from C04-08 to C04-14. The comparison also documents the superior quality of JPL-Comb2018, even though it has to be kept in mind that the solution processed at JPL is not updated routinely but instead processed at once for a fixed period of time. JPL-Comb2018 therefore should be regarded as the target accuracy that should be aimed at with any EOP solution processed in an operational setting.
The new experimental EOP solutions processed by DGFI-TUM and BKG in an operational setting agree well to the results obtained for the publicly available series. GAM from a combination of data from different geodetic space techniques at normal equation level that utilizes a priori coordinates as given in the SINEX files show spurious high-frequency signals and corresponding poor fits to the geophysical EAM. In the underlying EOP series the inconsistencies in the TRFs lead to high-frequency artifacts together with several jumps followed by short-lasting drifts that cannot be removed easily when combining EOP at the solution level. The quality of EOP obtained from a NEQ level combination drastically increases when a priori coordinates are harmonized to a consistent common reference frame. This solution generally even outperforms C04-14, thereby demonstrating that the operational setting with input data from independent sources combined at normal equation level, developed by DGFI-TUM and BKG, results in highly competitive EOP estimates. Furthermore, it demonstrates that a combination at normal equation level is preferable to a combination at parameter level.
From a theoretical perspective, a combination at observation level that utilizes space ties among the different geodetic techniques would be ideal for the processing of EOP. Available to us are a multi-GNSS solution processed by ESA as a contribution to the third reprocessing campaign of the IGS as well as preliminary results from a combination of Sentinel-3A and Sentinel-3B with GNSS processed at ESOC. EOP from these solutions are characterized by exceptionally low noise at the highest frequencies which lead to the best fit with the geophysical model for the equatorial components among all operational geodetic series considered. For the axial component, information from VLBI that is still missing in those solutions leads to a degraded quality with respect to the results of a NEQ level combination (including VLBI R1-, R4-, and intensive sessions) with ITRF2014 a priori coordinates. Nevertheless, the achieved results for the pole are very promising and efforts should be expedited to also include VLBI and other techniques into this solution type.
It should be emphasized that no additional smoothing has been applied to the EOP series specifically processed for this study. Spurious effects identified in either the time series or the spectral analysis as presented will now be analyzed further in order to identify possible causes for those artifacts. This might include the consequences of the selected accumulation length of 7 days; the regular schedule of the 24-hr sessions (which might be assessed by focusing on the epochs of the CONT campaigns, where significantly more VLBI data are available); or the impact of certain background model choices including the treatment of subdaily tidal signals.
On a final note, the demonstrated ability to reliably identify consequences of individual processing choices on geodetic data products with the geophysical model-based angular momentum functions demonstrate the tremendous improvement in accuracy in those models achieved in the more recent past. For low-frequency signals that allow for the accumulation of geodetic observations over long periods of time and thus abundant redundancy, geodetic estimates might be still safely regarded as a reference to benchmark numerical models against. For the higher frequencies with less observations and a relatively higher impact of systematic errors, however, it would be prudent to evaluate for each individual case if information readily provided by numerical models that incorporate information from various nongeodetic sources could be advantageously combined with data from space geodesy to finally arrive at products with better external accuracies.