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

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.

erence 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(IUGG , 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 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 ITRS/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).
All 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).

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.

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(Altamimi et al. , 2004Dermanis 2001Dermanis , 2004Sillard 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 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 combination 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 stations. The solutions include the full set of EOPs, in addition to station positions.

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 agencies operating ITRF colocation sites to provide relative coordinates (local ties) connecting the instrument reference points. Figure 1 illustrates the full ITRF2020 network, comprising 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 instruments 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 combination 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.

ITRF2020 analysis strategy
Among a number of other innovations and analysis improvements, 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 constraints 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 station position time series and determine parametric models for stations subject to post-seismic deformation.

Time series analysis
Although the full equations of the combination model implemented in the CATREF Software that is used for the ITRF computations were published in past articles, see for instance , they are detailed in the appendix of this stand-alone paper for the convenience of the reader.
In preparation for the global overall ITRF2020 combination described in Sect. 4, intermediate analysis steps and stacking of the individual technique time series were conducted, including the following:   -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: where i1 and i2 are two close-by stations or two consecutive segments (SOLNs) of a station with discontinuities.Ẋ i1 andẊ i2 are the corresponding velocities.
T are the coefficients of the cosine and sine functions describing the periodic station displacements at frequency j (see Appendix A for more details). σ v and σ 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.

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 longterm 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.

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.

Post-seismic deformation
Initiated at the occasion of the ITRF2014 , 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)arithmic, (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.   (Dziewonski et al. 1981;Ekström et al. 2012) and ITRF2020 sites (green) impacted by Post-Seismic Deformation 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.

ITRF2020 global combination
Past ITRF combinations up to ITRF2014 were operated using a two-step procedure (Altamimi et al. , 2011: (1) stacking the individual time series to estimate a longterm 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 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 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 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 (σ v and σ 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 pres- 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 Subsect. 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 σ v and σ f for the velocity and seasonal signals constraints, respectively, were chosen as described in Subsect. 4.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 σ 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 σ 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 σ 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 σ v of Eq. 1 was increased for sites where the velocity normalized residual exceeded the threshold of 3 in any of the three components.

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 semiannual) of the same chosen couples of GNSS and reference points of the colocated instruments that were constrained to have the same velocity. The initial σ 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, σ 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.

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.

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.

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.

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.

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.

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 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.
where the couple (X i (t 0 ),Ẋ i ) are the station position and velocity for segment (SOLN) i bounded by dates (t i , t i+1 ), δ X PSD (t) is the total sum of the PSD corrections-the PSD corrections of successive earthquakes are cumulative, and δ X i f (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 δ X PSD (t) and δ X i f (t).

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, δ X P S D (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.

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 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 δ X f (t) using Eq. 4 and add it to the station coordinates at epoch t, using Eq. 3.
where t is the decimal year and (a i x , a i y , a i z , b i x , b i y , b i z ) 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  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 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 CF/CM . The obtained ITRF2020 seasonal coefficients of vector ΔX CF/CM 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 CF/CM along x, is obtained by Eq. 5.
where t is the time in decimal years, A 1 and A 2 are the annual and semi-annual amplitudes, φ 1 and φ 2 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).

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 f (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 (interannual, 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 semiannual 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 Δ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 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- 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.

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 . 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).

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 stationdependent 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  where the agreement was at the level of 1.39 ppb (about 8 mm at the equator).  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.  (Bird 2003)

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 Antarctica. 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).  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.

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(Ray , 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-arcseconds, 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.

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).
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 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.

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 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 that 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 : 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 , 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.

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://creativecomm ons.org/licenses/by/4.0/. and L O D 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  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: where P i = [P i,1 , . . . , 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:

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 ).

C.1 Periodic signals: minimum constraints
In the stacking case, the equations of the minimum constraints approach take the form where A is the well known design matrix of partial derivatives of the 7 transformation parameters, and (a i x R , a i y R , a i z R , b i x R , b i y R , b i z R ) 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:

D Equations of post-seismic deformation models
After an earthquake, the position of a station during the postseismic trajectory at an epoch t is described by an additional term δ X P S D (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: where R is the transformation (Jacobian) matrix from topocentric to geocentric frame. For each component L ∈ {E,N,U}, the total PSD displacement in the local frame can be computed by: where: n l : Number of logarithmic terms of the parametric model n e : Number of exponential terms of the parametric model A l i : Amplitude of the ith logarithmic term A e i : Amplitude of the ith exponential term τ l i : Relaxation time of the ith logarithmic term τ e i : Relaxation time of the ith exponential term t l i : earthquake time(date) corresponding to ith logarithmic term t e i : earthquake time(date) corresponding to the ith exponential term The variance of δL(t) is given by: where θ is the vector of parameters of the post-seismic deformation model: θ = A l 1 , τ l 1 , ...., A l n l , τ l n l , A e 1 , τ e 1 , ...., A e n e , τ e n e Note that the values of var(θ ) of the parametric models are provided to the users in SINEX files available at the ITRF2020 web site.
The elements of the matrix C = ∂δL ∂θ are computed by the following formulas: ∂δL Note that the PSD models are determined and provided to the users per component L ∈ {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 var(δ E), var(δ N ), var(δU ) are computed, the transformation into cartesian coordinates is obtained by: