The expanding Earth at present: evidence from temporal gravity field and space-geodetic data

The Earth expansion problem has attracted great interest, and the present study demonstrates that the Earth has been expanding, at least over the recent several decades. Space-geodetic data recorded at stations distributed globally were used (including global positioning system data, very-longbaseline interferometry, satellite laser ranging stations, and stations for Doppler orbitography and radiopositioning integrated by satellite), which covered a period of more than 10 years in the International Terrestrial Reference Frame 2008. A triangular network covering the surface of the Earth was thus constructed based on the spherical Delaunay approach, and average-weighted vertical variations in the Earth surface were estimated. Calculations show that the Earth is expanding at present at a rate of 0.24 ± 0.04 mm/yr. Furthermore, based on the Earth Gravitational Model 2008 and the secular variation rates of the second-degree coefficients estimated by satellite laser ranging and Earth mean-pole data, the principal inertia moments of the Earth (A, B, C) and in particular their temporal variations, were determined: the simple mean value of the three principal inertia moments (i.e., [A+B+C]/3) is gradually increasing. This clearly demonstrates that the Earth has been expanding, at least over the recent decades, and the data show that the Earth is expanding at a rate ranging from 0.17 ± 0.02 mm/yr to 0.21 ± 0.02 mm/yr, which coincides with the space geodetic evidence. Hence, based on both space geodetic observations and gravimetric data, we conclude that the Earth has been expanding at a rate of about 0.2 mm/yr over recent decades.


Introduction
Whether the Earth is expanding or contracting at present is an interesting problem in science.Some scientists support the viewpoint of Earth expansion, and some are against this viewpoint.
The expanding Earth hypothesis is mainly supported by paleontology, paleomagnetism, paleoclimatology and geology data [e.g., Scalera 2003a].Wilson [1960] declared that the Earth is expanding based on geological evidence and the Wegener's continental drift hypothesis.By simulating the expansion process of the Earth, Creer [1965] concluded that the Earth radius R E was around 0.55R in the early Precambrian (ca.3,800 Myr ago), around 0.94R to 0.96R in the early Paleozoic (ca.544 Myr ago), and around 0.96R to 0.97R in the early Mesozoic (ca. 230 Myr ago), where R = 6371 km, the mean radius of the Earth at present.Based on various geological evidence, Dearnley [1965] showed that R E was around 4400 km before 270 Myr ago, and around 6000 km before 6.5 Myr ago, and that the mean radius has increased at a rate of about 0.6 mm/yr.Carey [1976Carey [ , 1988] ] concluded that the Earth is expanding within the ocean-floor expansion framework, and Chen [2000] reported that the Earth started to expand around 4300 Myr ago, with a radius of about 4651 km and at an increased rate of around 0.4 mm/yr at that time, compared with the radius increase at a rate of around 0.1 mm/yr at present.A study by Müller [2010, http://zeitexpansion/de] suggested that not only is the universe expanding, but also the Earth, and the radius of the Earth increases at a rate of about 0.6 mm/yr.
Following a series of studies of three paleogeographical reconstructions for the Paleocene, Cretaceous and Jurassic, Scalera [1998Scalera [ , 2001] ] carried out a further reconstruction for the Triassic period with the assistance of paleomagnetic data.He found that the data are best reconciled if they are treated in an expanding Earth framework, and the mean rate of increase of the Earth radius is 15 mm/yr.Based on analyses of the Laser Geodynamics Satellites and very-long baseline interferometry (VLBI) data for stable nonorogenic continental regions, Gerasimenko (1993) indicated an increase in the Earth radius of 4.15 ± 0.27 mm/yr.However, later studies of Gerasimenko [1996Gerasimenko [ , 1997] ] showed that the rate of increase of the Earth radius does not exceed 1 mm/yr, by analyzing the VLBI data.Further, Gerasimenko [2003] reprocessed the NASA VLBI data (NASA GSFC solution number 1122, June 1999) and obtained a possible radius increase of 0.2 mm/yr after removing the network points that influenced their heights by as much as ±4.0 mm/yr.

Subject classification:
Station coordinates and vertical velocities, Temporal gravity variation, Principal inertia moments, Earth expansion.
It appears that the increase in the Earth radius by about 15 mm/yr declared by Scalera [2001] is too large.Later, Scalera [2003aScalera [ ,b, 2006a] ] admitted the existence of a much slower expansion rate at present, such as a few millimeters or fractions of millimeters per year, based on the following reasoning.One of the main indications for a very low expansion rate at present is the global map of the spreading rates in geologic time [see McElhinny and McFadden 2000, Fig. 5.19: the half spreading map of the ocean floors], which shows a minimum ocean floor expansion at present.The other indication is the neotectonic period [Ollier and Pain 2001], starting from a few millions of years ago, which is in strong contradiction with the subductive processes that should be active for hundreds of millions of years.Scalera [2008] proposed a new orogenic model for more efficient extrusion of deep material onto the surface to constitute fold belts exactly in the periods of the minimum of this global expansion.This model has no large-scale subduction, but the phase changes towards a more unpacked lattice of the isostatically rising masses.Thus, this does not contradict the observed neotectonic period.Scalera [2003aScalera [ , 2006a] ] also suggested that in the expanding Earth frame, the polar motion parameters can be prolonged in the current geological past, up to 100 Myr ago, which reproduces with satisfactory approximation the more reliable true polar wander path and its slowing down around 50 Myr ago, which was considered as being detected by Besse andCourtillot [1991, 2002].Then, Scalera [2003aScalera [ ,b, 2006a] ] concluded that the expansion contributes partly to the observed secular decrease in the gravity zonal harmonics J 2 and the increasing length of day.However, the expansion induced =-5.53 ×10 -13 /yr does not coincide with the observed value = =-(2.60±0.3)×10 -1 /yr [e.g., Groten 2004], and the expansion-induced increase in length of day is not quite certain, since it is related to not only the variation in the Earth inertia moments that are also under argument, but also the tidal friction.
On the other hand, some scientists favor the viewpoint that during the geological time, the Earth radius has held invariant, or has been subject to very tiny changes that cannot be detected with the present observation accuracy and by present techniques.The major argument in favor of a constant radius is the belief that a substantial increase in size of the Earth should have had an observable effect on its rotational dynamics (mainly the length of day), which was considered as not being confirmed by the paleontologic data [Burša and Šidlichovský 1984].It has been demonstrated that the observed decrease in the second zonal harmonic of the geopotential field and that of the angular velocity of the Earth rotation affect the Earth rotation dynamics significantly.Burša [1985] further argued that the angular momentum equation of the Earth-moon system is well satisfied using recent data, and that there is no space for an expanding Earth hypothesis as well as for secular variations in the Newtonian gravitational constant [Dirac 1937].In the study of Burša and Hovorková [1994], which argues against the conclusion of Carey [1988], three different kinds of Earth density models were assumed to estimate the internal energy necessary for expansion from 0.6R to R, and they stated that no dynamic evidence exists for the origin of this energy over the last 450 Myr.Different from their opinions, we here remark on two points.First, the angular momentum of the Earth-moon system can still be held invariant if the expansion does exist.For instance, a small expansion of the Earth gives rise to a small slow-down of the Earth rotation, which causes a small expansion of the moon orbit around the Earth.Various observations by lunar laser ranging (LLR) show that the moon orbit is expanding, the origin of which is mainly due to the tidal friction, and might be partly due to the Earth expansion that causes the deceleration of the Earth rotation to some extent.Second, the expansion energy source might also include contribution from nonuniform expansion in the interior of the Earth.For instance, the loss of the potential energy arising from denser particles falling down towards the center of the Earth could compensate for the energy necessary for expanding the outer part of the Earth.
We note that some scientists favor the viewpoint that the Earth is contracting within a shorter period.The main evidences come from seismic and space-geodetic studies.Using seismic data, Lyttleton [1965Lyttleton [ , 1983] ] concluded that the radius of an original all-solid Earth would have been 370 km greater than the present radius.Based on the coordinates and velocities, and on their error estimations of the global positioning system (GPS), satellite laser ranging (SLR) and VLBI stations issued by the International Earth Rotation Service, Huang et al. [2002] and Sun et al. [2006] concluded that the Earth is contracting.Carey [1976] suggested that the Earth expansion hypothesis can be identified using the observed gravity variations.His suggestion, however, has not been well considered as the gravity data at that time were not accurate enough to detect relatively fine volume variations of the Earth.Now, precise gravimetric measurements have been developed, especially superconducting gravimetry and satellite gravimetry (e.g., the Challenging Mini-Satellite Payload (CHAMP), the Gravity Recovery and Climate Experiment (GRACE), and the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) sattelites), and these enable us to access temporal information in gravity fields with relatively high accuracy.Then, based on the recently released gravity field model, i.e., the Earth Gravitational Model 2008 (EGM2008) [Pavlis et al. 2008], the Earth (temporal) principal inertia moments (PIMs), A, B and C, can be determined, as they are closely related to the second-degree harmonic coefficients of the gravity field Heiskanen andMoritz 1967, Marchenko andAbrikosov 2001].Here, we note that EGM2008 did not provide new variation rates of the second-degree harmonic coefficients, and the values of the second-degree harmonic coefficients released by EGM2008 hold almost the same as those of EGM96, and consequently, for the purpose of the present study, using either EGM96 or EGM2008 does not make any difference.In the present study, the EGM2008 gravity field model is used.
The variations in A, B and C just reflect the variations in the Earth's figure and mass distribution with respect to time.As is well known, these PIMs (i.e., A, B and C) will become larger on average (the mean of the PIMs is just [A+B+C]/3 in this report) if the Earth expands.On the other hand, the rapid development of space-geodetic techniques (e.g., GPS, VLBI, SLR, Doppler orbitography and radiopositioning integrated by satellite (DORIS)) makes it possible to directly estimate the variations in the Earth volume (radius) over recent decades.In the present study, we demonstrate that the Earth has been expanding, at least over the recent decades, based on both the temporal PIMs of the Earth and the coordinates and vertical velocities of the globally distributed space-geodetic observation stations.
In summary, the topic of the Earth expansion is a relatively controversial subject.In the present study, we just focus on the variations in the Earth radius over the recent several decades (~30 yr).Here, it should be noted that in terms of the Earth expansion problem, there are still a lot of publications that are not cited in this paper, and interested readers can find them in the relevant publications.
The Earth expansion problem is related to different research branches and the determination of the expansion rate might have a significant influence in science.In Section 2, the Earth expansion rate is estimated based on the spacegeodetic data of 629 space-geodetic stations under the International Terrestrial Reference Frame 2008 (ITRF2008) (including GPS, SLR, VLBI and DORIS stations), covering a period longer than 10 yr.Then, in section 3, the temporal PIMs are determined based on the EGM2008 gravity field model, as well as the secular variation rates of the seconddegree harmonic coefficients that are provided by SLR observations [e.g., Cheng et al. 2011] and secular variations in the mean pole of the Earth.Furthermore, the Earth expansion rate is estimated here.In the last section, as well as our concluding remarks, we discuss some problems related to this Earth expansion and the mass growth of the Earth.

Evidence from space geodetic data
Investigations into the secular vertical displacements of the Earth crust are effective ways to judge whether the Earth is expanding.In this section, we estimate the Earth expansion (contraction) rate based on the space geodetic data covering more than 10 yr under ITRF2008.

Data and selection of stations
The spatial geodetic data used in this study are the ITRF solutions, which are reference frames that consist of sets of station positions and velocities with their variance/covariance matrices given.There are different versions of ITRFs [Altamimi et al. 2007[Altamimi et al. , 2011]].The most recently released version of ITRF is ITRF2008 [Altamimi et al. 2011, see also http://itrf.ensg.ign.fr/ITRF_solutions/2008/]. In this framework, there are 1572 recordings at present from various stations (including GPS, VLBI, SLR and DORIS).Due to discontinuities of many of these recordings, only 935 stations were available for the present study.In another aspect, stations that are very close to each other (e.g., the distance between any two stations is smaller than m; see Section 2.3) might be considered as one, and consequently they are merged into one station by weighted average in the present study.For instance, the stations with code numbers '21736S003' and '21736S005' are taken as one station.The same handling is applied to the stations with station codes '7810','7824' and 'KELY'.Thus, there are 845 stations left in all.
To calculate the expanding (or contracting) rate of the Earth, the most straightforward method is to estimate the vertical (or radial) velocities of every point on the surface of the Earth, and then to estimate the Earth expansion rate based on the averaged vertical velocities of all of the observation sites (stations), with their corresponding weights that are determined according to the corresponding observation accuracies.We constructed a spherical Delaunay triangular irregular network (TIN) from the stations provided by ITRF2008.The spherical Delaunay TIN (SDTIN) consists of spherical triangles that just cover the whole surface of the Earth, and the vertical velocity of a triangle can be expressed by a representative vertical velocity (see below).Then the expansion rates of the Earth can be estimated by the representative vertical velocities of all of the triangles that comprise the SDTIN.
The selection of the stations should be taken with great care.First, the distribution of the ITRF stations is not uniform.The number of stations in the northern hemisphere is much greater than that in the southern hemisphere.This means that not all of the stations have the same weight in calculating the Earth expansion rates.Second, the stations located in active tectonic zones (e.g., orogen belts or zones) should be removed from our calculations.
To find out which stations are located in active tectonic zones, we first examined the distributions of the orogen belts (figure 1).In the present study, we used the information of the plate boundaries and orogen belts based on the Peter Bird 2002 (PB2002) plate model [Bird 2003; see Figure 1 Bird [2003], the orogen belts 10 3 THE EXPANDING EARTH AT PRESENT are defined as the zones where the horizontal velocity at any point cannot be modeled from Euler poles.Then, we can reasonably assume that when the observation stations are located on orogen belts, the vertical velocities at these stations are not reliable.Hence, for our calculations, we removed the stations that are located in the orogen belts.
As shown in Figure 2, the absolute values of the vertical velocities of some stations are greater than 0.02 m/yr, and so large vertical movements of these kinds of stations are not related to the Earth expansion, but to tectonic activities like earthquakes and volcanoes, for example.Hence, these kinds of stations were not included in our calculations.After removing these stations that are located in orogen zones and those with vertical velocities greater than 0.02 m/yr, there are 629 stations left, and these are the stations used in our calculations.

Methods
To estimate the Earth expansion rate based on the vertical velocities of the chosen stations, we can use the TIN, which is constructed from data points (stations) spread over the surface of the Earth.The TIN is formulated using the Delaunay approach [e.g., Chew 1989, Park et al. 2001, Nico et al. 2005].The Delaunay TIN is a set that consists of triangular figures that satisfy the following conditions: (1) they are connected side-by-side to each other, and any two arbitrary figures do not overlap; and (2) the network just covers the whole surface of the Earth.
In earlier studies [e.g., Sun et al. 2006, Shen andZhang 2008], the TIN constructed on a sphere was mapped on a plane, which is constructed in two steps.First, the stations on the surface of the Earth are projected onto two planes.Then the SDTIN is constructed by setting up the TIN on these two  planes.For convenience, this is referred to as the pseudo-SDTIN here.However, the topologic relation on a sphere is different from that on a plane, and the major defect of pseudo-SDTIN lies in its torsions of the topological relations near the equator area, and additional interpolation stations will not reduce this defect in most cases.Consequently, in the present study, we set up the TIN directly on a sphere [e.g., Renka 1997], which is referred to as the SDTIN.
To estimate the Earth expansion rate, we need to carry out the following procedures.First, select a group of stations, and project them onto the Earth surface (a spherical surface with the Earth mean radius R = 6371 km), and set up a TIN on this spherical surface.Second, the vertical velocity of a triangle as a whole is expressed by the representative vertical velocity of the triangle, and the latter is determined by those at the three endpoints of this spherical triangle (see below).Finally, the Earth expansion rate is estimated by averaging the representative vertical velocities of all of the triangles (which construct the SDTIN), with their weights taken into account.
To determine the representative vertical velocity of a triangle, for the present purpose, it is accurate enough to take a spherical triangle as a plane triangle K. Let us use D, E, F to denote the three endpoints of the triangle K, and o D , o E , o F as the vertical velocities at D, E, F, respectively.Then, with equal weights, the representative vertical velocity of the triangle is: (1) To evaluate the accuracy of o tr , the error propagation law can simply be applied to Equation (1).
Hence, the averaged vertical velocity of a triangle is determined by the velocities at its endpoints.Then, the expansion rate of the Earth can be estimated based on the weighted and averaged vertical velocities of all of the triangles that just cover the whole surface of the Earth.

Results
In our calculations, we apply formulae with weights (see Equation 2 in the following) to calculate the expansion rate of the Earth from the vertical velocities of the triangles set up by the SDTIN, which is built from the ITRF stations.
To set up a TIN, the algorithm requires that two stations should not be too close to each other.Indeed, besides the discontinuity problem mentioned above, some groups of stations are very close to each other because different observation instruments are used in a very small area for the purpose of their collocation ties.For convenience and without loss of the generality, as mentioned in Section 2.1, in our calculations we choose a distance threshold of m.If the distance of two stations is less than this threshold, the two stations are considered as one new station or representative station, the determination of the position and velocity of which are described in the following.Here we note that the choice of the threshold is quite arbitrary, and it does not significantly influence the result.The basic requirement is that it should not be too small (otherwise it is still not convenient and certain algorithms would fail), nor too large (otherwise it is not necessary).When there are k stations in a very small region, we consider these k stations as one representative station.The position and velocity of the representative station are obtained by averaging the corresponding positions and velocities of the k stations, with their weights considered.
When the representative vertical velocities of each triangle are obtained based on Equarion (1), we can estimate the expansion rate of the Earth using the following formula, expressed as: (2) where o G is the expansion rate of the Earth, o i the representative vertical velocity of the i-th triangle, N the number of triangles, S i the spherical area of the i-th triangle, v 2 oi the variance of the i-th triangle, and P i the weight of o i .We note that P i is proportional to the area S i of the i-th triangle, and inversely proportional to the variance v 2 oi of o i .
The 629 stations are taken into our calculations in our final estimate of the Earth expansion rate.The absolute value of the vertical velocity of each station included in this study is smaller than 0.02 m/yr.The numerical results show that the variation rate of the Earth radius is (0.54 ±0.05) mm/yr, as listed in Table 1.
Post glacial rebound (PGR) has long been thought of as a reasonable candidate in explaining the uplift of the surface of the Earth [e.g., Lemoine et al. 2006].However, based on the data released from the website http://grace.jpl.nasa.gov/data/pgr/ [Paulson et al. 2007], the estimated average uplift rate o u of the surface of the Earth caused by PGR is only o u = = (5.69×10-6 ±1.09×10 -4 ) mm/yr.This means that the global contribution of PGR to the Earth volume variation can be neglected.However, in the present study, the expansion rate is estimated from 629 stations that are distributed over the 10 3    globe.Hence, to obtain a more reasonable indication of the influence of PGR on the uplift of the surface of the Earth, it is necessary to calculate the uplift rate for each spacegeodetic station.In the present study, we expand the grid data of uplift rates of the lithosphere into spherical harmonic coefficients, and then the uplift rate from PGR can be calculated according to the following equation: (3) where  3, and the contributions of PGR to the vertical motions of the stations are shown in Figure 4.The distributions of the vertical velocities with PGR effects filtered out are shown in Figure 5.
After subtracting the PGR contributions from each space-geodetic station and following the procedures described above, the estimated expansion rate of the Earth with PGR effects considered is (0.24 ± 0.04) mm/yr, as listed in Table 1.
Hence, based on the space geodetic data, we can conclude that the Earth expansion rate has been about 0.24 mm/yr over recent decades.
As a remark, we note that based on 736 GPS stations, Shen and Zhang [2008] concluded that the Earth expansion rate is 0.5 mm/yr, where the stations with vertical velocities larger than 0.05 m/yr were removed, the PGR effects were not considered, and the stations located in orogen belts were included in the calculations.We note that from Table 1, if the PGR effects are not considered, the expansion rate would be (0.54 ±0.05) mm/yr, which coincides with that of Shen and Zhang [2008].
The accuracy listed in Table 1 is optimistic compared to the accuracy of the techniques mentioned above [Altamimi et al. 2008].However, we note that the estimated magnitudes of the accuracy in the present report are consistent with the accuracy of the velocities of some stations listed in ITRF2008 (e.g.station with code 7207), and hence the order of the accuracy is acceptable.

Evidence from temporal gravity fields
In this section, we further improve the study of Shen et al. [2008], who provided a preliminary result that showed that the Earth is expanding at a rate of about 0.3 mm/yr, based on temporal gravimetric observations.

Secular trends of the Earth inertia moments
The parameters related to the temporal gravity field used in this study are listed in   .
In the Earth's principal inertia axis system [i.e., the coordinate axes coincide with the principal inertia axes of the Earth; e.g., see Munk and MacDonald 1960], the PIMs of the Earth (i.e., A, B, C) are linked to the second-degree harmonic coefficients (A _ 20 , A _ 22 , etc., also expressed in the principal inertia axis system) of the gravity field by the following equation [Marchenko andAbrikosov 2001, Marchenko andSchwintzer 2003]: (5) where M and a are the Earth mass and mean equatorial radius, respectively, and the dynamic flattening H D is defined as [e.g., Heiskanen and Moritz 1967]: In Equation ( 5

_
2m are just those provided by the EGM2008 gravity model (we note that using either EGM96 or EGM2008, there is no significant difference for the purpose of the present study).
To estimate the mean changes in the PIMs (i.e., A, B and C), it is convenient to introduce a mean rotational inertia I (a simple mean of the PIMs): (10) Based on Equations ( 5) and ( 10), this immediately gets: (11) or: (12) Then, from Equations ( 5) and ( 12) one has: (13) and its differential with respect to time reads: (14) Indeed, generally, the Earth does not expand uniformly in all directions, and now we would like to discuss some details about this.
Based on Equation ( 5), the PIM variation rates are: (15) Then, A, B and C and their variation rates A : , B : and C : can be determined, once the second-degree harmonic coefficients, their derivatives, and H D and H : D are given.The second-degree harmonic coefficients, their derivatives, and H D are listed in Table 2.As references, the H D values given by some previous studies are listed in Table 3.Then, based on Tables 2 and 3, A, B and C are determined, which are listed in Table 4 along with some previous results from other studies.We note that to determine A : , B : and C : , H : D should be given a priori.However, H : D is hard to obtain due to the inaccurate calculations of the dynamical flattening itself.Nevertheless, it was roughly estimated as about -1×10 -10 /yr to -7.4×10 -10 /yr [e.g., Dehant andCapitaine 1997, Bourda andCapitaine 2004] from the secular variation rate of J 2 , which could be quite precisely determined by long-term SLR observations [Cheng and Tapley 2004].Indeed, according to Bourda and Capitaine [2004], we have the following equation: ; It should be noted that the C _ : 20 appeared in Bourda and Capitaine [2004], i.e., the quantity as given above appeared in their Equation ( 8), should be referred to the Earth's principal inertia axis system, i.e., it should be replaced by A _ : 20 in the present study.Otherwise the dynamic flattening of H D defined by Bourda and Capitaine [2004]    This study (e) -7.84±0.88(a) Roughly estimated from space-gravimetric observations.(b) Estimated from space-gravimetric data by neglecting second-order terms.
(c) Estimated from space-gravimetric data by neglecting second-order terms based on assumption that the trace of the inertia moment tensor in the principal inertial axis system is invariant.(d) Estimated from space-gravimetric data by neglecting second-order terms based on the assumption that the trace of the inertia moment tensor in the principal inertia axis system is invariant.(e) Estimated from space-gravimetric data by neglecting second-order terms without assuming the invariance of the trace of the inertia moment tensor in the principal inertia axis system.which also lists the corresponding results from some previous studies.Now, based on Tables 2-5, the variation rates of the PIMs (i.e., A, B, C) can be determined, and the results here are listed in Table 6.
The PIMs are critique indicators of the Earth dynamics.Variations in both the volume and mass of the Earth can alter the PIMs.Hence, the increase in the PIMs on average might be caused by several factors: 1) the expansion of the Earth; 2) the increase in the Earth mass; and 3) the rise of the sea level.
The influence of a sea level rise can be eliminated, because even if under the assumption that the global mean sea level rise rate is as large as around 2 mm/yr [it is actually not so large; see e.g., Cabanes et al. 2001, Church 2001, Milly et al. 2003, Cazenave and Nerem 2004, Holgate and Woodworth 2004, Church et al. 2004, Alley et al. 2005, Church and White 2006, Shepherd and Wingham 2007], we can show that the contribution of the sea level rise to the mean increases in the PIMs is too small compared with the observed (calculated) mean increases in the PIMs.The reasoning here is as follows.
Assume that the sea level rise rate is 2 mm/yr (a value that is known to be overestimated), and assume that the whole surface of the Earth is covered by 3-km-thick sea water before the sea level rise.After the sea level rise, the thickness of the sea water becomes 3 km + 2 mm.Then, the mean changes in the PIMs caused by this sea level rise can be expressed as: (18) where R, R 0 , R 1 are respectively the radii of the solid Earth, the whole Earth (solid and ocean), and the whole Earth after the sea level rise, and t 0 and t 1 are the densities of the water before and after the sea level rise, which can be expressed as: where in the calculations we have used the following values: R = 6368 km, R 0 = 6371km [Dziewonski and Anderson 1981], and R 1 = 6371000.002m.According to Equation ( 18), taking into account Equations ( 19) and ( 20), we get: (21) As I = (A+B+C)/3 is (based on the results calculated in the present report, see Tables 2 and 4): the relative variation rate of I caused by the sea level rise is: This value is much smaller than the observed result of 5.24×10 -11 /yr (see Table 6).Hence, the contribution of a sea level rise to the observed average increases in the PIMs can be neglected.We note that in the present study the observed result of I : is (see Tables 2 and 6): (24) At present, there is no definite evidence showing that the mass of the Earth is increasing.Hence, we can assume that the Earth mass M holds invariant and only the Earth's volume changes with time in our present investigations.Now, after eliminating the influences of sea level rise and mass increase (or growth), only the Earth expansion is left as a valid candidate to explain the mean increases in the PIMs.
As a remark, we note that there is continuous gravitational differentiation and geothermal convection .
Table 6.The variation rates of the Earth principal inertia moments (PIMs) and relevant parameters determined in the present study.The calculations needed from the parameters in the present report come from Table 2, Table 3 and Table 5.
within the Earth.The lighter elements ascend, while the heavier ones go down towards the central parts of the Earth.
If the Earth's volume holds invariant, this process should cause a decrease in I, which appears to conflict with our calculated results based on observations.However, if the Earth as a whole is expanding during this period, this increase in I might be properly explained.Not only Earth expansion and mass growth, but also large-scale mass redistributions can give rise to the variations of the inertia moment I, such as PGR, ocean circulation and ice melting.However, the PGR cannot explain the true polar wander over the long geological times [see Scalera 2006a for a detailed discussion, and references therein], and ocean circulation can only cause periodic, and not long period, variations in I. Polar ice melting and the associated sea level rise will increase C but reduce A and B, which contradicts the corresponding values determined in the present study.Hence, with our present knowledge, under the assumption that the Earth mass holds invariant (in the opposite case, see Section 4), the increasing tendency of I can only be explained by the Earth expansion hypothesis.
As a remark, in our investigations, different from Marchenko and Schwintzer [2003], we do not assume that the trace of the inertia tensor is invariant, i.e., in the present study, the assumption of the invariance of the trace of the inertia tensor is not applied.Indeed, such an assumption does not hold.As can be imagined, whether the Earth expands or not cannot be assumed a priori, but needs to be observed.Then, if the Earth does expand or contract, the relation A+B+C = constant will fail.
In the following, we estimate the Earth expansion rate based on the calculated increased rate of I.

A simple model to explain the secular trend of the Earth's inertia moment
First, we roughly estimate the expansion rate by assuming that the Earth is a uniform sphere with the mean radius a = =6371 km and a mean density t .5.54 g cm -3 at an initial time.Then, 1 yr later, the Earth becomes a uniform sphere with the radius a' and the density t due to the expansion (note that it is assumed that the total mass of the Earth holds invariant).The inertial moments before and after the expansion are I = 2/5Ma 2 and I' = 2/5Ma' 2 , respectively.Then: (25) where m is the expansion ratio per year: (26) From this equation, a' -a=0.167mm/yr is arrived at, the accuracy of which will be considered later.

A more realistic expansion model
In the following, we estimate the expansion rate of the Earth by adopting a spherically stratified Earth model, the preliminary reference Earth model (PREM) [Dziewonski and Anderson 1981].Although the Earth has lateral density variations at all depths, as magnitude estimates, taking the Earth as a spherically stratified body as given by PREM is sufficient.
According to PREM (see Table 7, where x = r/R is the normalized radius and R = 6371 km is the mean radius of the Earth), we divide the Earth into N spherical layers.Each layer is uniform and different layers can have different densities, based on the density distribution given by PREM.Then, the expansion of the Earth is equivalent to the total expanding effect of the N layers.In the present study, we treat every 500 m in depth as one layer if r is smaller than 6291 km, and the depths ranging from r = 6291 km to r = =6356 km are divided into 650 layers, with the thickness of each layer as 100 m.The rest of the crust and the ocean are considered as two layers.So, the Earth is divided into 13234 layers.
A uniform spherical layer has the constant density t, and let us denote the density and radii of the inner and outer surfaces of the i-th layer by t i , r i-1 and r i , respectively, where i = 1, 2, ..., N. At the initial time, the volume of the i-th layer can be expressed as: (27) while at time t (in the present study, the unit of time is 1 yr), the volume of the i-th layer is: where, as noted before, r i and r i-1 are the radii of the outer and inner surfaces of the i-th layer, respectively, and at the same time they are the radii of the inner surface of the (i+1)th layer and the outer surface of the (i-1)-th layer, respectively.
The radial expansion ratio i i (i = 1, 2, ..., N) can be defined as: (29) and based on Equations ( 27)-( 29) we get the following expression: (30) Assume that the mass of the i-th layer holds invariant during the expansion, then taking into account Equation (30), the density before expansion, t i (0), and that after expansion, t i (t), meet the relation: The inertia moments of the i-th layer before and after the expansion can be expressed as: (32) where i = 1, 2, ..., N. Then the inertia moments of the whole Earth before and after the expansion can be expressed as: (33) and the annual variation rate of the Earth inertia moments can be expressed as (set t = 1 year): (34) Or, alternately: Defining l i and m i as follows: (36) (37) we obtain: (38) Generally, each layer has different chemical compositions and physical phases, and consequently the expansion rate of each layer might be different in the case where the Earth does expand.However, it is relatively difficult to estimate the expansion rate of every layer as we do not know exactly how the expansion rate of every layer responds to the composition, pressure, temperature, stress, etc., of every layer.Nevertheless, as quantity estimation and without loss of the generality, in the present study we just considered three simplified models, stated as follows.
Model 1: the expansion ratio is a linear function of r i :      (39) Model 2: the expansion ratio is a constant for all of the layers: (40) Model 3: the expansion ratio is a power function of r i : (41) In Equations.( 39) to ( 41), the parameters a, b, c are constants to be determined, which correspond to Models 1, 2 and 3, respectively.No matter in which case (Model 1, 2 or 3), Equation ( 35) can be written in the following form: (42) Generally, the above equation has no analytical solution for variable h (h=a, b, c).We use the Monte Carlo method to search for a certain h to meet the following condition: (43) where the positive real number f can be arbitrarily small.Then, we find values that best fit Equation ( 43).Using these values, we can obtain the increasing rates i i =i(r i ) (i = 1, 2…, N) for every r i according to Equations ( 29) to (31).Then, the Earth expansion rate can be determined based on Equation (29).

Estimates of the Earth expansion rates and their accuracies
To evaluate the accuracy of the estimated expansion rate, the following relation needs to be established: (44) where l i and m i are defined by Equations ( 36) and (37).Then, application of the error propagation law to Equation (44) leads to the needed accuracy evaluation.
With different models proposed in Sections 3.2 and 3.3, the Earth expansion rates can be estimated accordingly, and the corresponding results are listed in Table 8.For reference, some, but not all, of the results from previous studies related to the Earth expansion (or contraction) problem are also listed in Table 8.We note that in our investigations, the results from different models suggest that the Earth expansion rates range from 0.17 ± 0.02 mm/yr to 0.21 ± 0.02 mm/yr.These conclusions coincide not only with our results based on the space geodetic data (see Section 2), but also with the results provided by many previous studies [e.g., Creer 1965, Dearnley 1965, Carey 1976, Chen 2000, Scalera 2003b].

Discussion and conclusions
The topic of the Earth expansion is a relatively controversial subject that developed mainly in the first decades of the twentieth century [e.g., see Carey 1976].Nowadays, some scientists support and some are against the conclusion that the Earth is expanding.Nevertheless, huge amounts of geodetic data (e.g., deformation, gravity, and more) justify a renewed effort to examine the observables that support a varying Earth radius.
The present study uses space-technique-observed secular deformations and secular gravity field variations (second degree harmonic coefficients) that show that the Earth is expanding at the rates of 0.24 ± 0.04 mm/yr (see Table 1) and 0.18 ± 0.02 mm/yr (see Table 8), respectively.
If the Dirac's large numbers hypothesis is accepted [Dirac 1937[Dirac , 1974]], which states that the gravitational constant G decays with the age of the universe, then it might be concluded that the Earth is expanding, as a body will become looser when G decreases.However, it is almost impossible to calculate the Earth expansion rate based on this Dirac's large numbers hypothesis, due to the complicated structure and the relatively poorly determined geophysical parameters, such as density, viscoelastic coefficients, Love numbers, et cetera.
Based on various geological evidences, many scientists have concluded that the Earth is expanding at a rate ranging from 0.1 mm/yr to several millimeters every year.Our investigations show that the Earth is expanding, at least at present, at a rate of about 0.2 mm/yr, and the evidences come from both temporal gravity fields and space geodetic data.Based on the EGM2008 temporal gravity field model, calculations show that the PIMs are gradually increasing, at a mean rate of 1.73 ×10 -11 /yr (i.e., I : ), which demonstrates the present expansion of the Earth at a rate of about 0.18 mm/yr.From another aspect, based on the coordinates and vertical velocities of 629 space geodetic stations that are spread over the globe, we find that the average radius of the Earth is increasing at a rate of about 0.24 mm/yr.Whether our conclusion can be extended to the geological time scale is a problem to be investigated further.One possible explanation for this mean increase in the PIMs might be due to the growth of the Earth mass.If we assume that the mean increase in the PIMs is due to the growth of the mass of the Earth [Clark 2008] rather than the expansion of the Earth, what can we conclude?In the case where the growth is uniform, the mass M of the Earth increases at a rate of about 5.24 ×10 -11 /yr based on the mean increases in the PIMs given in the present report.In other cases, the rate of increase of M depends on how the Earth's mass gain is distributed inside the Earth.
If the mass of the Earth is growing uniformly at present, the global superconducting gravimeters (SGs) [Goodkind 1999] should record a common signal of gravity increase of Dg, which depends on the distribution of the growth of the Earth mass.The relevant calculations are stated as follows.
Assume that the radius of the Earth does not change and that the increases in the PIMs come from the Earth mass growth DM, then the gravity at the surface of the Earth will have a long-term trend that can be expressed as follows: (45) The growth DM will be determined based on the observed quantity I : /I = 5.24×10 -11 under the assumption that this quantity is caused by the growth of the Earth mass.
For different assumptions of the distribution of the Earth mass growth, the growth rate of the Earth mass determined will be different.We discuss the following cases.
Case I.The mass growth DM is uniformly distributed in the whole of the Earth.
In this case, we have: Noting that M = (5.972± 0.007)×10 24 kg and R = 6371 km, we have: (47) Case II.The mass growth DM is uniformly distributed in the mantle and crust of the Earth.
In this case, the density growth can be expressed as: (48) where R CMB = 3480 km is the radius of the core-mantle boundary [Dziewonski and Anderson 1981].The relation between Dt and DI can be expressed as: where R mantle = 6346.6km is the radius of the mantle [Dziewonski and Anderson 1981].
As a short summary, we conclude that if the continuous increase in the PIMs is due to the growth of the Earth mass, then corresponding to the cases of the uniform distributions of the whole Earth (Case I), the mantle and crust (Case II), and the crust (Case III), the gravity on the ground will increase by 42.6 nGal/yr, 37.5 nGal/yr and 25.6 nGal/yr, respectively.
However, if the expansion is only related to the volume, the gravity change on the Earth surface can be expressed as: (54) Assume that the Earth expansion rate is dR/dt = 0.15 mm/yr, then the long-term trend of gravity on the ground is: (55) which means that the gravity on the ground will decrease by as much as -46 nGal/yr if the Earth is expanding at a rate of 0.15 mm/yr.
Due to their high accuracy and sensitivity, SGs can be used to provide confirmation as to whether the mass of the Earth is growing and/or whether the Earth is expanding.Unfortunately, as the accuracy of the secular term in SG recordings is of the level of 100 nGal [e.g., see Crossley and Hinderer 2008], it is relatively difficult to provide confirmation using only one SG record.However, it might be possible to provide confirmation using the datasets continually recorded by about 30 SGs distributed over the globe.
Indeed, here we can indicate that our space-gravimetric observations do not support the mass growth hypothesis, as in the present study it holds that A : the estimated average uplift rate o u of the surface of the Earth caused by PGR is only (5.69×10 -6 ± 1.09×10 -4 mm/yr.Subtracting the PGR effects from the vertical velocities of our 629 stations, the estimated expansion rate from these 629 stations is 0.24 mm/yr (see Table 1).Finally, if the Earth expands, the following question might be raised: How is it possible, if the Earth radius increases, to evaluate the GPS (or VLBI, or other techniques) baseline variations using the same radius as today for previous years?Is a systematic error introduced into the computations?An answer might be that GPS and VLBI are geometric Cartesian methods in which a pure square Cartesian reference frame is used without the use of the Earth radius.However, as the Earth is like a rotating platform and the position in space is determined by electromagnetic signal speeds and times (based on some kind of time-keeping system, e.g., atomic clocks), then the distances (positions) determined are related to the radius of the Earth.This problem remains to be solved.

Figure 1 .
Figure 1.Distribution of the 935 geodetic stations listed in ITRF2008.Green lines, plate boundaries; yellow areas, orogen belts.Information of both plate boundaries and orogen belts are from Bird [2003].Red, 306 geodetic stations located in orogen belts (* Oro); blue, other 629 stations (* Used).Stations: dots, GPS; circles, DORIS; squares, SLR; pentagons, VLBI.Green stations, those used in estimating the expansion rate in the present study.

Figure 2 .
Figure 2. Vertical velocity and distance from nearest plate boundaries or orogen boundaries of the 629 stations outside the orogen belts.Vertical velocities of most stations are less than 0.02 m/yr as absolute values, no matter how far the stations are from nearest plate boundaries or orogen belts.

Figure 3 .
Figure 3. Vertical velocities of the 629 geodetic stations outside the orogen belts.PGR effects are not filtered out.Large uplift motion can be seen in northern America and northern Europe.

Figure 4 .
Figure 4. Uplift due to PGR for the 629 stations outside the orogen belts.Large uplift trends can be seen in northern America and northern Europe.Moderate uplift trends can be seen in Antarctica.

Figure 5 .
Figure 5. Vertical velocities of the 629 stations outside the orogen belts with PGR effects appropriately filtered out.Moderate uplift trends can still be seen in northern America and northern Europe, indicating large uplift trends in these area cannot be simply attributed to PGR.The uplift trends in Antarctica are not as large as indicated by Figure 3.
,1,2) are the fully normalized seconddegree harmonic coefficients of the gravitational potential in the International Terrestrial Reference System.Here, C _ 2m and S can vary with C _ 20 because the latter is related to the coordinate systems; this contradicts the lack fact that H : D does not depend on the choice of the coordinate systems.Hence, in the present study, Equation (16) should be written as: from the precession rate of the equator.(b) Estimated by fitting the rigid Earth nutation series to the observed data.(c) Estimated by fitting the nutation series to observations from VLBI.(d) Estimated by least-square estimation of the values from all of the H D -values corresponding to the same precession rate as given in Marchenko and Schwintzer [2003], with reference epoch 1997.0.(e) Estimated by least-square estimation of the same values as given in Marchenko and Schwintzer [2003], with reference epoch 2000.0.
Case III.The mass growth is uniformly distributed in the crust of the Earth.Similar to case II, we can directly write out the following result: > B : > C : , while if the mass growth hypothesis holds, the relation A : = B : = C : should be satisfied approximately.As a remark, PGR might be an explanation for the global uplift of the surface of the Earth.Large vertical velocities in Greenland and Antarctica can be explained by PGR.However, based on the data released from the website http://grace.jpl.nasa.gov/data/pgr/[Paulson et al. 2007],

Table 2
Marchenko and Abrikosov [2001;inty are provided by the International Earth Rotation Service 2010(Petit and Luzum,  2010).The mean pole data (ftp://tai.bipm.org/iers/conv2010/chapter7/annual.pole)areused to derive the time series of C _ 22 can be calculated by an iterative approach, as given byMarchenko and Abrikosov [2001; see Equations 84-90 therein].A mean value of the dynamic flattening H D is determined based on six different H D values as listed in Table3, which is referred to the epoch J2000.0, and H D reads:

Table 2 .
Relevant [Marchenko and Abrikosov 2001, Marchenko and Schwintzer 2003, Groten 2004roten 2004, Lemonine et al. 1998] used in the calculations in the present study.The second-degree harmonic coefficients of potential are from the EGM2008 gravity field model.

Table 3 .
Values of the dynamic flattening H D .

Table 5 .
The estimated values of H :D .

Table 8 .
Expansion rates of the Earth.