<figure><img src="image_1.png" /></figure>

HAL
open science

# ITRF2020: an augmented reference frame refining the modeling of nonlinear station motions

Zuheir Altamimi, Paul Rebischung, Xavier Collilieux, Laurent Métivier, Kristel Chanard

► To cite this version:

Zuheir Altamimi, Paul Rebischung, Xavier Collilieux, Laurent Métivier, Kristel Chanard. ITRF2020: an augmented reference frame refining the modeling of nonlinear station motions. Journal of Geodesy, 2023, 97, 10.1007/s00190-023-01738-w. insu-04155715

HAL Id: insu-04155715

https://insu.hal.science/insu-04155715

Submitted on 8 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ITRF2020: an augmented reference frame refining the modeling of nonlinear station motions

Zuheir Altamimi<fnref n="1" />,<fnref n="2" />* · Paul Rebischung<fnref n="1" />,<fnref n="2" /> · Xavier Collilieux<fnref n="1" />,<fnref n="2" /> · Laurent Métivier<fnref n="1" />,<fnref n="2" /> · Kristel Chanard<fnref n="1" />,<fnref n="2" />

Received: 23 September 2022 / Accepted: 12 April 2023 / Published online: 19 May 2023

© The Author(s) 2023

## Abstract

To better describe the shape of the constantly deforming Earth's surface, the ITRF2020 is provided as an augmented terrestrial reference frame that precisely models nonlinear station motions for both seasonal (annual and semi-annual) signals present in the station position time series and Post-Seismic Deformation (PSD) for sites impacted by major earthquakes. Reprocessed solutions in the form of station position time series and Earth Orientation Parameters using the full observation history provided by the four space geodetic techniques (DORIS, GNSS, SLR and VLBI) were used as input data, spanning 28, 27, 38 and 41 years of observations, respectively. The ITRF2020 long-term origin follows linearly with time the Earth's Center of Mass (CM) as sensed by SLR, based on observations collected over the time span 1993.0–2021.0. We evaluate the accuracy of the ITRF2020 long-term origin position and time evolution by comparison to previous solutions, namely ITRF2014, ITRF2008 and ITRF2005, to be at the level of or better than 5 mm and 0.5 mm/yr, respectively. The ITRF2020 long-term scale is defined by a rigorous weighted average of selected VLBI sessions up to 2013.75 and SLR weekly solutions covering the 1997.75–2021.0 time span. For the first time of the ITRF history, the scale agreement between SLR and VLBI long-term solutions is at the level of 0.15 ppb (1 mm at the equator) at epoch 2015.0, with no drift. To accommodate most of ITRF2020 users, the seasonal station coordinate variations are provided in the CM as well as in the Center of Figure frames, together with a seasonal geocenter motion model. While the PSD parametric models were determined by fitting GNSS data only, they also fit the station position time series of the three other techniques that are colocated with GNSS, demonstrating their high performance in describing site post-seismic trajectories.

**Keywords** Reference systems · Reference frames · Nonlinear station motions · ITRF

## 1 Introduction

Solid Earth science and operational geodesy applications require the availability of a long-term global terrestrial reference frame, as a unique standard reference, to ensure inter-operability and consistency of geodetic products and to adequately exploit various measurements collected by ground-based sensors, or via artificial satellites. The International Terrestrial Reference System (ITRS) and its realization, the International Terrestrial Reference Frame (ITRF) are recommended by a number of international scientific organizations, such as the International Association of Geodesy and the International Union of Geodesy and Geophysics (IUGG 2007, 2019).

The ITRF is a collaborative international effort, starting with the investments of national mapping agencies, space agencies and research groups in implementing and operating geodetic observatories, followed by Data Centers (DCs) who collect and archive geodetic observations, and then by Analysis Centers (ACs) and Combination Centers (CCs) of the four space geodetic techniques that contribute to the ITRF computation: Doppler Orbitography and Radiopositioning

<figure><img src="image_3.png" /></figure>

Zuheir Altamimi
zuheir.altamimi@ign.fr

Paul Rebischung
paul.rebischung@ign.fr

Xavier Collilieux
xavier.collilieux@ensg.eu

Laurent Métivier
laurent.metivier@ign.fr

Kristel Chanard
kristel.chanard@ign.fr

<footnote n="1">¹ Institut de physique du globe de Paris, CNRS, IGN, Université de Paris Cité, 75005 Paris, France</footnote>

<footnote n="2">² IGN, ENSG-Géomatique, 77455 Marne-la-Vallée, France</footnote>
Integrated by Satellite (DORIS), Global Navigation Satellite Systems (GNSS), Satellite Laser Ranging (SLR), and Very Long Baseline Interferometry (VLBI). These techniques are organized as scientific services and are part of the structure of the International Association of Geodesy (IAG) and known by the International Earth Rotation and Reference Systems Service (IERS) as Technique Centers (TCs): the International DORIS Service (IDS; Willis et al. 2016), the International GNSS Service, formerly the International GPS Service (IGS; Johnston et al. 2017), the International Laser Ranging Service (ILRS; Pearlman et al. 2019), and the International VLBI Service (IVS; Nothnagel et al. 2017).

By precisely combining individual reference frame solutions of the four techniques, the ITRF greatly benefits from their specific strengths while integrating all of them into a single reference frame with a unique origin, scale and orientation.

Each ITRF solution is demonstrated to be of higher quality than past versions, as new and extended data are used, but also thanks to improved analysis strategies at both the technique and ITRF combination levels.

The ITRF2020 provided the opportunity to improve the modeling of nonlinear station motions and, in turn, better characterize the shape of the deforming Earth's surface. It is an augmented reference frame where the temporal station positions are modeled by piecewise linear functions, parametric models describing annual and semi-annual deformation, as well as Post-Seismic Deformation (PSD) for stations subject to major earthquakes.

The ITRS Center of the IERS, hosted by the Institut National de l'Information Géographique et Forestière (IGN) France, is responsible for the maintenance of the ITRF/ITRF and official ITRF solutions. Two other ITRS combination centers are also generating combined solutions using ITRF input data: Forschungsinstitut at the Technische Universität München (DGFI-TUM; Seitz et al. 2022), and Jet Propulsion Laboratory (JPL; Abbondanza et al. 2017; Wu et al. 2015).

The ITRF2020 products presented and discussed in this paper are available at a dedicated website: https://itrf.ign.fr/en/solutions/ITRF2020 (Altamimi et al. 2022).

## 2 ITRF2020 input data

The following two subsections describe the two data sets used in the computation of the ITRF2020 that is fundamentally based on the availability of colocation sites where two or more space geodesy instruments were or are currently operating.

### 2.1 Space geodesy solutions

The ITRF2020 benefits from reprocessed solutions, using the full history of observations of the four space geodetic techniques (DORIS, GNSS, SLR and VLBI). Table 1 lists the ITRF2020 input data of the four techniques: weekly solutions for IDS/DORIS and ILRS/SLR, daily solutions for IGS/GNSS and session-wise solutions for IVS/VLBI. It specifies the solution type (Normal Equations or Variance-Covariance), the constraints applied for the reference frame definition (free, loose or minimum constraints), and the Earth Orientation Parameters (EOPs) provided in addition to station positions. The reader may refer to (Altamimi et al. 2002, 2004; Dermanis 2001, 2004; Sillard and Boucher 2001) or to Chapter 4 of the IERS Conventions (Petit and Luzum 2010, Ch. 4) for more details regarding reference frame constraints and the minimum constraints concept in general. The submitted time series per technique are combined solutions of the individual AC solutions of that technique.

The IDS/DORIS contribution contains 1456 weekly combined solutions involving four ACs (Moreaux et al. 2022), and including 201 stations located at 87 sites.

The IGS/GNSS contribution to ITRF2020 consists of 9861 daily combined terrestrial frame (SINEX) solutions from the third IGS reprocessing campaign (repro3; Rebischung 2022). Those were obtained by combining the contributions from 10 (ACs), and include 1905 stations, of which 1344, located at 1159 sites, have been retained in ITRF2020. The observations used cover the time period 1994.0–2021.0 for GPS, 2002.0–2021.0 for GLONASS and 2013.0–2021.0 for GALILEO. One of the novelties of the IGS submission is the fact that the scale of the IGS frame is based on the calibrated Galileo satellite phase center offsets (z-PCOs) published by the European GNSS Service Center (GSC 2022) and that GPS and GLONASS satellite z-PCOs were corrected to become consistent with the Galileo satellite z-PCOs (Rebischung 2020).

The ILRS/SLR contribution consists of 244 fortnightly solutions, with polar motion and Length of Day (LOD) estimated every three days for the period 1983.0–1993.0, using LAGEOS I satellite data, and 1459 weekly solutions with daily polar motion and LOD estimates afterwards, using data acquired on LAGEOS I and II and ETALON I and II satellites (Pavlis and Luceri 2022). The main feature of the ILRS submission worth mentioning here is the careful handling of station range biases, which has an impact on the scale of the SLR frame (Pavlis et al. 2021).

The IVS contribution to the ITRF2020 is a combination of session-wise solutions of 11 different IVS Analysis Centers (Hellmers et al. 2022a, b). The number of the IVS/VLBI
session-wise solutions included in the ITRF2020 combina-
tion is 6178, spanning in total 41 years of observations,
involving 154 stations in 117 sites, but with only 14 sites
in the southern hemisphere. The majority of sessions (86%)
includes a small number of stations, ranging between 3 and 9.
842 (13%) sessions involve 10–19 stations, 8 sessions involve
20 stations, while one session exceptionally includes 31 sta-
tions. The solutions include the full set of EOPs, in addition
to station positions.

## 2.2 Local ties at colocation sites

In order to integrate and connect the networks of the four
techniques together, and subsequently form the ITRF, it is
fundamental to have a sufficient number of well-distributed
sites where two or more independent geodetic instruments
are observing and where local surveys between instrument
measuring points are available. Local surveys are usually
conducted using terrestrial measurements (angles, distances,
and spirit leveling) or/and GNSS observations. Least squares
adjustments of local surveys are performed by national agen-
cies operating ITRF colocation sites to provide relative
coordinates (local ties) connecting the instrument reference
points.

Figure 1 illustrates the full ITRF2020 network, compris-
ing 1835 stations located in 1243 sites, where about 10% of
them are colocated with 2, 3, or 4 distinct space geodetic
instruments. Over the entire ITRF2020 observation history,
we used local ties available for 101 colocation sites (10 more
sites than in ITRF2014) with two or more technique instru-
ments which were or are currently operating.

All past and newly submitted local ties to the ITRF2020
were used. After the ITRF2014 release, 32 new surveys were
conducted and their results, including local tie SINEX files
were submitted and included in the ITRF2020 computation.
All the DORIS colocation site survey data were re-adjusted
by the IGN survey department in order to generate full
SINEX files, including the most recent surveys operated at
these sites. All the local tie SINEX files used in the ITRF2020
combination are available at https://itrf.ign.fr/en/local-ties.

We recall here that the local ties used in the ITRF
combination are provided in SINEX (Solution-INdependent
EXchange) format with known measurement epochs (with
the exception of few old ties), and 80% of them are available
with full variance-covariance information.

The number of colocation sites between DORIS, SLR and
VLBI, taken by pair, is insufficient to allow for a reliable com-
bination of these three techniques alone: 12 DORIS-SLR, 10
DORIS-VLBI and 8 SLR-VLBI. However, almost all SLR
and VLBI sites are colocated with GNSS, and a significant
number of well-distributed DORIS sites are colocated with
GNSS. We count in total 253 tie vectors between GNSS and
the three other technique reference points: 77 for VLBI, 53

for SLR and 123 for DORIS. 66 additional tie vectors were
also used between old and current DORIS beacon reference
points in DORIS-only sites, but almost half of them are the
results of repeated surveys.

There are six 4-technique colocation sites in the ITRF2020:
Metsahovi (Finland), Wettzell (Germany), Badary (Russia),
Hartebeesthoek (South Africa), Greenbelt (USA), and Yarra-
gadee (Australia).

## 3 ITRF2020 analysis strategy

Among a number of other innovations and analysis improve-
ments, for the first time of the ITRF history, the ITRF2020
accumulates the time series of the four space geodetic
techniques all together, adding local ties and co-motion con-
straints at colocation sites in a single process. In preparation
for such a global stacking, intermediates steps, as detailed in
the following subsections, were performed in order to mainly
identify and reject outliers, introduce discontinuities in sta-
tion position time series and determine parametric models
for stations subject to post-seismic deformation.

### 3.1 Time series analysis

Although the full equations of the combination model imple-
mented in the CATREF Software that is used for the ITRF
computations were published in past articles, see for instance
(Altamimi et al. 2016), they are detailed in the appendix of
this stand-alone paper for the convenience of the reader.

In preparation for the global overall ITRF2020 combi-
nation described in Sect. 4, intermediate analysis steps and
stacking of the individual technique time series were con-
ducted, including the following:

* Apply No-Net-Rotation constraints to all SLR loosely constrained solutions, thus preserving the origin and the scale of each weekly solution;

* Apply No-Net-Translation and No-Net-Rotation con-
straints to IVS solutions, thus preserving the scale of
each session-wise solution provided in the form of nor-
mal equations;

* Use as they are minimally constrained solutions: this is the case of IGS daily and IDS weekly solutions where the GNSS scale and DORIS origin and scale information are preserved;

* Determine PSD parametric models using GNSS/IGS data (see Subsect. 3.4 for more details), then apply them to all technique time series at colocation sites impacted by earthquake-induced displacements/deformation, prior to the next step;

* Form per-technique long-term solutions, by rigorously stacking the time series, solving for station positions,
<table><caption>**Table 1** Summary of submitted solutions to ITRF2020</caption><thead><tr><th>TC</th><th>Data-span</th><th>Sampling</th><th>Solution type</th><th>Constraints</th><th>EOPs</th></tr></thead><tbody><tr><td>IVS</td><td>1980.0–2021.0</td><td>Daily</td><td>Normal Eq</td><td>None</td><td>PM, PMr, LOD, UT1-UTC, δX, δY</td></tr><tr><td rowspan="2">ILRS</td><td>1983.0–1993.0</td><td>Fortnightly</td><td>Var-Cov</td><td>Loose</td><td>PM, LOD</td></tr><tr><td>1993.0–2021.0</td><td>Weekly</td><td>Var-Cov</td><td>Loose</td><td>PM, LOD</td></tr><tr><td>IGS</td><td>1994.0–2021.0</td><td>Daily</td><td>Var-Cov</td><td>Minimum</td><td>PM, PMr, LOD</td></tr><tr><td>IDS</td><td>1993.0–2021.0</td><td>Weekly</td><td>Var-Cov</td><td>Minimum</td><td>PM</td></tr></tbody></table>

PM Polar motion, PMr Polar motion rate, LOD Length of day, UT1-UTC UT1 minus UTC, δX and δY: Nutation offsets

<figure><img src="image_4.png" /><figcaption>**Fig. 1** ITRF2020 network highlighting VLBI, SLR, and DORIS sites colocated with GNSS</figcaption></figure>

velocities, EOPs and 7 transformation parameters for each input solution with respect to the per-technique cumulative solution. During the iterative stacking process:

- Annual and semi-annual signals were estimated for stations with time spans greater than 2 years;

- Discontinuities were introduced in the station position time series as detailed in Subsect. 3.3. Each of the successive time intervals bounded by discontinuities is labelled with an integer called solution number (SOLN);

- For the particular case of IGS time series, in addition to the annual and semi-annual signals, periodic signals at the first 8 GPS draconitic harmonics (Ray et al. 2008) were estimated and then removed from the IGS time series.

During the stacking process of either of the individual technique time series, as well as in the global ITRF2020 combination described in Sect. 4, nearby stations, as well as the successive segments (SOLNs) of a given station, are constrained to have the same velocity and periodic signals

using the following two equations:

$$ \dot{X}^{i1} = \dot{X}^{i2} \quad (\sigma_v) \quad (1) $$

$$ F_j^{i1} = F_j^{i2} \quad (\sigma_f) \quad (2) $$

where $i1$ and $i2$ are two close-by stations or two consecutive segments (SOLNs) of a station with discontinuities. $\dot{X}^{i1}$ and $\dot{X}^{i2}$ are the corresponding velocities. $F_j^{i1} = [ax_j^{i1}, ay_j^{i1}, az_j^{i1}, bx_j^{i1}, by_j^{i1}, bz_j^{i1}]^T$ and $F_j^{i2} = [ax_j^{i2}, ay_j^{i2}, az_j^{i2}, bx_j^{i2}, by_j^{i2}, bz_j^{i2}]^T$ are the coefficients of the cosine and sine functions describing the periodic station displacements at frequency $j$ (see Appendix A for more details). $\sigma_v$ and $\sigma_f$ are the uncertainties at which the constraints (1) and (2) are imposed, which are chosen to be $10^{-6}$ m/yr and $10^{-6}$ m in case of the stacking of the individual technique time series, and as a function of the level of agreement between techniques at colocation sites in case of the global ITRF2020 combination. For some stations, different velocities were estimated before and after the discontinuity events in case of apparent slope change (i.e., no velocity equality constraints were applied in those cases). This also holds for some (not all) stations with post-seismic deformation.
## 3.2 Periodic signals

The periodic signals observed in the station position time series are due not only to geophysical or climatological phenomena, e.g., loading (van Dam and Wahr 1998), or thermal effects (Xu et al. 2017), but also due to technique systematic errors, such as the draconitic signals for satellite techniques, as evidenced by Ray et al. (2008) for GPS.

In the particular case of the ITRF computation, it is necessary to estimate the periodic signals embedded in the input data that are in the form of time series. Doing so is expected to improve the linear station velocity estimation, especially for stations with large seasonal signals and/or relatively short time series (Blewitt and Lavallée 2002). Estimating the periodic signals during the stacking process, along with station position and velocities also helps detecting discontinuities in the time series more easily, and consequently improves the offset determination.

The fact that the ITRF combination / stacking model described in Appendix A includes a 7-parameter similarity transformation between each individual daily / weekly solution and the cumulative long-term solution offers a number of options for the definition / specification of not only the long-term linear reference frame, but also the seasonal signals, in origin, scale and orientation (Collilieux et al. 2018).

Appendix C provides the equations implemented in the CATREF software that are used for the definition / specification of the periodic signals, distinguishing between minimum and internal constraints. In the stacking process, we need to specify how to separate seasonal variations of the transformation parameters from seasonal variations of the station positions. To that end, we generally chose the internal constraint approach for the origin and scale components, that is, to constrain to zero the periodic signals embedded in the time series of the corresponding transformation parameters. This approach is justified by the fact that it avoids the absorption of part of the station motions by the time series of translation and scale factor parameters. In regard to the rotation parameters, we usually chose the minimum constraint approach, i.e., to impose no net periodic rotation conditions on a set of well-distributed reference stations. More details are provided in Sect. 4 regarding the options used in the ITRF2020 global combination.

## 3.3 Handling of discontinuities

To identify the discontinuities present in the GNSS station position time series provided by the IGS, we inspected them visually, as manual offset detection tends to perform better than automated methods (Gazeaux et al. 2013). A model composed of a piecewise linear function, annual and semi-annual signals and, when needed, functions describing post-seismic deformation (see following subsection) was

adjusted to each of the time series. Position and velocity discontinuities were iteratively added to the model until no more offset nor trend change was visible in the residual time series. To help the identification of offsets and relate them to known events, we used as auxiliary data the equipment changes recorded in the station log files, as well as a list of co-seismic offsets predicted using the approach of Métivier et al. (2014). In the 1344 GNSS station position time series used in ITRF2020, a total of 2939 position and/or velocity discontinuities were identified. 1505 could be related to equipment changes or malfunction, 677 to earthquakes, 6 to volcanic activity, 185 were introduced to describe trend changes (e.g., in response to current ice melting), while the remaining 566 have unknown causes. This corresponds to one discontinuity every 6.4 years in average.

For colocation sites impacted by earthquakes, and in order to ensure consistency between the four technique time series, we adopted the same discontinuity epochs for DORIS, SLR and VLBI as for GNSS. The number of position discontinuities due to earthquakes are 37 (out of 86), 29 (out of 41), and 46 (out of 50), for DORIS, SLR and VLBI, respectively. Discontinuities unrelated to earthquakes were identified by visual inspection of the station position residual time series, resulting from the individual per-technique stacking where seasonal signals were estimated.

## 3.4 Post-seismic deformation

Initiated at the occasion of the ITRF2014 (Altamimi et al. 2016), the parametric models describing PSD for sites impacted by major earthquakes were refined using the extended ITRF2020 input data. In addition to the four PSD models introduced with ITRF2014, a fifth model is now added: (1) (Log)arithmetic, (2) (Exp)onential, (3) Log+Exp, (4) Exp+Exp, and (5) Log+Log. Each PSD model may or may not be accompanied by a velocity change.

The ITRF2020 PSD models were adjusted to the input IGS station position time series after all offsets were identified. For every earthquake followed by visible PSD, a single exponential function and a single logarithmic function were first separately tested. The function yielding the smallest weighted sum of squared residuals was retained. In case residual PSD was still visible in the residuals from this first model, a second exponential or logarithmic function was added to the PSD model, according to the same criterion.

Figure 2 displays in red the locations, extracted from the global Centroid-Moment-Tensor catalog (Dziewonski et al. 1981; Ekström et al. 2012), of the 65 earthquake epicenters that caused significant PSD in the ITRF2020 input time series, and in green the 118 impacted ITRF2020 sites. Figure 3 shows the histograms of the residuals from the trajectory models adjusted to the IGS repro3 station position time series without (in black) and with (in red) PSD. From this figure,
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Distribution of earthquake epicenters (red), extracted from the global Centroid-Moment-Tensor catalog (Dziewonski et al. 1981; Ekström et al. 2012) and ITRF2020 sites (green) impacted by Post-Seismic Deformation</figcaption></figure>

we can see that the distributions are similar whether or not PSD models are part of the trajectory models, indicating that the presence of PSD models or not has little influence on the goodness of fit of the overall trajectory models.

The corrections predicted by the GNSS-fitted models were then applied not only to the IGS station position time series, but also to the colocated stations of the three other techniques, before stacking their respective time series. We verified that the ITRF2020 PSD parametric models determined using IGS/GNSS data fit also the station position time series of the three other techniques that are sharing the same colocation sites. Thus, we inspected the post-fit residuals of these colocated station position time series to verify that no residual post-seismic deformation was visible. As an example, Fig. 4 illustrates the performance of the PSD parametric model for the colocated GNSS and VLBI stations at Tsukuba (Japan). It shows in blue the raw station position time series, in green the piecewise linear trajectories given by the ITRF2020 coordinates and in red the trajectories obtained when adding the parametric PSD model. In that figure, one can see the remarkable fit of the PSD model, not only to the GNSS, but also to the VLBI data.

## 4 ITRF2020 global combination

Past ITRF combinations up to ITRF2014 were operated using a two-step procedure (Altamimi et al. 2007, 2011, 2016): (1) stacking the individual time series to estimate a long-term solution per technique comprising station positions at a reference epoch, station velocities and daily EOPs; and (2) combining the resulting long-term solutions of the four techniques together with the local ties at colocation sites.

For the ITRF2020 computation, we opted for a one-step procedure, that is, to accumulate the full time series of the four techniques all together. Seven transformation parameters are estimated for each single daily, weekly, or VLBI

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 3** Histogram of the WRMS of the residuals from the trajectory models adjusted to the IGS repro3 station position time series without (in black) and with (in red) PSD</figcaption></figure>

session-wise solution with respect to the ITRF2020 global combined solution.

This strategy has multiple advantages, including the ease of imposing co-motion constraints between points of different techniques at colocation sites, but also the ease of transferring the instantaneous Center of Mass information as sensed by SLR to the three other technique networks, using the internal constraints approach. However, we also operated the two-step procedure and verified that we obtain the same results with both approaches, as it should mathematically be the case, provided that the same weighting is used for both space geodesy solutions and local ties. Indeed, the weighting is an important aspect in the ITRF combination analysis, as we deal with observations and constraints: global space geodesy solutions, local ties and co-motion constraints at colocation sites (cf. Eqs. 1 and 2). In addition, there
<figure><img src="image_9.png" /><figcaption>**Fig. 4** Trajectory of Tsukuba (Japan) site. Left GNSS and right VLBI. In blue raw data, in green the piecewise linear trajectories from ITRF2020 coordinates, and in red the trajectories obtained when adding the parametric PSD model</figcaption></figure>

are significant tie, velocity and seasonal signal discrepancies between technique solutions at a number of colocation sites that necessitate an iterative combination and empirical weighting process.

For the global ITRF2020 combination and to determine the relative weighting of the individual space geodetic technique solutions and the local ties, as well as the uncertainties of the velocity and seasonal signal constraints ($\sigma_v$ and $\sigma_f$ of Eqs. 1 and 2), we operated the following consecutive steps:

1. The weighting of the space geodesy solutions was determined by re-scaling the individual variance-covariance matrices by the following two factors:

* The a posteriori variance factors obtained from the individual stacking of each technique station position time series, whose square roots are: 2.92, 2.60, 3.93 and 11.56 for DORIS, GNSS, SLR and VLBI, respectively;

* Corrective factors aiming at accounting for the presence of time-correlated noise in the station position time series. While geodetic station position time series are known to exhibit time-correlated noise (Abban-danza et al. 2015; Williams et al. 2004; Williams and Willis 2006), the CATREF stacking model indeed assumes variable white noise only, and the above

variance factors are consequently under-estimated. To obtain corrective factors, we adjusted the ITRF2020 kinematic model of every station to its position time series together with either a variable white noise (only) model, or a variable white noise + power-law noise model. The parameters of the noise models (variances and spectral index) were estimated by restricted maximum likelihood, an unbiased alternative to classical maximum likelihood (Gobron et al. 2022; Harville 1977). The ratio between the formal velocity uncertainty obtained with the white + power-law noise model and the formal velocity uncertainty obtained with the white noise only model was then computed for each station and component. The corrective factors we applied in the ITRF2020 computation are the per-technique medians of the obtained ratios: 7.75 for GNSS, 2.12 for SLR and 2.69 for VLBI. For DORIS, although the median ratio is 2.00, we found empirically that it was necessary to use a corrective factor of 5.00. We noticed that selecting a corrective factor smaller than 5 would lead to an increase of normalized residuals exceeding the adopted threshold of 3 for a significant number of DORIS colocation sites.

2. For the purpose of equating station velocities at colocation sites, an iterative combination process of the four
technique velocity fields was operated as detailed in Sub-
sect. 4.1.

3. Using the weighting determined in steps 1 and 2, an iterative combination process of the four long-term solutions, adding local tie SINEX files at colocation sites, was performed. During this step, the uncertainties of local ties were re-evaluated and discrepant ties were down-weighted as described in Subsect. 4.2. Once the combination was stabilized, i.e., no normalized residual was larger than 3 in any component of both space geodesy estimates and local ties, the EOPs were then added to the combination. At this stage, the combination was iterated as necessary, and at each iteration EOP estimates were rejected if their normalized residuals were larger than 3. We also faced the situation where, when adding EOPs to the combination, small distortions of the SLR and VLBI frames appeared, namely an increase of position residuals for some stations, up to 4 mm, with normalized residuals exceeding the threshold of 3. We suspect this behavior to be caused by some mismatch between station position and polar motion estimates. In order to get the best of SLR and VLBI, and avoid these small distortions in the two corresponding frames, we down-weighted the ILRS and IVS polar motion estimates by a factor of 2. Doing so actually stabilized the two frames within the combination, and the small distortions were mitigated.

4. The global and final ITRF2020 combination was operated by accumulating the full time series of the four techniques all together, adopting the weighting for space geodesy time series and local ties as determined in the previous steps, and adding co-motion constraints where the $\sigma_v$ and $\sigma_f$ for the velocity and seasonal signals constraints, respectively, were chosen as described in Subsect. 4.1.

**4.1 Co-motion constraints at colocation sites**

**4.1.1 Equating station velocities at colocation sites**

In order to evaluate the level of the velocity agreement between techniques at colocation sites, and consequently adopt the appropriate $\sigma_v$ for the velocity constraint (Eq. 1) between colocated reference points, we conducted an iterative combination process of the four technique velocity fields embedded in the long-term solutions resulting from the individual stacking of the corresponding time series. In this analysis, we have chosen to constrain the couples of GNSS and reference points of the colocated instruments of the three other techniques that are available at these sites. For sites with multiple points (or segments in case of discontinuities) of the same technique, we constrained one couple of GNSS and other technique points, given the fact that the velocities of these multiple points (with some exceptions) are constrained technique per technique. Moreover, we

have chosen the longest and most stable time series of these constrained points by inspecting their respective detrended residuals. The initial $\sigma_v$ was chosen to be 0.1 mm/yr per component between GNSS/SLR and GNSS/VLBI colocated points, while 1 mm/yr was selected for the GNSS/DORIS pairs. Using $\sigma_v$ smaller than 1 mm/yr for GNSS/DORIS pairs leads to an increase of the normalized residuals exceeding the adopted threshold of 3 for a significant number of DORIS colocation sites. At each iteration of this combination of the velocity fields, the $\sigma_v$ of Eq. 1 was increased for sites where the velocity normalized residual exceeded the threshold of 3 in any of the three components.

**4.1.2 Equating station seasonal signals at colocation sites**

In the ITRF2020 global stacking of the four technique time series, we equated the seasonal signals (annual and semi-annual) of the same chosen couples of GNSS and reference points of the colocated instruments that were constrained to have the same velocity. The initial $\sigma_f$ of Eq. 2 was set to 0.1 mm per coefficient of the cosine and sine functions. In case of seasonal signal discrepancies between techniques at colocation sites larger than 1 mm in amplitude, $\sigma_f$ was increased to the level of their agreements.

In order to evaluate the benefit of estimating the seasonal signals, we operated an ITRF2020-like global combination without estimating the seasonal signals. We found that the horizontal velocity differences are almost all negligible for all stations and are less than 0.1 mm/yr, while the largest vertical velocity differences range between 0.5 and 1 mm/yr for about 20 poorly determined stations, with short time span of observations. The reduction of the velocity formal errors of about 10% when estimating the seasonal signals could be regarded as an improvement of the velocity determination.

**4.2 Handling of local ties**

In the same way as we use space geodesy solutions, local ties at colocation sites are used as observations, with full variance-covariance information provided in SINEX format. Note that the local ties are introduced in the combination with their corresponding survey epochs, as provided in their respective SINEX files. For the ITRF2020 computation, the variances of the used local ties were re-evaluated, in such a way that the uncertainty in any of the three components (East, North and Vertical) should not be below 2 mm (versus 3 mm in ITRF2014) using two scaling factors: one for the horizontal and one for the vertical components. The local tie lower bound uncertainty of 2 mm was adopted after a number of test combinations where we varied that threshold between 1 and 4 mm and found more coherent results with better agreement between space geodesy estimates and terrestrial ties, compared with ITRF2014 results as discussed at the end of
this subsection. The reason why we adopted two factors to re-scale the local tie variances is that the local tie vertical component is usually more precise than the horizontal ones, as it is determined using spirit leveling methods, independently from the horizontal components (Poyard et al. 2017).

Considering the number of station position discontinuities of the four techniques at colocation sites, but especially in case of GNSS time series, we had to adopt a comprehensive approach in selecting the segments to be tied together. For that purpose and as we did for the ITRF2014, we assigned the tie vectors to the longest and most stable segments, judging from their time series of residuals free from the linear trend and seasonal signals.

During the combination process of the long-term solutions of the four techniques, as described in Step 3 of Sect. 4, we iteratively down-weighted discrepant ties whose normalized residuals were exceeding the threshold of 3 in any of the three components. As in previous ITRF solutions, we chose to down-weight the local ties rather than space geodesy estimates in order to mainly maximize the consistency between individual technique solutions and ITRF2020.

Table 2 summarizes the main statistics of the local tie discrepancies in the ITRF2020 combination, in comparison with ITRF2014 results. Note that “Tie discrepancy” means difference between terrestrial ties and space geodesy estimates. As it can be seen in that table, ITRF2020 includes more tie vectors (from GNSS to other techniques reference points) and performs better than ITRF2014, as the percentage of tie discrepancies less than 5 mm is higher for all three techniques tied to GNSS. The reader can access the full list of tie discrepancies available at the ITRF2020 website.

## 4.3 ITRF2020 frame and seasonal signals definition

In the global ITRF2020 combination, the frame and seasonal signals defining parameters were specified as follows, noting that we used the concept of internal constraints for the origin and the scale, and the concept of minimum constraints for the orientation.

### 4.3.1 Origin

The origin of the ITRF2020 long-term frame is defined in such a way that there are zero translation parameters at epoch 2015.0 and zero translation rates between the ITRF2020 and the ILRS SLR long-term frame over the 1993.0–2021.0 time span, using the concept of internal constraints.

The origin of the ITRF2020 seasonal signals expressed in the CM (as sensed by SLR) frame is defined in such a way that there is no seasonal translation between the ITRF2020 seasonal signals and the input SLR solutions over the 1993.0–2021.0 time span, using the concept of internal constraints.

<figure><img src="image_10.png" /><figcaption>Fig. 5 Location of the reference frame sites used in the estimation of the 14 transformation parameters between ITRF2020 and ITRF2014 and their orientation alignment</figcaption></figure>

### 4.3.2 Scale

The scale of the ITRF2020 long-term frame is determined using the revised concept of internal constraints (see Appendix B for more details), in such a way that there are zero scale and scale rate between ITRF2020 and the scale and scale rate averages of VLBI selected sessions up to 2013.75 and SLR weekly solutions covering the 1997.75–2021.0 time span.

The scale of the ITRF2020 seasonal signals is determined using internal constraints, in such a way that the average of the seasonal scale variations between the selected SLR solutions and ITRF2020 is zero.

### 4.3.3 Orientation

The orientation of the ITRF2020 long-term frame is defined in such a way that there are zero rotation parameters at epoch 2015.0 and zero rotation rates between the ITRF2020 and ITRF2014. These two conditions were applied over a set of 131 reference stations located at 105 sites as illustrated by Fig. 5.

The orientation of the ITRF2020 seasonal signals is defined in such a way that there is no net seasonal rotation of that same core network.

## 5 ITRF2020 results

### 5.1 Station position kinematic model

The ITRF2020 is an augmented terrestrial reference frame comprising in addition to the classical parameters (station positions at a reference epoch $t_0 = 2015.0$ and station velocities), parametric models describing the trajectory of stations subject to post-seismic deformation, as well as annual and
<table><caption>**Table 2** Statistics of local tie discrepancies</caption><thead><tr><th rowspan="2">GNSS to</th><th colspan="2"># of tie vectors</th><th colspan="2"># of tie discrepancies &lt; 5 mm</th><th colspan="2">% of tie discrepancies &lt; 5 mm</th></tr><tr><th>ITRF2014</th><th>ITRF2020</th><th>ITRF2014</th><th>ITRF2020</th><th>ITRF2014</th><th>ITRF2020</th></tr></thead><tbody><tr><td>VLBI</td><td>60</td><td>77</td><td>27</td><td>38</td><td>45</td><td>50</td></tr><tr><td>SLR</td><td>49</td><td>53</td><td>14</td><td>19</td><td>29</td><td>36</td></tr><tr><td>DORIS</td><td>103</td><td>123</td><td>23</td><td>39</td><td>22</td><td>32</td></tr></tbody></table>

semi-annual periodic terms expressed in the CM (as sensed by SLR) frame, and in the Center of Figure (CF) frame.

Depending on the user need, the coordinates of a given station at epoch *t* can be computed using Eq. 3.

$$X(t) = X^i(t_0) + \dot{X}^i(t - t_0) + \delta X_{PSD}(t) + \delta X_f^i(t) \\ \forall t \in [t_i, t_{i+1}] \qquad (3)$$

where the couple $(X^i(t_0), \dot{X}^i)$ are the station position and velocity for segment (SOLN) *i* bounded by dates $(t_i, t_{i+1})$, $\delta X_{PSD}(t)$ is the total sum of the PSD corrections—the PSD corrections of successive earthquakes are cumulative, and $\delta X_f^i(t)$ is the total sum of the seasonal signal corrections for segment (SOLN) *i*. The following two subsections provide more details on how to compute $\delta X_{PSD}(t)$ and $\delta X_f^i(t)$.

## 5.2 Post-seismic deformation models

The full equations of the parametric models and their variance propagation used in the ITRF2020 generation are provided in Appendix D of this paper. While the ITRF2020 solution provides the usual/classical estimates (station positions at epoch 2015.0, station velocities and EOPs), the PSD models are also part of the ITRF2020 products. The users should be aware that they must compute the total sum of PSD corrections, $\delta X_{PSD}(t)$, using the equations supplied in Appendix D, to be then added to the ITRF2020 coordinates using Eq. 3 if the epoch of the needed position occurs during the post-seismic relaxation period. Failing to do so would introduce position errors at the decimeter level for many stations impacted by PSD. Full information regarding the PSD models fitted for stations subject to major earthquakes, and some useful subroutines in Fortran are available through the ITRF2020 website: https://itrf.ign.fr/en/solutions/ITRF2020.

## 5.3 ITRF2020 seasonal signals and geocenter motion

The seasonal (annual and semi-annual) signals were estimated for stations with time spans longer than 2 years, together with station positions, velocities and EOPs in the global ITRF2020 combination and expressed in the CM frame of SLR.

Figure 6 illustrates the ITRF2020 station annual signals, following East, North and Up components, respectively. They are expressed in the SLR-based CM frame where, as expected, geographically common mode patterns are evidenced, primarily due to the geocenter motion. In addition, it is well known that station seasonal displacements are spatially correlated (Collilieux et al. 2007; Dong et al. 2002).

The main motivation of estimating the seasonal signals consistent with the ITRF2020 whose long-term origin is inherited from SLR data is to provide to users, e.g., those dealing with satellite precise orbit determination, a means to refer the ITRF2020 coordinates to the quasi-instantaneous Earth CM, a point around which an artificial satellite is naturally orbiting.

To use this new feature, user should compute the seasonal coordinate variations $\delta X_f(t)$ using Eq. 4 and add it to the station coordinates at epoch *t*, using Eq. 3.

$$\delta X_f(t) = \sum_{i=1}^{2} \begin{pmatrix} a_x^i \\ a_y^i \\ a_z^i \end{pmatrix} \cos(2i\pi \cdot t) + \begin{pmatrix} b_x^i \\ b_y^i \\ b_z^i \end{pmatrix} \sin(2i\pi \cdot t) \quad (4)$$

where *t* is the decimal year and $(a_x^i, a_y^i, a_z^i, b_x^i, b_y^i, b_z^i)$ are the cosine and sine coefficients at frequency *i* (annual: *i* = 1; semi-annual: *i* = 2) provided in ITRF2020 for the given station.

Note that the seasonal motion of an ITRF2020 station provided in the CM frame is the sum of a net translation term that affects all stations equally and the station individual seasonal motion with respect to the center of figure of the solid Earth (CF). Recall that, by definition, the integral of seasonal deformation over the Earth's surface with respect to CF is zero (Blewitt 2003) and the net translation term is the opposite of the geocenter motion, namely the motion of CF with respect to CM.

Based on these concepts, it follows that from the ITRF2020 CM station seasonal signals, one can derive both seasonal geocenter motion and the CF station seasonal signals. To do so, we started by selecting a subset of ITRF2020 CM seasonal station motion coefficients. The selected subset excludes stations with poorly determined and abnormal seasonal coefficients and retains one single point per site with multiple stations that were constrained to have the same seasonal coefficients: the one for which the trace of the formal
<figure><img src="image_11.png" /><figcaption>**Fig. 6** ITRF2020 station annual signals following the three components: East (top), North (middle) and Up (bottom). The amplitudes (A) and phases (φ) are represented by the arrow lengths and orientations (counted from the north), respectively, following the model A. cos(2π.t - φ), with t in decimal year</figcaption></figure>

covariance matrix of the seasonal coefficients is minimal. The final selection contains a cleaned ensemble of 1124 stations. We then computed a weighted sum of the selected CM station seasonal signals aiming at approximating the integral of the Earth's surface seasonal deformation with respect to CM. The adopted weights are the areas of the stations' spherical Voronoï cells (Caroli et al. 2010), spatially smoothed so that nearby stations have similar weights. From the arguments above, it follows that this weighted sum approximates the seasonal motion of CF with respect to CM, hereafter noted ΔX<sub>CF/CM</sub>.

<table><caption>**Table 3** Seasonal displacement of CF with respect to CM = opposite of the seasonal geocenter motion</caption><thead><tr><th></th><th>Amplitude<br>mm</th><th>Phase<br>degrees</th></tr></thead><tbody><tr><td>Annual x</td><td>1.23</td><td>-123.2</td></tr><tr><td>±</td><td>0.16</td><td>7.2</td></tr><tr><td>Annual y</td><td>3.48</td><td>152.9</td></tr><tr><td>±</td><td>0.15</td><td>2.5</td></tr><tr><td>Annual z</td><td>2.76</td><td>-139.5</td></tr><tr><td>±</td><td>0.33</td><td>6.8</td></tr><tr><td>Semi-annual x</td><td>0.49</td><td>107.2</td></tr><tr><td>±</td><td>0.15</td><td>18.1</td></tr><tr><td>Semi-annual y</td><td>0.22</td><td>1.6</td></tr><tr><td>±</td><td>0.15</td><td>39.0</td></tr><tr><td>Semi-annual z</td><td>1.19</td><td>30.5</td></tr><tr><td>±</td><td>0.33</td><td>15.5</td></tr></tbody></table>

The obtained ITRF2020 seasonal coefficients of vector ΔX<sub>CF/CM</sub> are given in Table 3. Subtracting this model from the ITRF2020 CM station seasonal signals yields the ITRF2020 CF station seasonal signals, which are also provided at the ITRF2020 website. More details on the derivation of the ITRF2020 geocenter motion model will be given in a follow-up paper.

Following Table 3, the seasonal motion of CF with respect to CM, along any x, y, z axis, say Δx<sub>CF/CM</sub> along x, is obtained by Eq. 5.

$$ \Delta x_{CF/CM}(t) = A_1 \cos \left( 2\pi t - \frac{\pi}{180} \cdot \phi_1 \right) + A_2 \cos \left( 4\pi t - \frac{\pi}{180} \cdot \phi_2 \right) \quad (5) $$

where *t* is the time in decimal years, A₁ and A₂ are the annual and semi-annual amplitudes, φ₁ and φ₂ are the annual and semi-annual phases expressed in degrees, respectively, as listed in Table 3. This equation is valid for each one of the three axes x, y, z.

The reader should note that by convention, the IERS Conventions (Petit and Luzum 2010, Ch. 4) defines seasonal geocenter motion as the motion of the Earth's CM with respect to CF. To obtain seasonal geocenter motion as defined in the IERS Conventions, the user can either add 180 degrees to the phase values listed in Table 3, or change sign of Eq. (5).

## 5.4 Usage of the ITRF2020 seasonal signals

When using the ITRF2020 station position kinematic model (3), the user should decide, depending on his application, if the seasonal signals, i.e., adding δX<sub>f</sub>(t), should be considered, and if so, whether the seasonal signals referred to CM
or CF should be used. In this regard, it should be noted that,
in principle:

* The origin of the *linear* ITRF2020 coordinates reflects CM on secular time scale, but CF on shorter (inter-annual, seasonal, sub-seasonal, etc.) time scales, as a consequence of the linear kinematic model used to represent station motions. Linear coordinates can indeed not account for net nonlinear translations of the Earth's surface, and CF is by definition the point with respect to which there is no net translation of the Earth's surface (see full development of this argument in (Dong et al. 2003));

* If the ITRF2020 seasonal signals referred to CF are added, the origin of the resulting coordinates remains the same as in the purely linear case: it still reflects CM on the long term, but CF at both the annual and semi-annual periods (by construction of the ITRF2020 CF seasonal signals) and all other frequencies (since they are not represented in the kinematic model). However, compared to the purely linear case, the resulting coordinates additionally account for the individual seasonal station motions with respect to CF;

* On the other hand, if the ITRF2020 seasonal signals referred to CM are added, the origin of the resulting coordinates changes compared to the purely linear case: it now reflects CM on both the long term and at the annual and semi-annual periods, but still reflects CF on all other time scales (e.g., inter-annual, sub-seasonal).

In practice, different usage scenarios could be the following:

* *Precise orbit determination:* In this application, the reference station positions should describe as well as possible the instantaneous shape of the Earth, and be referred as closely as possible to the instantaneous CM. Hence, using as reference positions the linear ITRF2020 coordinates augmented with the seasonal signals referred to CM, rather than the linear ITRF2020 coordinates only, can be expected to be beneficial. This would indeed make the reference positions fit the Earth's shape and follow CM not only on the long term, but also at the annual and semi-annual periods.

* *Alignment of global solutions:* When aligning series of global instantaneous solutions derived from space geodetic data to a long-term linear reference frame by means of Helmert transformations, part of the nonlinear station motions aliases into the estimated Helmert parameters (Collilieux et al. 2012). The incorporation of annual and semi-annual signals in the target reference frame does however resolve this aliasing at the annual and semi-annual frequencies. Hence, aligning global instantaneous

solutions to the linear ITRF2020 frame augmented with seasonal signals (referred to either CM or CF), rather than to the linear ITRF2020 frame only, can be expected to be beneficial and result in seasonal station motions more accurately retained in the aligned solutions. This holds equally true whether the ITRF2020 seasonal signals referred to CM or CF are used. Employing the ITRF2020 seasonal signals referred to one or the other origin would only change the origin of seasonal station motions in the aligned solutions. Namely, if the ITRF2020 seasonal signals referred to CM are used, the seasonal station motions in the aligned solutions will also be referred to CM, and will include a seasonal translational motion common to all stations (i.e., seasonal geocenter motion). On the other hand, if the ITRF2020 seasonal signals referred to CF are used, the seasonal station motions in the aligned solutions will also be referred to CF, and will be free of this global seasonal translational common mode. Whether the ITRF2020 seasonal signals referred to CM or CF should be used does thus depend on the usage intended for the aligned solutions. For instance, if the station positions in the aligned solutions are intended to serve as reference station positions for precise orbit determination, then the ITRF2020 seasonal signals referred to CM should be preferred (see also previous item). On the other hand, if the aligned solutions are mainly intended to study the Earth's surface deformation, it seems more natural to have them free of a global seasonal common mode, hence to use the ITRF2020 seasonal signals referred to CF. Note finally that if an instantaneous solution was aligned to "linear ITRF2020 + CM seasonal signals", it is straightforward to align it to "linear ITRF2020 + CF seasonal signals" instead: for that purpose, the offset $\Delta x_{CF/CM}(t)$ given by Eq. (5) should simply be removed from all station positions. The inverse transformation is achieved by adding this offset instead of removing it.

* *Alignment of local or regional GNSS solutions:* In this case, since seasonal displacements are spatially correlated (Collilieux et al. 2007; Dong et al. 2002), an alignment to the linear ITRF2020 frame augmented with seasonal signals (referred to either CM or CF) would introduce a seasonal common mode in the resulting station position time series. If the user wishes to avoid such a regional, seasonal common mode, he should use as a reference the linear ITRF2020 frame without adding seasonal signals.

It is also worth mentioning here that non-tidal loading corrections were not applied to the ITRF2020 input data, nor during the ITRF2020 computation, so that part of the ITRF2020 seasonal signals reflects seasonal loading deformation. As a consequence, station displacements predicted by a non-tidal loading model should not be added to the lin-
<table><caption>**Table 4** Translation parameters from ITRF2020 to past ITRF solutions at epoch 2015.0 and their rates in mm and mm/yr</caption><thead><tr><th>Solution</th><th>TX</th><th>TY</th><th>TZ</th></tr></thead><tbody><tr><td>ITRF2014</td><td>-1.4</td><td>-0.9</td><td>1.4</td></tr><tr><td>Rates</td><td>0.0</td><td>-0.1</td><td>0.2</td></tr><tr><td>ITRF2008</td><td>0.2</td><td>1.0</td><td>3.3</td></tr><tr><td>Rates</td><td>0.0</td><td>-0.1</td><td>0.1</td></tr><tr><td>ITRF2005</td><td>2.7</td><td>0.1</td><td>-1.4</td></tr><tr><td>Rates</td><td>0.3</td><td>-0.1</td><td>0.1</td></tr></tbody></table>

ear ITRF2020 coordinates augmented with seasonal signals, unless mean annual and semi-annual variations are removed from the predicted loading displacements beforehand. Otherwise, seasonal loading deformation would be accounted twice.

## 5.5 ITRF2020 long-term origin

As stated earlier in this article, the long-term origin of the ITRF2020 is defined by the ILRS SLR input data only, and satisfies the conditions of no translation and no translation rate between ITRF2020 and the SLR long-term frame. Figure 7 illustrates the full time series of the ILRS SLR origin components with respect to the ITRF2020, as result of the global stacking. From that figure, one can see no discernable offset or drift, nor seasonal signals as these were constrained to zero (over the 1993.0–2021.0 time span), using the internal constraints approach, in order to express the ITRF2020 station seasonal signals in the CM frame as sensed by SLR.

One way to evaluate the accuracy and performance of the ITRF2020 origin is to estimate the translation components with respect to past ITRF solutions, namely ITRF2014, ITRF2008 and ITRF2005 whose origins were also defined by SLR data submitted in the form of time series. Table 4 lists the translation components and their rates from ITRF2020 to the past three frames. In that table, we can see variations between -1.4 and 3.3 mm for the translations at epoch 2015.0, and between -0.1 and 0.3 mm/yr for the translation rates. In other words, we can stipulate that the origins of the successive ITRF solutions stand roughly within a sphere of 5 mm in diameter and 0.5 mm/yr for their time variation. These figures are an indication of the level of the ITRF origin intrinsic accuracy achievable today using SLR data, dominated by LAGEOS I and II observations (Pavlis and Luceri 2022) and are consistent with previous independent analyses (Collilieux et al. 2014).

Figure 8 displays the time series of the DORIS geocenter components with respect to ITRF2020, as obtained from the ITRF2020 global stacking of all technique time series. The estimated DORIS long-term origin component offsets with respect to SLR are small: -2.7, 1.1 and 5.4 mm, at epoch 2015.0, for X, Y and Z, respectively, with almost zero drifts.

However, we can clearly see in Fig. 8 significant seasonal variations in X, but more pronounced in Y, suggesting that SLR (which defines the ITRF2020 origin on the seasonal term) and DORIS seasonal geocenter motions are inconsistent for these two components. The DORIS Z-geocenter component is very scattered and still suffering from systematic errors, such as solar radiation pressure mis-modeling (Gobinddass et al. 2009).

## 5.6 ITRF2020 long-term scale

This is the first time of the ITRF history where we are in presence of independent and competitive scales stemming from the 4 techniques (VLBI, SLR, GNSS and DORIS) that contributed input data for the ITRF2020. Note that as stated earlier, the IGS/GNSS scale is based on the calibrated z-Phase Center Offsets (z-PCOs) of the Galileo satellites published by GSC (GSC 2022), consistently to which the GPS and GLONASS satellite z-PCOs have been re-evaluated (Rebischung 2020). This is also the first time where the ILRS did a great effort in handling and estimating the station-dependent range biases (Luceri et al. 2019) that greatly impact the scale of the SLR frame by up to 1 ppb (about 6 mm at the equator), as demonstrated by (Appleby et al. 2016).

The left panel of Fig. 9 illustrates in orange and red the full time series of VLBI session-wise scales and in light and dark blue the fortnightly/weekly SLR scales with respect to ITRF2020. As it can be seen in that figure, we detected an unexpected VLBI scale drift after around epoch 2014.0, but also an offset in the SLR scale before epoch 1997.75, as well as poor and scattered scale determinations before 1993.0 when only LAGEOS I satellite data were available. The reasons for these VLBI scale drift and SLR scale offset are not clear to the authors and are under investigation by the IVS and ILRS. In order to use the best and most stable scale time series of both techniques for the ITRF2020 scale definition, we retained the ILRS weekly solutions from 1997.75 on (dark blue dots at Fig. 9), and considered only the VLBI sessions before 2013.75 for which station networks have convex hulls with volumes larger than $10^{19}$ m$^3$. To empirically localize the epoch of the VLBI scale drift, we operated 8 test combinations involving the 4 technique long-term solutions where we varied the starting epoch of the VLBI scale drift between 2012.0 and 2021.0 and ended up by finding that assigning the start of the VLBI scale drift at epoch 2013.75 gives the best scale agreement between SLR and VLBI, namely a 0.15 ppb offset (equivalent to 1 mm at the equator) at epoch 2015.0 and zero drift. This is the first time of the ITRF history where the scale agreement between the two techniques reaches that level, compared with the ITRF2014 results (Altamimi et al. 2016) where the agreement was at the level of 1.39 ppb (about 8 mm at the equator).
<figure><img src="image_12.png" /><figcaption>**Fig. 8** Time series of DORIS geocenter components with respect to ITRF2020, in mm</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 9** SLR and VLBI scale time series with respect to ITRF2020 (left), adding DORIS and GNSS (right), in mm. Light blue and orange: all SLR and VLBI solutions, dark blue and red: selected solutions whose average defines the ITRF2020 scale, black: DORIS and green: GNSS</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Table 5** Scale offsets and rates with respect to ITRF2020</figcaption></figure>

<table><thead><tr><th>Solution</th><th>Scale at 2015.0<br>ppb</th><th>Scale rate<br>ppb/yr</th></tr></thead><tbody><tr><td>IDS/DORIS</td><td>1.386</td><td>0.028</td></tr><tr><td>±</td><td>0.037</td><td>0.003</td></tr><tr><td>IGS/GNSS</td><td>0.682</td><td>0.018</td></tr><tr><td>±</td><td>0.018</td><td>0.001</td></tr><tr><td>ILRS/SLR</td><td>-0.075</td><td>0.000</td></tr><tr><td>±</td><td>0.038</td><td>0.004</td></tr><tr><td>IVS/VLBI</td><td>0.075</td><td>0.000</td></tr><tr><td>±</td><td>0.040</td><td>0.003</td></tr></tbody></table>

Note that using the full scale time series of SLR and VLBI would introduce a scale discrepancy between these two techniques of 0.3 ppb at epoch 2015.0 and a significant rate of 0.03 ppb/yr.

The right panel of Fig. 9 displays the scale time series of IGS/GNSS in green and IDS/DORIS in black, respectively, superimposed over the SLR and VLBI time series, where one can see their behaviors, offsets and rates with respect to ITRF2020. The estimated scale offsets at epoch 2015.0 and rates with respect to ITRF2020 of the four technique solutions are listed in Table 5. The DORIS and GNSS scale time series were not used in the ITRF2020 scale definition, given their significant offsets and rates (see the values listed in Table 5 and Fig. 9) compared with the high level agreement between SLR and VLBI.

<figure><img src="image_15.png" /><figcaption>**Fig. 10** ITRF2020 horizontal site velocities with formal error less than 1 mm/yr. Major plate boundaries are shown according to (Bird 2003)</figcaption></figure>

## 5.7 ITRF2020 velocities

The ITRF2020 site velocities with formal errors less than 1 mm/yr are illustrated by Fig. 10 for the horizontal and by Fig. 11 for the vertical components, respectively. Compared with the ITRF2014 horizontal velocity field, ITRF2020 contains more sites in some large tectonic plates, such as North and South Americas, Australia, and to some extent Antarc- tica. The angular velocities of these plates are expected to be improved in the ITRF2020 Plate Motion Model which is under preparation at the time of writing this article, compared to the one published in (Altamimi et al. 2017).
<figure><img src="image_16.png" /><figcaption>Fig. 11 ITRF2020 vertical site velocities with formal error less than 1 mm/yr</figcaption></figure>

Inspecting Fig. 11, one can easily distinguish regional patterns in the ITRF2020 vertical velocities, especially in North America, Greenland, Fennoscandia and Antarctica, reflecting the uplift caused by not only Glacial Isostatic Adjustment, but also recent ice melting, in particular in Greenland (Métivier et al. 2020). Thanks to the density of ITRF2020 sites in North America, the vertical velocity pattern clearly shows the significant subsidence belt, right south of the Canadian uplift area, consistent with the results discussed in (Kreemer et al. 2018).

In order to assess the overall velocity agreement between techniques at colocation sites, Fig. 12 shows the histograms of the velocity differences between GNSS and the three other techniques, following the three components (East, North and Up). As a statistical estimator, we computed the standard deviation (SD) of the different distributions illustrated in Fig. 12. For the horizontal components, we found that the SD is about 0.7 mm/yr for VLBI and SLR and 1.2 mm/yr for DORIS, respectively. The SD values for the vertical components range between 1.3 mm/yr and 1.6 mm/yr for the three techniques (VLBI, SLR, DORIS), compared with GNSS.

## 5.8 ITRF2020 earth orientation parameters

The ITRF2020 global stacking of the four technique time series rigorously combines station positions, velocities and seasonal signals, together with EOPs: Polar Motion (PM) components and their daily rates (PM-rate), as well as Universal Time (UT1-UTC) and Length of Day (LOD). The EOPs are thus well integrated and rigorously referenced in the ITRF2020 frame. ITRF2020 PM components are the results of the combination of the four technique time series, PM-rate values are determined by combining GNSS and VLBI results, while UT1-UTC and LOD are determined using VLBI data uniquely, in order to avoid possible biased estimates from satellite techniques (Ray 1996, 2009). The full consistency between the ITRF2020 frame and EOPs is ensured via the three rotation parameters appearing in the first three lines

of Eq. 8 of Appendix A. In particular, PM coordinates play a significant “colocation” role in linking the four technique frames, via the rotation parameters around the X and Y axes.

The level of agreement between the four techniques in Polar Motion estimates is illustrated by the post-fit residuals in x and y components shown in Fig. 13, where one can see that the GNSS PM series dominates the other technique estimates. The WRMS values computed over the post-fit residuals between the combined and the individual polar motion time series are (for the couple x and y components): (20, 20), (196, 204), (260, 251), (305, 294) in micro arc-seconds, for GNSS, VLBI, SLR and DORIS, respectively. The reader should note that these WRMS values are computed over the full time series of the four techniques which makes SLR and VLBI values larger than what is observed during the last two decades.

## 5.9 Transformation parameters between ITRF2020 and ITRF2014

For a number of geodetic applications, the ITRF users need the transformation parameters between the new frame and past ITRF versions. For that purpose, the same 131 stations displayed in Fig. 5 that were used to ensure the alignment of the ITRF2020 orientation and orientation rate to the ITRF2014, were also used to estimate the transformation parameters between the two frames. The main criteria for the selection of these 131 stations are (1) to have the best possible site distribution; (2) to involve as many as possible VLBI, SLR, GNSS and DORIS stations and (3) to have the best agreement between the two frames in terms of post-fit residuals of the 14-parameter transformation. Regarding this third criteria, the WRMS values of the 14-parameter similarity transformation fit are 1.6, 1.8 and 3.4 mm in position (at epoch 2015.0) and 0.2, 0.1, 0.3 mm/yr in velocity, in East, North and Vertical components, respectively. Table 6 lists the transformation parameters from ITRF2020 to ITRF2014, to be used with the transformation formula given by Eq. (6).

$$ \left\{ \begin{aligned} \begin{pmatrix} x \\ y \\ z \end{pmatrix}_{i14} &= \begin{pmatrix} x \\ y \\ z \end{pmatrix}_{i20} + T + D \begin{pmatrix} x \\ y \\ z \end{pmatrix}_{i20} + R \begin{pmatrix} x \\ y \\ z \end{pmatrix}_{i20} \\ \begin{pmatrix} \dot{x} \\ \dot{y} \\ \dot{z} \end{pmatrix}_{i14} &= \begin{pmatrix} \dot{x} \\ \dot{y} \\ \dot{z} \end{pmatrix}_{i20} + \dot{T} + \dot{D} \begin{pmatrix} x \\ y \\ z \end{pmatrix}_{i20} + \dot{R} \begin{pmatrix} x \\ y \\ z \end{pmatrix}_{i20} \end{aligned} \right. \quad (6) $$

where $i14$ designates ITRF2014 and $i20$ ITRF2020, $T$ is the translation vector, $T = (T_x, T_y, T_z)^T$, $D$ is the scale factor
<figure><img src="image_17.png" /><figcaption>Fig. 13 ITRF2020 post-fit residuals of polar motion in milli-arc-seconds</figcaption></figure>

<figure><img src="image_18.png" /></figure>

and R is the matrix containing the rotation angles, given by

$$R = \begin{pmatrix} 0 & -R_z & R_y \\ R_z & 0 & -R_x \\ -R_y & R_x & 0 \end{pmatrix}$$

The dotted parameters designate their time derivatives. The values of the 14 parameters are those listed in Table 6. Note that the inverse transformation from ITRF2014 to ITRF2020 follows by interchanging (i20) with (i14) and changing the sign of the transformation parameters.

## 6 Conclusion

For the first time of the ITRF history, the ITRF2020 was generated by rigorously stacking the full time series of the four techniques all together, adding local ties and co-motion constraints at colocation sites. In addition to station posi-

tions and velocities and EOPs, refined parametric functions were also determined to accurately describe nonlinear station motions for both seasonal signals and post-seismic deformation, yielding an augmented ITRF2020 frame.

Using the concept of internal and minimum constraints, the ITRF2020 frame and seasonal signal parameters are comprehensively defined, thus preserving the ITRF2020 network shape which approximates the shape of the constantly deforming Earth’s surface.

The ITRF2020 seasonal signals are provided in both the CM (as sensed by SLR) and the CF frames, together with a consistent new seasonal geocenter motion model, offering to the users several options to choose from, allowing to model seasonal deformations. We provided in Sect. 5.4 guidance to users to help choosing which options are the most suitable for their applications.

We verified that the ITRF2020 PSD parametric models determined using IGS/GNSS data fit also the station position time series of the three other techniques that are sharing the
<table><caption>**Table 6** Transformation parameters at epoch 2015.0 and their rates from ITRF2020 to ITRF2014, to be used with Eq. (6)</caption><thead><tr><th></th><th>Tx<br>mm</th><th>Ty<br>mm</th><th>Tz<br>mm</th><th>D<br>ppb</th><th>Rx<br>mas</th><th>Ry<br>mas</th><th>Rz<br>mas</th></tr><tr><th></th><th>T&#x0307;<br>mm/y</th><th>T&#x0307;<br>mm/y</th><th>T&#x0307;<br>mm/y</th><th>D&#x0307;<br>ppb/y</th><th>R&#x0307;<br>mas/y</th><th>R&#x0307;<br>mas/y</th><th>R&#x0307;<br>mas/y</th></tr></thead><tbody><tr><th>±</th><td>-1.4<br>0.2</td><td>-0.9<br>0.2</td><td>1.4<br>0.2</td><td>-0.42<br>0.03</td><td>0.00<br>0.07</td><td>0.00<br>0.06</td><td>0.00<br>0.07</td></tr><tr><th>±</th><td>0.0<br>0.2</td><td>-0.1<br>0.2</td><td>0.2<br>0.2</td><td>0.00<br>0.03</td><td>0.00<br>0.07</td><td>0.00<br>0.06</td><td>0.00<br>0.07</td></tr></tbody></table>

same colocation sites. This is a key indicator of the high performance of the PSD parametric models in describing the trajectory of ITRF2020 sites impacted by major earthquakes.

By estimating the origin translational components between ITRF2020 and past ITRF solutions determined via station position time series, namely ITRF2005, ITRF2008 and ITRF2014, we conclude that the current accuracy and stability of the ITRF long-term origin achievable today using SLR data is at the level of or better than 5 mm in position and 0.5 mm/yr in time variation.

By filtering the scale time series of both SLR and VLBI solutions, and therefore using the best fit of the two series, we ended with a scale agreement at the level of 0.15 ppb (≈ 1 mm at the equator) at epoch 2015.0, with no significant rate, a key finding of the ITRF2020 analysis.

Compared to past ITRF solutions, ITRF2020 results show modest, but noticeable improvements in terms of agreement between terrestrial ties and space geodesy estimates at colocation sites. We still believe that most (if not all) of the tie discrepancies are caused by technique systematic errors rather than local ties. Improving the consistency between the two ensembles would not be possible without improving the geodetic infrastructure by investing in new generation of SLR and VLBI technologies, with a better distribution between the northern and southern hemispheres, but also by finding a mechanism to reduce, if not to eliminate, the IGS/GNSS station position discontinuities. Finally, as we stated in the ITRF2014 article (Altamimi et al. 2016): *VLBI in particular needs to evolve toward more frequent global observational session schedules, with increased number and well-distributed stations.*

Despite the high level results on reference frame determination presented in this paper, thanks to improved input data by the 4 techniques and the combination strategy, we are still far from the stringent science requirement of 1 mm accuracy and 0.1 mm/yr stability of the ITRF defining parameters. A factor of at least 5 will be needed to reach that goal.

As a final and important note, we repeat the last sentence of the ITRF2014 article (Altamimi et al. 2016), that is (with slight modification): *Improving the geodetic infrastructure is the prerequisite for the long-term sustainability and accuracy of the ITRF, as recognized by the United Nations General*

*Assembly resolution (2015) on the global geodetic reference frame for sustainable development.*

**Acknowledgements** The International Terrestrial Reference Frame (ITRF) is the result of a global collaboration of hundreds of institutions around the worlds: from the build-up of geodetic observatories, satellite missions, data collection, analysis and combination, to the ITRF generation, thanks to the investment of national mapping agencies, space agencies and universities. We are indebted to all Analysis and Combination Centers of the IAG Technique Services of the four techniques who constantly provide data for ITRF solutions which would not exist without their valuable contributions. We acknowledge valuable comments and suggestions provided by Rick Bennett and two anonymous reviewers which improved the content of this article.

**Author contributions** ZA conceptualized and designed the initial CATREF combination model used in the ITRF computation and wrote the full software code, with the help of XC for its updates. XC and PR greatly contributed to the CATREF equation improvements and updates. PR fitted the PSD parametric models using IGS/GNSS data, designed geocenter motion model estimation, and greatly contributed to the ITRF2020 analysis and validation of the results. LM developed the method to detect co-seismic discontinuities in station position time series and participated in the time series analysis. ZA operated the final ITRF2020 global combination and wrote the first draft of the paper. All authors contributed to regular discussions throughout the ITRF2020 project, read and commented on the final paper and provided contributions to the interpretation of the results.

**Funding** The ITRF activities are funded by the Institut National de l'Information Géographique et Forestière (IGN), France, and partly by the Centre National d'études Spatiales (CNES), under TOSCA Grant.

**Data availability** The input data used in the ITRF2020 computation are available at the NASA Crustal Dynamics Data Information System (CDDIS), https://cddis.nasa.gov, or through the Analysis Center Coordinators of the IAG Technique Services: http://ids-doris.org, http://www.igs.org, http://ilrs.gsfc.nasa.gov, http://ivscc.gsfc.nasa.gov. The local ties used in ITRF2020 are available at the ITRF2020 web site.

**Code availability** CATREF Software is available free of charge to research groups within institutions committed to sign an agreement with IGN for research, non-commercial activities.

## Declarations

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

## A Combination model

The main two equations of the combination model implemented in CATREF Software are given below. They involve a 14-parameter similarity transformation, station positions and velocities, periodic signals and EOPs and are written as:

$$ 
\left\{ 
\begin{aligned}
X_s^i &= X_c^i + (t_s^i - t_0) \dot{X}_c^i \\
&\quad + T_k + D_k X_c^i + R_k X_c^i \\
&\quad + (t_s^i - t_k) \left[ \dot{T}_k + \dot{D}_k X_c^i + \dot{R}_k X_c^i \right] \\
&\quad + \sum_{j=1}^{n_f} a_j^i \cos[\omega_j(t_s^i - t_0)] + b_j^i \sin[\omega_j(t_s^i - t_0)] \\
\dot{X}_s^i &= \dot{X}_c^i + \dot{T}_k + \dot{D}_k X_c^i + \dot{R}_k X_c^i
\end{aligned}
\right. 
\quad (7) 
$$

$$ 
\left\{ 
\begin{aligned}
x_s^p &= x_c^p + R2_k \\
y_s^p &= y_c^p + R1_k \\
UT_s &= UT_c - \frac{1}{f} R3_k \\
\dot{x}_s^p &= \dot{x}_c^p \\
\dot{y}_s^p &= \dot{y}_c^p \\
LOD_s &= LOD_c
\end{aligned}
\right. 
\quad (8) 
$$

where for each point i, $X_s^i$ (at epoch $t_s^i$) and $\dot{X}_s^i$ are positions and velocities of technique solution s and $X_c^i$ (at epoch $t_0$) and $\dot{X}_c^i$ are those of the combined solution c. For each individual frame k, as implicitly defined by solution s, $D_k$ is the scale factor, $T_k$ the translation vector and $R_k$ the rotation matrix. The dotted parameters designate their derivatives with respect to time. The translation vector $T_k$ is composed of three origin components, namely $T_x, T_y, T_z$, and the rotation matrix of three small rotation parameters: $R_x, R_y, R_z$, following the three axes, respectively, X, Y, Z. $t_k$ is a conventionally selected epoch of the 7 transformation parameters. In addition to station positions and velocities, Eq. 7 involves also the periodic terms modeled as a sine wave function, adding the total sum of the contributions of all the frequencies considered, $n_f$ is the number of frequencies, $\omega_j = \frac{2\pi}{\tau_j}$, where $\tau_j$ is the period of the jth frequency, e.g., annual, semi-annual, etc.

In addition to the six station position and velocity parameters, each frequency adds six parameters per station, i.e., $(a_x^i, a_y^i, a_z^i, b_x^i, b_y^i, b_z^i)^T$.

EOPs are added by Eq. 8, making use of pole coordinates $x_s^p, y_s^p$ and universal time $UT_s$ as well as their daily rates $\dot{x}_s^p, \dot{y}_s^p$ and $LOD_s$, where $f = 1.002737909350795$ is the conversion factor from UT into sidereal time. The link between the combined frame and the EOPs is ensured via the three rotation parameters appearing in the first three lines of Eq. 8.

Note that Eq. 7 uses the linearized form of the general similarity transformation formula, neglecting second and higher order terms (Petit and Luzum 2010, Ch. 4), (Altamimi and Dermanis 2012).

The first line of Eq. 7 and the entire Eq. 8 are used to estimate long-term solutions, by accumulating (rigorously stacking) time series of station positions and EOPs. They are also the ones used, together with local ties and co-motion constraints, in the global ITRF2020 combinations. The entire two equations are used to combine long-term solutions, e.g., in the step 3 described in Sect. 4.

## B Concept of internal constraints revisited

We have introduced the concept of internal constraints at the occasion of the ITRF2005 (Altamimi et al. 2007) where we started to use as input data time series of station positions and EOPs. The fundamental property of the internal constraints is to preserve the intrinsic long-term frame physical parameters such as the origin and scale of SLR, and the scale of VLBI when stacking their corresponding time series, individually. We have extended the concept of internal constraints for the purpose of the ITRF2020 global combination, where the time series of several techniques are stacked all together, in order to ensure that the average of the selected SLR and VLBI scale time series is strictly zero. Denoting $P_i$ the time series of a transformation parameter P (such as the scale) of technique i, and n the number of techniques involved, the generalized internal constraints condition can be written as:

$$ \sum_{i=1}^{n} (A_i^T A_i)^{-1} A_i^T P_i = 0 \quad (9) $$

where $P_i = [P_{i,1}, \dots, P_{i,n_i}]^T$, $n_i$ is the dimension of the time series of the parameter P of technique i, and $A_i$ takes the form:

$$ A_i = \begin{bmatrix} 1 & (t_{i,1} - t_0) \\ \vdots & \vdots \\ 1 & (t_{i,n_i - t_0}) \end{bmatrix} $$
## C Periodic signals

The inclusion of the transformation parameters in the combination model offers a number of options for the reference frame definition, but also for the periodic signals, such as minimum or internal constraint approaches similar to those developed in (Altamimi et al. 2007).

### C.1 Periodic signals: minimum constraints

In the stacking case, the equations of the minimum constraints approach take the form

$$ (A^T A)^{-1} A^T \begin{pmatrix} \delta a_x^i \\ \delta a_y^i \\ \delta a_z^i \end{pmatrix} = (A^T A)^{-1} A^T \begin{pmatrix} a_{xR}^i - a_x^i \\ a_{yR}^i - a_y^i \\ a_{zR}^i - a_z^i \end{pmatrix} = 0 \quad (10) $$

$$ (A^T A)^{-1} A^T \begin{pmatrix} \delta b_x^i \\ \delta b_y^i \\ \delta b_z^i \end{pmatrix} = (A^T A)^{-1} A^T \begin{pmatrix} b_{xR}^i - b_x^i \\ b_{yR}^i - b_y^i \\ b_{zR}^i - b_z^i \end{pmatrix} = 0 \quad (11) $$

where A is the well known design matrix of partial derivatives of the 7 transformation parameters, and $(a_{xR}^i, a_{yR}^i, a_{zR}^i, b_{xR}^i, b_{yR}^i, b_{zR}^i)^T$ are the reference values for each frequency $i$ that could be taken from an external loading model, or as zeros over a set of well distributed stations.

### C.2 Periodic signals: internal constraints

In the stacking case, one possibility is to constrain to zero, internally, the periodic signals embedded in the time series of a given transformation parameter $[P(t_1), ..., P(t_K)]$. In such a case the internal constraint equation takes the form:

$$ (B^T B)^{-1} B^T [P(t_1), \dots, P(t_K)]^T = 0 \quad (12) $$

where K is the number of the individual daily or weekly solutions, and B is the matrix of partial derivatives given by:

$$ B = \begin{pmatrix} \cos(\omega_i t_1) & \sin(\omega_i t_1) \\ \vdots & \vdots \\ \cos(\omega_i t_K) & \sin(\omega_i t_K) \end{pmatrix} \quad (13) $$

## D Equations of post-seismic deformation models

After an earthquake, the position of a station during the post-seismic trajectory at an epoch $t$ is described by an additional term $\delta X_{PSD}(t)$ in Eq. 3.

ITRF2020 provides the PSD parameters in the local frame (Easting, Northing, Up) at epoch $t$. We remind that displacements in the local frame can be converted to displacements in the global frame with the following equation:

$$ \begin{bmatrix} \delta x \\ \delta y \\ \delta z \end{bmatrix} = R \cdot \begin{bmatrix} \delta E \\ \delta N \\ \delta U \end{bmatrix} \quad (14) $$

where R is the transformation (Jacobian) matrix from topocentric to geocentric frame.

For each component $L \in \{E,N,U\}$, the total PSD displacement in the local frame can be computed by:

$$ \delta L(t) = \sum_{i=1}^{n'} A_i^l \log \left( 1 + \frac{t - t_i^l}{\tau_i^l} \right) + \sum_{i=1}^{n'} A_i^e \left( 1 - e^{-\frac{t - t_i^e}{\tau_i^e}} \right) \quad (15) $$

where:

$n^l$: Number of logarithmic terms of the parametric model

$n^e$: Number of exponential terms of the parametric model

$A_i^l$: Amplitude of the $i$th logarithmic term

$A_i^e$: Amplitude of the $i$th exponential term

$\tau_i^l$: Relaxation time of the $i$th logarithmic term

$\tau_i^e$: Relaxation time of the $i$th exponential term

$t_i^l$: earthquake time(date) corresponding to $i$th logarithmic term

$t_i^e$: earthquake time(date) corresponding to the $i$th exponential term

The variance of $\delta L(t)$ is given by:

$$ \mathrm{var}(\delta L) = C \cdot \mathrm{var}(\theta) \cdot C^T \quad (16) $$

where $\theta$ is the vector of parameters of the post-seismic deformation model:

$$ \theta = \begin{bmatrix} A_1^l, \tau_1^l, \dots, A_{n^l}^l, \tau_{n^l}^l, A_1^e, \tau_1^e, \dots, A_{n^e}^e, \tau_{n^e}^e \end{bmatrix} $$

Note that the values of $\mathrm{var}(\theta)$ of the parametric models are provided to the users in SINEX files available at the ITRF2020 web site.
The elements of the matrix $C = \frac{\partial \delta L}{\partial \theta}$ are computed by the following formulas:

$$ \frac{\partial \delta L}{\partial A_i^l} = \log \left( 1 + \frac{t - t_i^l}{\tau_i^l} \right) \qquad (17) $$

$$ \frac{\partial \delta L}{\partial \tau_i^l} = - \frac{A_i^l (t - t_i^l)}{(\tau_i^l)^2 \left( 1 + \frac{t - t_i^l}{\tau_i^l} \right)} \qquad (18) $$

$$ \frac{\partial \delta L}{\partial A_i^e} = 1 - e^{-\frac{(t-t_i^e)}{\tau_i^e}} \qquad (19) $$

$$ \frac{\partial \delta L}{\partial \tau_i^e} = - \frac{A_i^e (t - t_i^e) e^{-\frac{(t-t_i^e)}{\tau_i^e}}}{(\tau_i^e)^2} \qquad (20) $$

Note that the PSD models are determined and provided to the users per component $L \in \{E,N,U\}$, independently, and so there are NO cross-terms (or correlations) between components. However, cross-terms between amplitude and relaxation time for each LOG or/and EXP term should be taken into account in the variance calculation of Eq. (16). As an example, if for a given station there are 3 earthquakes that were taken into account in the estimation of the PSD models of its component $E$, and it has one EXP for the first EQ, one EXP for the 2nd EQ and LOG+EXP for the 3rd EQ, the one line matrix $C$ for component $E$ in Eq. (16) will have 8 terms.

Once the variances $\text{var}(\delta E)$, $\text{var}(\delta N)$, $\text{var}(\delta U)$ are computed, the transformation into cartesian coordinates is obtained by:

$$ \text{var} \begin{bmatrix} \delta x \\ \delta y \\ \delta z \end{bmatrix} = R \cdot \begin{bmatrix} \text{var}(\delta E) & 0 & 0 \\ 0 & \text{var}(\delta N) & 0 \\ 0 & 0 & \text{var}(\delta U) \end{bmatrix} \cdot R^T \quad (21) $$

## References

Abbondanza C, Altamimi Z, Chin TM, Gross RS, Heflin MB, Parker JW, Wu X (2015) Three-Corner Hat for the assessment of the uncertainty of non-linear residuals of space-geodetic time series in the context of terrestrial reference frame analysis. J Geodesy 89(4):313–329. https://doi.org/10.1007/s00190-014-0777-x

Abbondanza C, Chin TM, Gross RS, Heflin MB, Parker JW, Soja BS, van Dam T, Wu X (2017) JTRF2014, the JPL Kalman filter and smoother realization of the International Terrestrial Reference System. J Geophys Res Solid Earth 122(10):8474–8510. https://doi.org/10.1002/2017JB014360

Altamimi Z, Dermanis A (2012) The Choice of Reference System in ITRF Formulation. In: Sneeuw N, Novák P, Crespi M, Sansò F (eds) VII Hotine-Marussi Symposium on Mathematical Geodesy. Springer, Berlin, Heidelberg, International Association of Geodesy Symposia, pp 329–334. https://doi.org/10.1007/978-3-642-22078-4_49

Altamimi Z, Boucher C, Sillard P (2002) New trends for the realization of the international terrestrial reference system.

Adv Space Res 30(2):175–184. https://doi.org/10.1016/S0273-1177(02)00282-X

Altamimi Z, Sillard P, Boucher C (2004) ITRF2000: From Theory to Implementation. In: Sansò F (ed) V Hotine-Marussi Symposium on Mathematical Geodesy. Springer, Berlin, Heidelberg, International Association of Geodesy Symposia, pp 157–163. https://doi.org/10.1007/978-3-662-10735-5_21

Altamimi Z, Collilieux X, Legrand J, Garayt B, Boucher C (2007) ITRF2005: A new release of the international terrestrial reference frame based on time series of station positions and earth orientation parameters. J Geophys Res Solid Earth. https://doi.org/10.1029/2007JB004949

Altamimi Z, Collilieux X, Métivier L (2011) ITRF2008: an improved solution of the international terrestrial reference frame. J Geodesy 85(8):457–473. https://doi.org/10.1007/s00190-011-0444-4

Altamimi Z, Rebischung P, Métivier L, Collilieux X (2016) ITRF2014: a new release of the International Terrestrial Reference Frame modeling nonlinear station motions. J Geophys Res Solid Earth 121(8):6109–6131. https://doi.org/10.1002/2016JB013098

Altamimi Z, Métivier L, Rebischung P, Rouby H, Collilieux X (2017) ITRF2014 plate motion model. Geophys J Int 209(3):1906–1912. https://doi.org/10.1093/gji/ggx136

Altamimi Z, Rebischung P, Collilieux X, Métivier L, Chanard K (2022) ITRF2020 [Data set]. IERS ITRS Center Hosted by IGN and IPGP. https://doi.org/10.18715/IPGP.2023.LDVIOBNL

Appleby G, Rodríguez J, Altamimi Z (2016) Assessment of the accuracy of global geodetic satellite laser ranging observations and estimated impact on ITRF scale: estimation of systematic errors in LAGEOS observations 1993–2014. J Geodesy 90(12):1371–1388. https://doi.org/10.1007/s00190-016-0929-2

Bird P (2003) An updated digital model of plate boundaries. Geochem Geophys Geosyst. https://doi.org/10.1029/2001GC000252

Blewitt G (2003) Self-consistency in reference frames, geocenter definition, and surface loading of the solid Earth. J Geophys Res Solid Earth. https://doi.org/10.1029/2002JB002082

Blewitt G, Lavallée D (2002) Effect of annual signals on geodetic velocity. J Geophys Res Solid Earth 107(B7):ETG 9-1-ETG 9-11. https://doi.org/10.1029/2001JB000570

Caroli M, de Castro PMM, Loriot S, Rouiller O, Teillaud M, Wormser C (2010) Robust and Efficient Delaunay Triangulations of Points on Or Close to a Sphere. In: Festa P (ed) Experimental Algorithms. Springer, Berlin, Heidelberg, Lecture Notes in Computer Science, pp 462–473. https://doi.org/10.1007/978-3-642-13193-6_39

Collilieux X, Altamimi Z, Coulot D, Ray JR, Sillard P (2007) Comparison of very long baseline interferometry, GPS, and satellite laser ranging height residuals from ITRF2005 using spectral and correlation methods. J Geophys Res Solid Earth. https://doi.org/10.1029/2007JB004933

Collilieux X, van Dam T, Ray J, Coulot D, Métivier L, Altamimi Z (2012) Strategies to mitigate aliasing of loading signals while estimating GPS frame parameters. J Geodesy 86(1):1–14. https://doi.org/10.1007/s00190-011-0487-6

Collilieux X, Altamimi Z, Argus DF, Boucher C, Dermanis A, Haines BJ, Herring TA, Kreeemer CW, Lemoine FG, Ma C, MacMillan DS, Mäkinen J, Métivier L, Ries J, Teferle FN, Wu X (2014) External Evaluation of the Terrestrial Reference Frame: Report of the Task Force of the IAG Sub-commission 1.2. In: Rizos C, Willis P (eds) Earth on the Edge: Science for a Sustainable Planet. Springer, Berlin, Heidelberg, International Association of Geodesy Symposia, pp 197–202. https://doi.org/10.1007/978-3-642-37222-3_25

Collilieux X, Altamimi Z, Rebischung P, Métivier L (2018) Coordinate kinematic models in the International Terrestrial Reference Frame releases. In: Quod Erat Demonstrandum—in quest of the ultimate geodetic insight, Special issue for Professor Emeritus Athanasios Dermanis
Dermanis A (2001) Establishing global reference frames. nonlinear, temporal, geophysical and stochastic aspects. In: Gravity, Geoid and Geodynamics 2000. Springer, pp 35–42. https://doi.org/10.1007/978-3-662-04827-6_6

Dermanis A (2004) The rank deficiency in estimation theory and the definition of reference systems. In: Sansò F (ed) V Hotine-Marussi Symposium on Mathematical Geodesy. Springer, Berlin, Heidelberg, International Association of Geodesy Symposia, pp 145–156. https://doi.org/10.1007/978-3-662-10735-5_20

Dong D, Fang P, Bock Y, Cheng MK, Miyazaki S (2002) Anatomy of apparent seasonal variations from GPS-derived site position time series. J Geophys Res Solid Earth 107(B4):ETG 9-1-ETG 9-16. https://doi.org/10.1029/2001JB000573

Dong D, Yunck T, Heflin M (2003) Origin of the international terrestrial reference frame. J Geophys Res Solid Earth. https://doi.org/10.1029/2002JB002035

Dziewonski AM, Chou T-A, Woodhouse JH (1981) Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. J Geophys Res 86:2825–2852. https://doi.org/10.1029/JB086IB04p02825

Ekström G, Nettles M, Dziewonski AM (2012) The global CMT project 2004–2010: centroid-moment tensors for 13,017 earthquakes. Phys Earth Planet Inter 200–201:1–9. https://doi.org/10.1016/j.pepi.2012.04.002

Gazeaux J, Williams S, King M, Bos M, Dach R, Deo M, Moore AW, Ostini L, Petrie E, Roggero M, Teferle FN, Olivares G, Webb FH (2013) Detecting offsets in GPS time series: first results from the detection of offsets in GPS experiment. J Geophys Res Solid Earth 118(5):2397–2407. https://doi.org/10.1002/jgrb.50152

Gobinddass ML, Willis P, de Viron O, Sibthorpe A, Zelensky NP, Ries JC, Ferland R, Bar-Sever Y, Diament M (2009) Systematic biases in DORIS-derived geocenter time series related to solar radiation pressure mis-modeling. J Geodesy 83(9):849–858. https://doi.org/10.1007/s00190-009-0303-8

Gobron K, Rebischung P, de Viron O, Demoulin A, Van Camp M (2022) Impact of offsets on assessing the low-frequency stochastic properties of geodetic time series. J Geodesy 96(7):46. https://doi.org/10.1007/s00190-022-01634-9

GSC (2022) Galileo Satellite Metadata. European GNSS Service Centre. https://www.gsc.europa.eu/support-to-developers/galileo-satellite-metadata

Harville DA (1977) Maximum likelihood approaches to variance component estimation and to related problems. J Am Stat Assoc 72(358):320–338. https://doi.org/10.1080/01621459.1977.10480998

Hellmers H, Modiri S, Bachmann S, Thaller D, Bloßfeld M, Seitz M, Gipson J (2022) Combined IVS contribution to the ITRF2020. Int Assoc Geod Symp Series. https://doi.org/10.5194/egusphere-egu21-10678

Hellmers H, Modiri S, Thaller D, Gispon J, Bloßfeld M, Seitz M, Bachmann S (2022b) The IVS contribution to ITRF2020. Tech. rep., available at the ITRF2020 website https://itrfr.ign.fr/en/solutions/ITRF2020

IUGG (2007, 2019) IUGG resolutions. https://iugg.org/meetings/iugg-general-assemblies/#5776d8e80445797e0

Johnston G, Riddell A, Hausler G (2017) The International GNSS Service. In: Teunissen PJ, Montenbruck O (eds) Springer Handbook of Global Navigation Satellite Systems, Springer Handbooks. Springer International Publishing, Cham, pp 967–982. https://doi.org/10.1007/978-3-319-42928-1_33

Kreemer C, Hammond WC, Blewitt G (2018) A robust estimation of the 3-D intraplate deformation of the North American plate from GPS. J Geophys Res Solid Earth 123(5):4388–4412. https://doi.org/10.1029/2017JB015257

Luceri V, Pirri M, Rodríguez J, Appleby G, Pavlis EC, Müller H (2019) Systematic errors in SLR data and their impact on the ILRS products. J Geodesy 93(11):2357–2366. https://doi.org/10.1007/s00190-019-01319-w

Métivier L, Collilieux X, Lercier D, Altamimi Z, Beauducel F (2014) Global coseismic deformations, GNSS time series analysis, and earthquake scaling laws. J Geophys Res Solid Earth 119(12):9095–9109. https://doi.org/10.1002/2014JB011280

Métivier L, Altamimi Z, Roubay H (2020) Past and present ITRF solutions from geophysical perspectives. Adv Space Res 65(12):2711–2722. https://doi.org/10.1016/j.asr.2020.03.031

Moreaux G, Stepanek P, Capdeville H, Lemoine FG, Otten M (2022) The DORIS contribution to ITRF2020. Tech. rep., available at the ITRF2020 website https://itrfr.ign.fr/en/solutions/ITRF2020

Nothnagel A, Artz T, Behrend D, Malkin Z (2017) International VLBI service for geodesy and astrometry. J Geodesy 91(7):711–721. https://doi.org/10.1007/s00190-016-0950-5

Pavlis E, Luceri V, Basoni A, Sarrocco D, Kuzmicz-Cieslak M, Evans K, Bianco G (2021) ITRF2020: The International Laser Ranging Service (ILRS) Contribution. ISSN: 1050-9208 Section: Geodesy

Pavlis EC, Luceri V (2022) The ILRS contribution to ITRF2020. Tech. rep., available at the ITRF2020 website https://itrfr.ign.fr/en/solutions/ITRF2020

Pearlman MR, Noll CE, Pavlis EC et al (2019) The ILRS: approaching 20 years and planning for the future. J Geod 93:2161–2180. https://doi.org/10.1007/s00190-019-01241-1

Petit G, Luzum B (2010) IERS conventions (2010)(v1.3.0). Verlag des Bundesamts für Kartographie und Geodäsie Frankfurt 2010

Poyard JC, Collilieux X, Muller JM, Garayt B, Saunier J (2017) IGN best practice for surveying instrument reference points at ITRF co-location sites (No. IERS-TN-39). Verlag des Bundesamts für Kartographie und Geodäsie Frankfurt 2017

Ray JR (1996) Measurements of length of day using the global positioning system. J Geophys Res Solid Earth 101(B9):20141–20149. https://doi.org/10.1029/96JB01889

Ray JR (2009) A quasi-optimal, consistent approach for combination of UT1 and LOD. In: Drewes H (ed) Geodetic reference frames: IAG symposium Munich, Germany, 9-14 October 2006, International Association of Geodesy Symposia. Springer, Berlin, Heidelberg, pp 239–243. https://doi.org/10.1007/978-3-642-00860-3_37

Ray JR, Altamimi Z, Collilieux X, van Dam T (2008) Anomalous harmonics in the spectra of GPS position estimates. GPS Solut 12(1):55–64. https://doi.org/10.1007/s10291-007-0067-7

Rebischung P (2020) IGS Reference Frame Working Group Coordinator Report (2019). In: Villiger A, Dach R (eds) International GNSS service: technical report 2020, Bern Open Publishing. IGS Central Bureau and University of Bern, Bern. https://doi.org/10.48350/156425

Rebischung P (2022) The IGS contribution to ITRF2020. Tech. rep., available at the ITRF2020 website https://itrfr.ign.fr/en/solutions/ITRF2020

Seitz M, Bloßfeld M, Angermann D, Seitz F (2022) DTRF2014: DGFI-TUM’s ITRS realization 2014. Adv Space Res 69(6):2391–2420. https://doi.org/10.1016/j.asr.2021.12.037

Sillard P, Boucher C (2001) A review of algebraic constraints in terrestrial reference frame datum definition. J Geodesy 75(2):63–73. https://doi.org/10.1007/s001900100166

van Dam TM, Wahr J (1998) Modeling environment loading effects: a review. Phys Chem Earth 23(9):1077–1087. https://doi.org/10.1016/S0079-1946(98)00147-5

Williams SDP, Willis P (2006) Error analysis of weekly station coordinates in the DORIS network. J Geodesy 80(8–11):525–539. https://doi.org/10.1007/s00190-006-0056-6
Williams SDP, Bock Y, Fang P, Jamason P, Nikolaidis RM, Prawirodirdjo L, Miller M, Johnson DJ (2004) Error analysis of continuous GPS position time series. J Geophys Res Solid Earth. https://doi.org/10.1029/2003JB002741

Willis P, Lemoine FG, Moreaux G, Soudarin L, Ferrage P, Ries J, Otten M, Saunier J, Noll C, Biancale R, Luzum B (2016) The International DORIS Service (IDS): Recent Developments in Preparation for ITRF2013. In: Rizos C, Willis P (eds) IAG 150 Years. Springer International Publishing, Cham, International Association of Geodesy Symposia, pp 631–640. https://doi.org/10.1007/1345_2015_164

Wu X, Abbondanza C, Altamimi Z, Chin TM, Collilieux X, Gross RS, Heflin MB, Jiang Y, Parker JW (2015) KALREF-A Kalman filter and time series approach to the international terrestrial reference frame realization. J Geophys Res Solid Earth 120(5):3775–3802. https://doi.org/10.1002/2014JB011622

Xu X, Dong D, Fang M, Zhou Y, Wei N, Zhou F (2017) Contributions of thermoelastic deformation to seasonal variations in GPS station position. GPS Solut 21(3):1265–1274. https://doi.org/10.1007/s10291-017-0609-6