1 Introduction

Modern geometric reference frames (RFs) are realized through the epoch-specific, time-dependent coordinates assigned to the reference stations used to define the frame. Those coordinates establish, by implication, the scale, orientation, and origin of the frame. Geometric global RFs (GRFs), such as the International Terrestrial Reference Frame 2014 (ITRF, Altamimi et al. 2016), can be accessed and utilized for regional geodesy. Yet, there are many advantages to densifying global frames on a regional basis. Regional frames can incorporate many additional continuous operating reference stations (and passive benchmarks) within each region of interest, allowing surveyors and others to position over much shorter baselines. Denser regional reference frames also allow scientists and geodetic authorities to refine their models for the crustal velocity field pervading each region. This is particularly useful in regions affected by the earthquake deformation cycle. However, regional reference frames (RRFs) are far more useful if they are geometrically consistent with a standard GRF, and when both frames are as consistent as possible with physical reality. Attaching many dense regional or national RFs to a standard GRF provides us with a distributed processing strategy for accommodating very large numbers of geodetic reference stations worldwide.

RRFs are usually regarded as ‘densifications’ or regional ‘implementations’ of the GRF (International Organization for Standardization 2020). Although the use of Precise Point Positioning (PPP) techniques is on the rise, this work discusses a densification method that uses more accurate differential processing. Using differential or ‘network’ solutions, one can achieve this densification by ‘aligning’ the daily or weekly regional solutions onto the GRF’s epoch-specific coordinates, or by aligning a long-term stacked regional solution with the GRF. We may view a RF as a network trajectory model consisting of a suite of station trajectory models, one for each station in the reference network (Bevis and Brown 2014). Geometrical inconsistencies between the trajectory models used to define a frame and the reference network’s actual evolving geometry constitute one of the major sources of reference frame realization error. Thus, in attaching a RRF to a GRF, the trajectory models associated with the global and regional RFs should be as consistent as possible.

To date, one particular inconsistency in RRF realization is the failure to incorporate periodic motions that are entirely consistent with the displacement oscillations observed in the GRF. Station oscillations (i.e., seasonal displacements) are produced by effects such as changes in the position of the center of mass of the Earth (relative to the crust), elastic deformation of the crust and the entire Earth driven by seasonal changes in the loads imposed by the hydrosphere, atmosphere, and cryosphere—and pole tides. Also, apparent periodic motions can be generated by non-physical errors such as processing artifacts. RRFs are often (if not always) realized using regional GPS or GNSS-only differential solutions that cannot sense the absolute values of displacement oscillations. Differential GNSS solutions can sense periodic displacements only in relative terms, wherein a common-mode oscillation is ignored (Blewitt 2003; Collilieux et al. 2010, 2012). As a result, the ‘disconnect’ between the periodic station displacements in the RRF and the GRF yields a regional realization that is likely to have periodic biases relative to the GRF.

For many years various versions of ITRF were defined and realized using trajectory models that did not include annual or semiannual periodic displacements. As oscillatory displacements were resolved with greater accuracy, it was recognized that adding seasonal terms to station trajectory models increases the stability and accuracy of station velocity estimates (e.g., Blewitt and Lavallée 2002; Dong et al. 2002; Gómez et al. 2015). This stability and accuracy increase led Bevis and Brown (2014) to suggest that new classes of trajectory models that incorporated annual periodic displacements or oscillations (and other features) should be used to define RFs, not just to characterize the displacements of stations expressed in those frames. The greater geometrical consistency of RFs defined using modern, generalized trajectory models prompted Altamimi et al. (2016) to add seasonal terms to the equations underlying ITRF 2014. Although the parameters or coefficients of the periodic components of their trajectory models were not published along with all the other trajectory model parameters, these Fourier coefficients are available from the authors. We used those coefficients in this study.

ITRF 2014 provides a framework that promotes greater consistency between global and regional RFs. Before ITRF included periodic terms, several authors augmented previous versions of the ITRF (e.g., Zou et al. 2014) or tested the impact of ignoring the periodic motions on regional and global reference frame alignment (e.g., Collilieux et al. 2010, 2012). Yet, these studies proposed to mitigate the effects of periodic motions using dynamic approaches by introducing loading model corrections or estimates of displacements using other techniques (e.g., GRACE). Our goal is to realize a RRF that does not suffer from biases due to unmodeled periodic motions using a purely kinematic approach that does not require any a priori knowledge of, or geophysical models for periodic motions. We describe a method that ‘propagates’ the regional periodic motions incorporated in ITRF 2014 into the GPS-only RRF Posiciones Geodésicas Argentinas (POSGAR), realized at the National Geographic Institute of Argentina (Instituto Geográfico Nacional, IGN), without using geophysical models or non-GNSS observation techniques. We will show that the ‘propagation’ of periodic terms can be thought of in the same way as a position and velocity space alignment of a stack of polyhedrons. Application of our technique increased the consistency between the POSGAR time series and those defining ITRF 2014 by reducing the magnitude of the periodic component discrepancy of several stations from ~ 4 mm down to less than 1 mm in the vertical component. To do this we must extend the concept of alignment to include periodic coefficients and behavior. To formalize and generalize this extended alignment concept, we will borrow some terminology from object-oriented programming (OOP) to extend the International Terrestrial Reference System ISO standard (ISO 19161–1:2020). Our goal is to improve the terminology that describes the hierarchical relationship between RFs, which, in our opinion, is only loosely defined in the geodetic community.

As we explain below, we ensure that a RRF is optimally attached to a GRF by using compatible trajectory models at the set of reference stations that are common to the global and regional networks. In effect, the station trajectory models specified in the GRF are imposed on or ‘inherited by’ those same stations in the geodetic analysis that realizes or defines the RRF. There are multiple sources of RRF realization error, including the so-called ‘network effect’ (Legrand and Bruyninx 2009). We ensure that such problems, both physical and procedural, do not unduly affect our alignment of the regional and global RFs by posterior comparison of all common station trajectories.

2 Hierarchical relationships between geodetic reference frames

RFs have a hierarchy (sometimes implicit), which is well known to geodesists. This hierarchy is described in a recent document released by the International Organization for Standardization, ‘The International Terrestrial Reference System (ITRS) ISO standard 19,161–1:2020.’ Nevertheless, the relationship between two RFs can lead to some confusion within the user community. Therefore, this section provides a few definitions to extend the jargon associated with geometric RFs.

For example, to be useful to the engineering or scientific communities, a national geodetic RF most likely depends on another frame, higher in the hierarchy and, thus, of lower order. ISO 19161–1:2020 describes these two RFs, low- and high-order, as a ‘primary’ and ‘secondary’ realization of the ITRS. Using the ISO terminology, a national RF (i.e., secondary) can be referred to as a ‘densification’ or ‘implementation’ of a primary frame. According to the ISO standard, one achieves such implementation by aligning the secondary onto the primary. There are at least two methods to achieve this alignment. One method is to repeatedly use a six- or seven-parameter Helmert transformations of the secondary’s daily (or weekly) solutions at each and every epoch of interest. We use H to denote Helmert, followed by the number of transformation parameters that define it, in this case, H6 or H7. Another methodology is to project an extended time series ‘stack’ of regional solutions (polyhedrons), already aligned in space and time onto the model trajectories of the primary RF using an appropriate generalized Helmert transformation (e.g., H12, H14, or higher, as we will show). In this paper, we concentrate on the second method.

To allow more flexibility in describing the relationship between frames, we propose an augmentation to ISO 19161–1:2020 with terms and concepts commonly used in object-oriented programming (OOP). These concepts came to us while attempting to realize a secondary RF at IGN. We seemed to be describing relationships between frame parameters that had never been described before, at least to the best of our knowledge. The application of the OOP concept of inheritance to geometric geodesy seemed to fit the things that we were doing.

In OOP, a ‘child’ (or secondary) object can inherit the properties and functions (or methods) of a ‘parent’ (or primary) object, allowing modifications (called ‘overrides’ in OOP) to the properties and functions in the child object. Therefore, we conceptualize a ‘parent’ or primary geodetic RF as an object with a ‘child’ or secondary RF that inherits all (or some) of its properties or parameters. It is easy to see that in geometric geodesy, the concept of inheritance may be used interchangeably with the notion of alignment, either in space (3D) or in space and time (4D). Nevertheless, we will describe how the concept of inheritance allows us to be more flexible in describing the parameter relationship between primary and secondary RFs.

Current geometric geodesy jargon implies that RF parameters can be inherited from the position and velocity spaces of a parent RF. The position space contains vectors representing spatial coordinates or reference positions at some reference epoch. These coordinates imply a specific scale, origin, and orientation. The linear component of the temporal evolution of these coordinates constitute station velocities, and these velocity vectors reside in velocity space. Position and velocity space inheritance means that the secondary RF will inherit all these parameters: orientation, origin, scale, and velocity from the primary RF.

As the ISO standard describes, the RF (global or regional) spatial scope does not necessarily define a primary–secondary relationship. One example is realizing a ‘pure-GNSS’ global RF based on the ITRF. In such a case, position and velocity inheritance is invoked to align the secondary (pure-GNSS) RF with the primary. Conceptually speaking, inheritance with an override allows the secondary RF to have different trajectory models for its stations (although this may not be ideal) while aligning the position and velocity parameters as close as possible to those of the primary. By way of example, a RRF definition might modify a station trajectory model to include a logarithmic transient that mimics a postseismic transient produced by an earthquake that occurred after the GRF was formulated and published.

Now that we have presented the concept of inheritance, we expand it by introducing more properties or behaviors not covered by the ISO standard and, to our knowledge, not discussed within the geodetic community: those developed in oscillation, frequency, or Fourier space. As discussed above, behavior represented by coefficients in velocity space is inherited by, or imposed on a secondary frame, using H transformations. These H transformations are designed to constrain the station velocities referred to the secondary RF to lie as closely as possible to those expressed in the primary RF (at all common stations). Similarly, we can inherit periodic behavior using H transformation parameters that constrain oscillations expressed in the secondary RF to match those expressed in the primary RF by minimizing a penalty function that gauges disagreement between the two sets of Fourier coefficients.

For better readability, we want to clarify a few aspects for readers unfamiliar with stack realizations performed using H transformations in position and velocity spaces, which we are now extending to include frequency space. Each space can be accessed only through the parameters of the stations that realize the frame. Because each space (position, velocity, and frequency) is orthogonal to the others, the H transformations used to inherit their parameters can be applied sequentially. We should clarify that the H transformations operating in velocity and frequency space are essentially independent if the stations’ time series do not contain large unmodeled jumps, given a time series of sufficient duration (~ 3 years). When these conditions are satisfied, the Fourier coefficients and the velocity coefficients for each station are very weakly correlated. If we wish to inherit the position and velocity parameters simultaneously, we can use a H14 or H12 transformation. However, we may also apply two separate H7 or H6 transformations, one for the position inheritance and the other for the velocity inheritance. Thus, later on, when we work with frequency space H transformation, we will work only with the Fourier coefficients necessary to drive this class of inheritance, temporarily ignoring the velocity and position coefficients.

Ignoring periodic behavior when realizing a RRF can have undesirable consequences, e.g., periodic RRF realization biases. The following section shows the consequences of disregarding periodic behavior in the most recent ITRF and the Geodetic Reference System for the Americas (Sistema de Referencia Geodésico para las Américas, SIRGAS) realizations. Then, in Sect. 4, we analyze a GPS solution stack obtained using GAMIT/GLOBK (Herring et al. 2010) for part of the SIRGAS continuous GNSS (CGNSS) reference network (SIRGAS-CON, Sánchez and Drewes 2020) in Central and South America, Antarctica, and the Caribbean.

3 ITRF common-modes in South America and the SIRGAS periodic bias

A RRF realized solely using GNSS observations is incapable of capturing any effects considered common mode in the network. Such effects can be large-scale crustal elastic deformation, pole tide effects, displacements of the center of mass of the Earth if we wish to realize a RRF tied to such center (see Blewitt 2003 for a review on RF origins), or even apparent displacements due to processing artifacts. This common-mode omission occurs because we use differential processing techniques to produce the GNSS polyhedrons. Thus, the GNSS solution will not discern or resolve any displacement component that is common to all stations.

Periodic variations of station coordinates can be introduced to the station trajectory models by incorporating the following term (Bevis and Brown 2014):

$$s=\sum_{k=1}^{{n}_{f}}[{s}_{k}\mathrm{sin}({\omega }_{k}t)+{c}_{k}\mathrm{cos}({\omega }_{k}t)],$$

where \({n}_{f}=2\) for annual and semiannual components and \({s}_{k}, {c}_{k}\) are the corresponding amplitudes of the \({k}^{th}\) frequency component. We fitted the SIRGAS time series to obtain the \({s}_{k}, {c}_{k}\) periodic amplitudes because the SIRGAS realization does not include periodic components. Although the effect of ignoring periodic motions in RF realization has been investigated before (Collilieux et al. 2010, 2012; Zou et al. 2014), we analyzed this effect in the context of a RRF in our region of interest. To determine if any biases are present in the \({s}_{k}, {c}_{k}\) SIRGAS coefficients (using the multi-year SIR17P01 solution), we compared our fit of \({s}_{k}, {c}_{k}\) of the SIRGAS-CON time series against those of ITRF 2014 at some of the CGNSS stations that correspond to the intersection between the two frames. Since SIR17P01 is aligned to ITRF 2014 through IGS14 (see Sánchez and Drewes 2020), a RRF user would expect that their periodic terms would be nearly identical. Any periodic term difference larger than zero reveals a discrepancy of the periodic behavior between SIR17P01 and ITRF 2014.

Figure 1a shows a trajectory fit including periodic terms in SIR17P01 to station CHPI (Cachoeira Paulista, Brazil) and the magnitude of the difference between the periodic amplitudes of ITRF 2014 and SIR17P01 in both the vertical and horizontal components (Fig. 1b). For convenience, instead of showing the difference of \({s}_{k}, {c}_{k}\) separately, Fig. 1b shows the magnitude of the difference between the annual parameters in the north, east, and up, computed as

Fig. 1
figure 1

a Trajectory model of SIRGAS time series for CHPI incorporating seasonal terms. b The seasonal terms adjusted for the SIRGAS-CON RF and ITRF 2014 at eight stations common to both frames

Full size image
$$m=\sqrt{{\left({s}_{k}^{ITRF}-{s}_{k}^{SIRGAS}\right)}^{2}+{\left({c}_{k}^{ITRF}-{c}_{k}^{SIRGAS}\right)}^{2}}$$
(1)

The value computed using (1) is only zero when both the magnitude and phase of the periodic motions in the RFs are equal. For CHPI, m reaches almost 1 mm, while other stations show biases above 2 mm as observed in the vertical of stations BRAZ (Brasilia, Brazil) and CRO1 (Saint Croix, Virgin Islands). In turn, these results support previous studies showing that unmodeled seasonal signals are aliased into the H transformation parameters during RRF (and GRF) realizations (e.g., Collilieux et al. 2012; Zou et al. 2014).

The magnitude of the bias shown here ranges from ~ 0.5 to 2 mm even after processing hundreds of stations well-distributed across Central and South America, Antarctica, and the Caribbean. We should mention that increasing the extent of the network should minimize the periodic bias due to local loading effects, as long as the common stations (between SIR17P01 and ITRF 2014) are not located within the same area. For example, if a network is limited to the Amazon basin, then the bias’ amplitude would be much larger and more harmful to the users trying to access the ITRF using this hypothetical regional realization (see the discussion in Zou et al. 2014). In this example, the periodic bias is mainly due to the well-known periodic loading and unloading of the crust in the Amazon (Bevis et al. 2005; Galván et al. 2018).

To better understand the effect of common-mode oscillations, we analyzed the periodic terms of ITRF at a set of 27 sites located in our region of interest. We do this by estimating the parameters of an H7 transformation that operates in frequency space, i.e., we find the transformation parameters that minimize the amplitudes of \({s}_{k}, {c}_{k}\) in ITRF (for the annual and semiannual components) for the 27 sites. If the transformation parameters yield nonzero values, then common-mode periodic signals in the ITRF are present at the network-scale level (for the specific set of stations).

Figure 2 shows the translation and scale components of the previously described H7 transformation results for the annual and semiannual components using phasorial plots. We do not show rotations since we found this common-mode signal less significant because the periodic rotation components at the Earth’s surface were < 1 mm. To demonstrate that the common modes obtained in this analysis are not an artifact due to a trajectory fit at a noisy station, we computed 27 H7 transformations by removing one station at a time. We also computed the transformation that includes all stations, shown in blue in Fig. 2. The 27 H7 transformations show a phase variation of about ± two weeks, but no significant amplitude difference (standard deviation < 0.2 mm).

Fig. 2
figure 2

Phasorial plots of common modes observed in ITRF 2014 using 27 Central, South America, Antarctica, and Caribbean stations. Green and red arrows represent the 27 common-mode components (annual and semiannual), leaving one station out of the analysis at the time; blue arrows represent the phasor computed by using all stations. A significant component exists in the X and Z components of annual translation. An apparent seasonal scale change is also observed for the annual component

Full size image

As shown in Fig. 2, all signals related to semiannual components are < 0.2 mm. Therefore, we will not discuss semiannual signals in any detail, but we show how our method improves the RRF-ITRF agreement of this frequency as well. As we would expect for the annual component, which shows summer–winter differences in surface loading (that have opposing signs across the equator), there is a significant periodic translation in the Z direction with its peak between December and January. The translation is also significant in the X-axis, with a peak between July and August. The periodic scale change in Fig. 2d should be considered apparent since the network only spans a portion of the Earth. Therefore, some translation signal leaks into the RRF network’s scale changes (due to the data noise). Nevertheless, we cannot ignore the possibility of an apparent periodic scale change due to variations from unmodeled tropospheric water vapor. Whether we include this scale change or not (i.e., use an H7 or H6 transformation during inheritance) depends on the application we are pursuing: In this case we wish to match ITRF as closely as possible, so we include it during the inheritance process.

We have established and quantified a bias in the SIRGAS-CON RF that is aligned with, or attached to ITRF 2014 without considering the periodic variations of station coordinates. The following sections will describe the RRF realized at the IGN and the method we used to inherit the frequency behavior from ITRF 2014.

4 Description of the IGN solution stack

The IGN is responsible for defining Argentina’s official geodetic RF, known as POSGAR. This section will discuss the most recent realization that we call POSGAR07b or just P07b. The United Nations Global Geospatial Information Management (UN-GGIM) subcommittee on geodesy recommends that ITRF be the RF for geospatial and scientific applications. With this in mind, all the available GPS data since the creation in 1998 of IGN’s CGNSS network, RAMSAC (Piñón et al. 2018), was reprocessed to realize P07b by inheriting key properties from appropriate ITRF trajectory parameter spaces.

Figure 3 shows the ~ 1100 stations processed from 1998 to 2019, representing a total of ~ 1.5 million station-days (RINEX files). The network by design included all stations in Argentina and many neighboring countries, with a particular interest in Chile (Báez et al. 2018) and Brazil (Fortes et al. 2012). The westernmost part of Argentina is affected by the non-steady-plate boundary deformation associated with the earthquake cycle near the Nazca-South America plate boundary, and slower, permanent tectonic deformation of the entire Andes. Thus, the Chile stations provide data about the near-field deformation related to past and future earthquakes at the plate boundary. In contrast, Brazil stations provide data from the South American plate’s stable core, constraining any jumps due to seismic events in or near the Andes.

Fig. 3
figure 3

Map of the stations (continuous and campaign) used to realize POSGAR07b processed by the Scientific Processing Center – Argentina (Centro de Procesamiento Científico – Argentina) IGN. Red dots represent the reference stations common to POSGAR, SIRGAS, and the ITRF reference frames

Full size image

To process all the data shown in Fig. 1, we developed a parallelized Python wrapper for GAMIT/GLOBK (Parallel.GAMIT, available through GitHub). We used a three-node Linux cluster with 112 cores, which allows us to process GPS data ~ 100 × faster than serial processing. This reprocessing of GPS data used the latest orbits and antenna calibration parameters available from the IGS, the Vienna Mapping Functions (Boehm et al. 2006) to estimate the atmospheric delay, and the ocean tide loading model FES2014b (Lyard et al. 2020).

To realize P07b, rather than just aligning each daily polyhedron to the primary RF, we realized our stack using extended trajectory models (ETMs) as described by Bevis and Brown (2014) and Bevis et al. (2019), which among other terms, incorporates periodic variations of station coordinates. As in the latest version of the ITRF, we adopted one annual and one semiannual seasonal mode at all stations with sufficient data to constrain these coefficients. We note that our stacking method does not require any constraints at this point (as in the normal equation stacking method) because it uses the iterative approach described by Bevis and Brown (2014). The weighted root mean square (wrms) misfit between the stacked solutions of P07b (we call this stack SP07b, before inheriting parameter space behavior from ITRF 2014) and the trajectory models show a scatter that varies as a function of time, as shown in Fig. 4. There are two prominent variations in the wrms: (1) a visible reduction in the mean wrms as a function of time, described by an adjusted polynomial in Fig. 4; (2) a periodic variation of the wrms that peaks in the southern hemisphere summer season. These two effects have different origins. The reduction of the wrms is mainly due to incorporating more stations into the network. However, some reduction during the early 2000s is related to GPS improvements (more satellites and improvements to receiver technology). Regarding the seasonal effect, the polyhedrons’ alignment to the trajectory models tends to show a larger disagreement during the southern hemisphere’s summer season due to a larger scatter in the GPS solutions. This additional scatter is mostly caused by unmodeled propagation delays due to water vapor, which are both larger and more variable in summer.

Fig. 4
figure 4

Total weighted root mean square (wrms) of the stack with respect to the trajectory models; a as a function of time also showing the adjusted curve (solid red line) to the wrms, and several stations as a function of time (solid orange line); b histogram of wrms showing the peak wrms at ~ 2.6 mm

Full size image

We also show the wrms as a function of the component (north, east, and up) to demonstrate the performance of SP07b (Fig. 5). As expected, the vertical shows a larger average wrms, typical of GNSS observations. As shown in Figs. 4 and 5, SP07b shows an excellent wrms scatter resulting from the consistency of the models used to realize the RF with respect to station coordinates’ observed behavior. This high consistency is reflected by the small overall variation of wrms scatter (~ 2 mm peak-to-peak after 2006) shown in Fig. 4.

Fig. 5
figure 5

North (a), east (b), and up (c) weighted root mean square (wrms) of the stack with respect to the trajectory models

Full size image

The high consistency and low wrms of SP07b, shown in Fig. 4, are achieved by analyzing the stack and searching for metadata blunders and other undetected mechanical jumps in the time series. Our stacking algorithm also incorporates any earthquakes that might have introduced geophysical jumps into the models. If we later determine that the earthquake had no significant effect on a specific station, we remove that jump from that station’s model. In this way, we can assure that our GPS solutions’ inner geometry is always honored by the resulting stack, avoiding biases and jitter during the frame realization (Bevis and Brown 2014). The stack’s position and velocity space parameters are arbitrary in terms of their absolute values (i.e., their ‘outer geometry’ is arbitrary) because we stack the polyhedrons iteratively onto the stations’ trajectory models without any prior constraints. Only their relative or ‘inner geometry’ is meaningful at this stage. Bevis and Brown (2014) describe the trajectory modeling, inner and outer geometry, and the stacking process in greater detail.

We first discuss how the P07b realization inherits the position and velocity parameters of ITRF 2014 using the set of stations common to both frames to perform the corresponding H14 transformation. We discuss the frequency space inheritance in the next section. We removed some stations during the inheritance process due to poor quality time series or large discrepancies in trajectory fits between SP07b and ITRF 2014. Including these stations would have degraded our final solution. Therefore, we used 24 stations to inherit the position parameters of ITRF 2014 and 18 to inherit the velocity parameters. Figure 6 shows the residuals after inheritance. The total wrms misfit of the position inheritance was 1 mm, while the wrms misfit of the velocity inheritance was 0.20 mm/yr. Some stations show velocity discrepancies above 0.4 mm/yr (e.g., KOUR and PARC), probably from a weaker determination of the polyhedron geometry at the edge of the network. These discrepancies are similar to those described by Legrand and Bruyninx (2009) as a ‘network effect,’ although the magnitude of the differences found by us are well below 1 mm/yr.

Fig. 6
figure 6

North, east, and up residuals in position (a) and velocity space (b)

Full size image

Like the position and velocity parameters, the frequency parameters of SP07b have an arbitrary alignment due to the lack of constraints. To avoid large spurious oscillations in our RF, which make the visual inspection of the stack hard, we develop our working or provisional solutions by removing common-mode oscillations from SP07b and applying the internal constraints described in Appendix B of Altamimi et al. (2016) for the ITRF 2014 realization. The goal of the common-mode removal is purely operational and meant to obtain time series with reasonably sized oscillations. Regardless, when applying internal constraints, the results will have biases with respect to the GRF because of the network’s regional aperture. The use of inheritance reduces these biases and achieves a better agreement between the frequency parameters of P07b and ITRF 2014.

5 Inheriting parameters in frequency space

Inheriting the position and velocity parameters through an H7 or H14 transformation (H6 or H12 if no scale factor is used) is well understood and widely practiced. To inherit the frequency parameters of ITRF 2014, we first determine the difference between the periodic coefficients of SP07b (resulting from the internal constraints) and ITRF 2014 at the common stations (Fig. 7a and b). It should be noted, however, that the common-mode removal step (internal constraints) can be omitted if a ‘provisional’ stack is not needed. As in Fig. 1, we show these residuals (the difference between the periodic coefficients of SP07b and ITRF 2014) as periodic displacement magnitudes in the north, east, and up (NEU) directions. In some cases, residuals reach ~ 4 mm in the vertical (Fig. 7a). These NEU \({s}_{k}, {c}_{k}\) residuals are converted to the Earth Centered-Earth Fixed (ECEF) coordinate system to find the H7 transformation that minimizes the secondary RF periodic residuals, expressed as geocentric \({s}_{k}, {c}_{k}\) (for both the annual and semiannual components). To invert for the parameters of the H7 transformation than transforms \({s}_{k}\) of SP07b (\({s}_{k}^{SP07b}\)) into \({s}_{k}\) of ITRF (\({s}_{k}^{ITRF})\), thus minimizing their difference, we use least squares:

$$H7({s}_{k}^{SP07b} |{s}_{k}^{ITRF})={{(A}^{t}PA)}^{-1}{A}^{t}PL$$

where \(A\) is the design matrix for the linearized Helmert transformation, \(P\) is the weights matrix, and \(L\) contains the differences between \({s}_{k}^{SP07b}\) and \({s}_{k}^{ITRF}\) expressed in ECEF. This H7 transformation aims to make these two coefficient vectors, \({s}_{k}^{SP07b}\) and \({s}_{k}^{ITRF}\), as close as possible to equal. The same process is repeated for \({c}_{k}\). The resulting H7 transformations are then applied to SP07b by subtracting the amplitudes to effectively achieve inheritance as:

$${X}_{i}=X- \mathrm{sin}\left({\omega }_{k}t\right) A Hs7({s}_{k}^{SP07b} |{s}_{k}^{ITRF})-\mathrm{ cos}\left({\omega }_{k}t\right) A Hc7({c}_{k}^{SP07b} |{c}_{k}^{ITRF})$$

where \(X\) are the original coordinates of SP07b, \({X}_{i}\) are the coordinates after inheritance, Hs7 is the H7 transformation for \({s}_{k}\), and Hc7 is the H7 transformation for \({c}_{k}\). This process is repeated for each k frequency. The inversion for the H7 transformation parameters uses a robust least-squares algorithm so as to avoid undesired distortions due to major outliers.

Fig. 7
figure 7

Annual and semiannual residuals between SP07b and ITRF 2014. a Residuals before inheriting the frequency space of ITRF 2014. b Same as a for semiannual. c Residuals after inheriting the seasonal coefficients of ITRF 2014. d Same as c for semiannual

Full size image

Figures 7c and d show the results of the process described above. In most cases, NEU magnitudes computed using Eq. (1) are well below 1 mm for the annual component and below 0.4 mm for the semiannual, although some outliers can be observed. We investigated these and determined that for Bahía Blanca, Argentina (VBCA), the ITRF 2014 solution shows large excursions in the east component, which we successfully filtered in our processing. In other words, the outliers in VBCA altered the periodic coefficients in ITRF 2014, leading to the large residuals observed in Fig. 7a–d. The same occurs with Fortaleza, Brazil (BRFT)’s vertical component, as shown in Fig. 7c and d. Thanks to the robust least-squares algorithm, these two outliers have almost no effect in the final solution, as shown in Fig. 7.

We take as a case study stations BRAZ and CHPI, where we showed that the SIRGAS-CON RF has a vertical residual in frequency space of up to 2 mm with respect to ITRF 2014. Our results show a residual below 0.5 mm for the same two sites after inheritance, i.e., after alignment of our RRF and ITRF in Fourier space. Consequently, a user accessing the ITRF 2014 through P07b should observe the same vertical and horizontal periodic component, using current state-of-the-art differential processing.

6 Discussion and conclusions

We have shown significant periodic biases between a primary and secondary RF when using trajectory models that do not match the observed station coordinate behavior. This periodic bias further emphasizes the need to use and publish periodic components in all future versions of the ITRF, to enable more consistent secondary realizations. We have shown periodic motion amplitude inconsistencies of up to 2 mm for the SIRGAS-CON RF and up to 4 mm in SP07b (applying internal constraints), which can be reduced to well below 1 mm using state-of-the-art differential solutions. In this work, we have discussed the signals that can be modeled using four Fourier coefficients, and how an RRF can inherit behavior manifest in the parent’s frequency space. Of course, the Fourier coefficients do not account for interannual variations in GNSS periodic motions, which remains a topic for future research (Zou et al. 2014).

We have extended the inheritance process to work within frequency space. We note that inheritance can also extend to other parameters affecting RFs, such as ‘jump’ or ‘relaxation’ space, characterizing step functions and logarithmic transients that occur after a large earthquake. If all the stations in a local network were strongly affected by co- and postseismic deformation, these relevant trajectory parameters (i.e., the coefficients describing the amplitudes of the jump and the transient) are impossible to constrain without adding stations that are outside of the earthquake’s influence area. However, a local RF can inherit the jump or relaxation behavior from a parent GRF or RRF, thereby avoiding the need to expand the aperture of the network.

So far, inheritance can be invoked in the following parameter spaces: the frequency space (amplitude of sine and cosine terms for each frequency), the velocity space (the station velocity components), and the position space (scale, origin, and orientation). Emphasizing the use of ‘partial’ rather than ‘full’ inheritance, being equivalent to a partial alignment of two frames, we recognize that this notation is not in everyday use. As we previously described, not all state space parameters need to be inherited from a primary frame, and not all parameters being inherited need to come from a fixed set of stations. Thus, it is possible to select the best fitting stations for each space or parameter type. For example, for VBCA, one might use this station during velocity space inheritance but not for periodic space because of noisy data in the parent frame. On the other hand, another station might be useful for periodic space inheritance but not for velocity space.

Using long-term stack inheritance rather than traditional densification via H transformations of individual solutions presents another significant advantage. Stack inheritance allows modern RF realizations where the trajectory models’ parameters are not static after the frame is realized (also known as a dynamic or kinematic reference frame, KRF). After a new release of the ITRF, mapping agencies may be required to maintain consistency with legacy versions of the ITRF. They can use inheritance to initialize their realizations as closely as possible to the desired version of ITRF and then use a kinematic approach to maintain their RF stack. Although KRFs are outside of this work’s scope, P07b was conceived as a KRF and used inheritance to achieve a high consistency with ITRF 2014.