Accuracy of the International Terrestrial Reference Frame origin and Earth expansion

1. Introduction

[2] The coordinate origin of the ITRF is theoretically defined as the mean CM of the total Earth system about which all Earth satellites orbit. Since other geodetic techniques are either insensitive to or inaccurate about the location of CM, the ITRF2008 origin is realized by averaging weekly instantaneous frame origins spanning 26 years of SLR observations. The scale of ITRF2008, on the other hand, is realized by averaging scales provided by VLBI, spanning 29 years of observations, and SLR [Altamimi et al., 2011]. Since their definitions can only be accessed through inexact measurements and dynamic models, realized origin and scale implied by coordinates and velocities of an ITRF solution will be slightly different from the true mean CM and metric scale. In other words, the coordinates and velocities of the ITRF network will contain net translational and dilatational errors. Also, Helmert transformations between different ITRF realizations will result in origin and scale offsets and drifts. Hence, such net errors are often referred to as offset and drift errors in the realized origin and scale (the net translational error in the geocentric coordinates of the network is equivalent to an error in the realized origin with a negative sign). Drifts of ITRF origins at the level of 1 mm yr⁻¹ away from the mean CM caused by measurement errors and deficiencies in the data analyses have been suspected from successive ITRF realizations [Altamimi et al., 2007]. Such drifts indicated by successive realizations, however, only reflect internal inconsistency of the ITRF. A more independent assessment of the ITRF origin stability is desirable but difficult to achieve at this level of accuracy.

[3] An origin drift will contaminate estimates of any process that is determined from the velocities of geodetic stations that are given with respect to that reference frame. The stability of the ITRF origin has thus become a major concern for measurements of sea level rise, present-day surface mass trend (PDMT), glacial isostatic adjustment (GIA), and geocenter velocity (between CM and the geometric center-of-figure of the solid Earth surface (CF)), where 1 mm yr⁻¹ is important compared to the signals. For example, the resulting orbit drift has significant effects on direct altimetric measurements of regional sea levels including the non-uniform spatial pattern (fingerprint) of sea level change due to melt and runoff water [Mitrovica et al., 2001]. The effect on global mean sea level is smaller but non-negligible [Morel and Willis, 2005]. Moreover, altimeters may contain instrument drifts themselves that are validated and monitored by global tide-gauge measurements [Mitchum, 1998]. The validation procedures may use GPS derived velocities to convert relative sea levels to the geocentric frame and so will be affected by possible scale rate error and origin drift because of the hemispherical station distribution bias [Collilieux and Woppelmann, 2011].

[4] Various theories and hypotheses have been proposed about the expansion of the solid Earth and debated throughout the plate tectonic revolution [e.g., Carey, 1976]. Even though there is no proven theory or hypothesis of an expanding Earth, legitimate concerns and speculations persist about a possible mean rate of change in the radius of the solid Earth on accretionary, geothermal, climate change, cosmological, magmatic [Mjelde et al., 2010], and other physical grounds, despite an indication to the contrary from records of historical moment of inertia [Williams, 2000]. Although attempts have been made to actually measure the expansion using space geodetic techniques [e.g., Heki et al., 1989], the confounding effects of limited spatial distribution of geodetic sites and the presence of other major geophysical processes have limited the accuracy of these attempts.

[5] Recently, a somewhat independent estimation of geocenter velocity has been obtained using relative surface velocities with respect to the center-of-network (CN) and other data [Wu et al., 2010]. The global simultaneous inverse approach also provides a suitable platform for accurate parameter estimation with much reduced contaminations. Here, we extend that study to invert multiple geodetic data sets of 3-dimensional absolute ITRF2008 [Altamimi et al., 2011] surface velocities, gravity trends measured by the Gravity Recovery and Climate Experiment (GRACE) mission [Tapley et al., 2004], and linear trends from two ocean bottom pressure (OBP) models. ITRF2008 origin drift components and a mean solid Earth expansion rate are estimated and resolved simultaneously with rigid plate motions, PDMT and GIA from the data combination and a dynamically constructed a priori GIA model with a full covariance matrix.

2. A New Global Simultaneous Estimation

[6] The most recently realized ITRF2008 geocentric surface velocities (with respect to CM) at suitably located sites (see below) r_i observed by VLBI, SLR or GPS can be written in the center-of-mass of the solid Earth (CE) coordinate system [Farrell, 1972; Wu et al., 2010] as:

where the first three terms on the right-hand side are the components of an origin drift vector d from the mean CM to the ITRF2008 origin. _j is the unit coordinate vector along a Cartesian or spherical polar axis. The sum of the next three terms is the site velocity in the CE frame _i caused by radial expansion, plate motion, PDMT loading and GIA. is the constant mean radial expansion rate of the solid Earth. As described by the supplementary notes of Wu et al. [2010], ω_p is the rotation vector of the pth plate, and _nmq, _nmq^v,h, and _nmq^v,l are normalized spherical harmonic coefficients of present-day surface mass density trend, vertical and horizontal GIA velocities respectively. Y_nmq are normalized spherical harmonic functions of colatitude ϑ and longitude ϕ. Note that the last term _cm is a function of _1mq only [Blewitt et al., 2001; Wu et al., 2010]. h′_n and l′_n are fixed load Love numbers.

[7] A simple estimation of geocenter velocity (or radial expansion) used in the geodetic community is based on the average velocity (or the average radial velocity) of the tracking network. As station velocities are given with respect to CM, geocenter velocity is the negative velocity of CF in this frame, which is often approximated by the negative velocity of CN. Equivalently, the same estimates can be obtained by using reduced observation equations with equal weighting: V_i = , or V_i = _r, where and are the net translational velocity and radial expansion rate of the network respectively. For example, the estimate of geocenter velocity is − = − ΣV_i/N = −Σ_i /N + _cm + d, where N is the number of sites. If the site distribution were sufficiently global and dense, then Σ_i/N = _cn ≈ _cf, and − would be the estimate of geocenter velocity biased only by the origin drift. The plate motion, PDMT and GIA processes contain no degree-0 or net expansion term, and would not contaminate the average estimate of radial expansion with the ideal site distribution. However, in reality, the geodetic sites are sparse, limited largely to easily accessible land areas, and have a strong distribution bias toward the northern hemisphere and certain tectonic plates. Consequently, _cn ≠ _cf, and the simple estimates of geocenter velocity (or origin drift if geocenter velocity is known) and radial expansion are significantly contaminated by the geophysical processes that do not average out for the limited network.

[8] To accurately determine the origin drift and radial expansion, we combine absolute ITRF2008 velocities with GRACE gravity and OBP models to estimate all parameters (except the assumed known and fixed Earth properties) on the right-hand side of equation (1), and spherical harmonic coefficients of GIA-induced gravitational geoid trend simultaneously. The spherical harmonic series are truncated at degree and order 60. With equation (1) and detailed observation equations for GRACE and OBP described by Wu et al. [2010], the three data types can be symbolically described as Y in the following vector equation:

where A through G are coefficient matrices. , ^v,l, and ^v,h are parameter vectors containing spherical harmonic coefficients up to n = 60 for PDMT, horizontal and vertical GIA velocities. ^v,k is the parameter vector of spherical harmonic coefficients up to degree 7 for GIA geoid trend. Higher-degree GIA geoid trend coefficients are approximately related to GIA vertical coefficients and are substituted by them [Wu et al., 2010]. ω contains the rotation vectors of 15 major tectonic plates. These parameter vectors are estimated simultaneously by the method of least squares with reduced a priori information [Wu et al., 2006].

[9] ITRF2008 geocentric velocities from a global network of 233 SLR, VLBI and GPS sites distributed as evenly as possible are selected and supplemented by the published and transformed (to ITRF2008) continuous and episodic GPS velocities at 448 sites in North America, Fennoscandia, Alaska, Antarctica and Greenland [Wu et al., 2010 and references therein] (see Figure 1 for site distribution). They are not located in orogenic or local tectonic areas and are at least 200 km away from plate boundaries. Sites with suspected sediment loading or man-made ground water extractions are also excluded. Only variances of velocities are used with empirical colored noise components [Dixon et al., 2000] since covariances are largely not available. Systematic errors and model deficiencies including those in non-gravitational force models, atmospheric delays, station biases and weighting, and source/antenna geometries tend to be the dominant source for origin and scale drift errors but are not properly represented in the covariance matrix. These require independent assessments such as introducing the origin drift error vector in our estimation, and comparing scales provided by SLR and VLBI. The same linear trends in the spherical harmonic coefficients of the GRACE gravitational geoid derived by Wu et al. [2010] are modeled and used here without modification. A full covariance matrix is constructed from a calibrated (using a simple scaling factor) formal covariance matrix, and an empirical covariance matrix based on temporal fit residuals. We compare and use the linear trends from Jet Propulsion Laboratory (JPL)'s ECCO data-assimilating OBP model [Fukumori et al., 1999] and the Ocean Model for Circulation and Tides (OMCT) [Dobslaw and Thomas, 2007] averaged on a 3° × 3° grid between 73°S to 73°N (Figure 1). OMCT is a global simulative ocean circulation model being used as a de-aliasing product for the GRACE data analysis. Assuming equal uncertainties for the two models, a baseline uncertainty for each model is derived by dividing the RMS model differences at all grid points by , and used as the model uncertainty at grid points where model differences are less than or equal to the RMS difference. At grid points where model differences are larger than the RMS difference, model differences divided by are used as the model uncertainties. While such assessed diagonal covariance matrix reduces weighting from sampling points with large model differences indicative of potential model defects, spatially coherent errors in the models are not well represented. We will conduct inversions using ECCO and OMCT separately. The difference in the results should include a measure of the effect of the coherent errors.

[10] The kinematic vector observation equation (2) contains various rank deficiencies and underdeterminedness. A key ingredient for a successful inverse resolution is the use of a dynamically constructed a priori GIA model (see Wu et al. [2010] for details). The a priori mean values of all GIA parameters are predicted by standard ice history and Earth rheology models ICE-5G/IJ05/VM2 [Peltier, 2004; Ivins and James, 2005] using a spectral method [Vermeersen and Sabadini, 1997]. Deglaciation history over confirmed or suspected regions described in the literature and mantle viscosities are perturbed by large amounts and propagated to construct a full covariance matrix. Although conservatively assessed, such a priori GIA information contributes significantly to the inversion by retaining intrinsic relations among various GIA signatures in the data combination and ensuring dynamic consistency of the estimated GIA parameters.

3. Results

[11] Table 1 lists origin drift and Earth expansion results from the new and expanded inversions. The magnitudes of origin drifts along the coordinate axes are fairly small at the level of 0.3 mm yr⁻¹ or less when ECCO is used for data combination (Table 1, third column). Replacing ECCO by OMCT results in a larger-magnitude drift component of −0.6 mm yr⁻¹ along the Z-axis without noticeably changing the drift values along the X- and Y-axes (Table 1, fourth column). For the drift components, the formal uncertainties from the inversions are all around 0.1 mm yr⁻¹. The listed total uncertainties in Table 1 include additional components that reflect the differences in results by the two OBP models with the effect of spatially coherent model errors. The drift results using the data-assimilating ECCO model are particularly encouraging since they are statistically not very far from zero considering the uncertainties. We adopt the inversion using a simple average of ECCO and OMCT in the data combination as our main case study with the results listed in Table 1 (fifth column).

Table 1. Estimated Global Parameters and Uncertainties^a

Component	Symbol	ITRF2008+GRACE+ECCO	ITRF2008+GRACE+OMCT	ITRF2008+GRACE+ECCO+OMCT
Origin Drift	dx(mm yr−1)	−0.3 ± 0.1	−0.4 ± 0.1	−0.4 ± 0.1
	dy(mm yr−1)	−0.2 ± 0.1	−0.1 ± 0.1	−0.2 ± 0.1
	dz(mm yr−1)	−0.3 ± 0.2	−0.6 ± 0.2	−0.5 ± 0.2
Expansion	(mm yr−1)	0.1 ± 0.2	0.1 ± 0.2	0.1 ± 0.2

• a The uncertainties listed include components derived from differences in results from the two OBP models and a systematic scale drift (for the Earth expansion rate), see text for discussion. The formal inversion uncertainties for the origin drift components and radial expansion are all about 0.1 mm yr⁻¹.

[12] The combination of relative ITRF velocities, GRACE, and an OBP model with the a priori GIA model information resolves geocenter velocity contributions from both PDMT and from GIA. These allow the new inversion of the combination of absolute ITRF velocities, GRACE and OBP models to accurately determine the origin drift simultaneously with parameters of the geophysical processes. The contribution of relative ITRF velocities to the geocenter velocity determination comes from degree-1 deformation [Blewitt et al., 2001; Wu et al., 2010], which is different from the translational motion [Lavallée et al., 2006] and unlikely to be affected in the same way by the systematic errors that cause the origin drift. Therefore, we argue that the inversion provides a valuable quasi-independent validation of the ITRF origin accuracy.

[13] Different OBP models in the data combination have no material effect on the mean expansion rate of the solid Earth. The formal uncertainties for the expansion are all smaller than 0.1 mm yr⁻¹. The estimates, however, may be subject to an additional systematic drift error in the ITRF2008 scale, which is realized through the average scale information from VLBI and SLR. As discussed before, this systematic error is not properly represented in the covariance matrix. Currently, the uncertainty for the scale drift is estimated to be 0.025 × 10⁻⁹ yr⁻¹ by comparing information from the two independent techniques [Altamimi et al., 2011]. The total uncertainties for the expansion rate listed in Table 1 are thus about 0.2 mm yr⁻¹. In either case, the tiny estimated expansion rates of about 0.1 mm yr⁻¹ are not statistically different from zero.

4. Discussion and Conclusions

[14] The origin accuracy of the ITRF is critically important in monitoring global change processes such as sea level rise, ice mass imbalance, and GIA. By exploiting a recently demonstrated inverse sensitivity of the combination of relative surface velocities, GRACE and the ECCO OBP model to the total geocenter velocity [Wu et al., 2010], we have found a quasi-independent way to quantitatively evaluate the origin drift in an ITRF solution. The global simultaneous inversion of the data combination also provides an effective platform to mitigate the confounding effect of a sparse ground geodetic network and the presence of the multiple geophysical processes. While the technique depends on the OBP models to some extent, we conclude that the ITRF2008 origin realization is consistent with the mean CM at a level of 0.5 mm yr⁻¹. This complements the results of an internal comparison of the two most recent ITRF realizations [Altamimi et al., 2011]. Currently, the ITRF realization assumes a linear motion model. Some geodetic velocities span longer periods than the GRACE and OBP data. Therefore, a significant portion of the determined drift may result from interannual variations and inhomogeneous temporal data distribution.

[15] A possible secular change in the mean radius of the solid Earth has long been suspected from various scientific disciplines. The problem was not settled due to lack of a proven theory and the inherent difficulty in actually determining such a change. By using the data combination and the simultaneous global inverse approach, the issue of aliasing from the geophysical processes due to network sparseness is largely overcome here. In fact, our uncertainty in the mean expansion rate determination is dominated by an ITRF scale drift uncertainty of 0.16 mm yr⁻¹ estimated from comparing VLBI and SLR. In conclusion, no statistically significant present expansion rate is detected by our study within the current measurement uncertainty of 0.2 mm yr⁻¹.