**GNSS Observations of Crustal Deforamtion: A Case Study in East Asia**

### Shuanggen Jin

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/51536

### **1. Introduction**

The East Asia is located in a complex convergent region with several plates, e.g., Pacific, North American, Eurasian, and Philippine Sea plates. Subduction of the Philippine Sea and Pacific plates and expulsion of Eurasian plate with Indian plate collision [22, 36, 15, 18, 11] make the East Asia as one of the most active seismic and deforming regions (Fig. 1). Current deformation in East Asia is distributed over a broad area extending from the Tibet in the south to the Baikal Rift zone in the north and the Kuril-Japan trench in the east, with some ambiguous blocks, such as South China (Yangtze), Ordos and North China blocks, and pos‐ sibly the Amurian plate, embedded in the deforming zone. The inter-plate deformation and interaction between the blocks are very complicated an active, such as the 1978 Mw=7.8 Tangshan earthquake located between North China and Amurian blocks. Since the 1960s when the theory and models of seafloor spreading and plate tectonics were established, the large Eurasian plate has been considered as an independent rigid plate, such as RM2, P071 and NNR-NUVEL1A [6], and even the present-day global plate motion models of ITRF se‐ quences based on the space geodetic data [10]. In fact, East Asia is very complicated defor‐ mation zone with several possible rigid blocks [22]. [36] proposed that East Asia consists of several micro-plates. However, because of low seismicity and there being no clear geograph‐ ical boundary except for the Kuril-Japan trench and the Baikal rift zone, it is difficult to accu‐ rately determine the geometry and boundary of sub-plates in these areas. The borderlines of the micro-plates are not visible, especially in the complicated geological tectonic regions in North China far away from the Qinghai-Tibet plateau, and convergent belts of Eurasia, North America and the Pacific plate. Therefore, it is difficult to confirm the tectonic features and evolution of the deformation belts in East Asia.

© 2013 Jin; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Jin; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

In the northern part of East Asia, [36] first proposed the existence of the Amurian plate based on the clear geographical boundary of the Kuril-Japan trench and the Baikal rift zone, but it becomes diffuse throughout continental East Asia. The proposed Amurian plate (AM) in East Asia is of special interest to constrain the relative motion of the major and minor plate in East Asia and provides a rigorous framework for interpreting seismicity and the kinematics, especially for seismically active Japan. However, the southern boundary of the suggested Amurian plate is poorly understood. Figure 1 shows some recent proposed boun‐ daries. The upper blue solid line is from [4], the middle black dash boundary line is from [34], and the bottom red dash line is from [8].

**Figure 1.** Tectonics setting in East Asia. The un-continuous lines are the main fault lines and the dash lines are the undefined plate boundary. Circles are the earthquakes from Harvard CMT catalogue (1976-2005, Mw >5.0). The upper blue solid line is from [4], the middle black dash boundary line is from [34], the bottom red dash line is from [8], and the upper green solid+dash lines is this study.

Therefore, the existence of the Amurian plate and its boundary geometry remain controver‐ sial. Over the last two decades a number of investigations of the micro-plates tectonics in East Asia have been conducted using geologic, seismological and geodetic data [36, 34, 8, 27, 5, 2]. However all of these studies have suffered from limited data quality and quantity, re‐ sulting in ambiguous conclusions. For instance, [34] estimated the AM motion with earth‐ quake slip vectors and found that the spreading rates of the Baikal Rift are with the order of < 1mm/yr. [8] assumed Amuria and North China as an independent Amurian plate, and es‐ timated spreading rates at the Baikal Rift were about 10 mm/yr. [29] used only four GPS sites in or around the AM plate and concluded that AM can neither be resolved nor exclud‐ ed as a separate plate.[5] postulated that North China (including the possible AM) and South China could be a single rigid block using data from only 9 GPS sites. [2] claimed the existence of the Amurian plate, but still could not determine the location of the southern Amurian plate boundary. These investigations of the tectonics in this region are not conclu‐ sive because of the sparse and limited data that were used. Hence, there is much debate sur‐ rounding the nature of microplate and its boundaries in East Asia.

With more GPS observations in East Asia, it provides a new ways to clearly monitor the large-scale crustal motion and to distinguish the possible blocks, e.g., the national projects "Crustal Movement Observation Network of China (CMONC) with more than 1000 GPS sites and the Japanese GPS Earth Observation Network (GEONET) with more than 1000 continuous GPS sites. These dense GPS observations will obtain a more accurate estimate of plate geometry and its interior crustal deformation in East Asia. In this Chapter, we present new dense geodetic results for East Asia from about 1000 GPS sites in China, South Korea and Japan for the period 1998-2005, as well as combining recently published velocities for the Bailkal Rift and Mongolia [5]. The possibility of microplate motion independent of the Eurasian plate is tested using GPS derived velocities and its boundary and kinematics are further discussed [13].

In addition, accurate measurements of crustal strain accumulated energy rates will contrib‐ ute to understand tectonic features and to evaluate the earthquake potential. Now the high precision space geodesy techniques, especially the low-cost and all weather GPS, play a key role in monitoring the crustal strain state and accumulated energy variation. Meanwhile, the significant strain accumulation caused by the tectonic activities (such as earthquakes) will provide an essential constraint on the physical processing of earthquakes. Therefore, moni‐ toring the spatial variation of the strain and comprehensive understanding of strain accu‐ mulation pattern are beneficial to reveal the physical process of crustal tectonic activities and to evaluate the earthquake risk [14]. Here, the denser GPS velocity filed is used to esti‐ mate the strain rates and strain energy density rates, in an attempt to know largely aseismic areas and to assess the future earthquake risk potential in East Asia.

### **2. GPS observations and results**

In the northern part of East Asia, [36] first proposed the existence of the Amurian plate based on the clear geographical boundary of the Kuril-Japan trench and the Baikal rift zone, but it becomes diffuse throughout continental East Asia. The proposed Amurian plate (AM) in East Asia is of special interest to constrain the relative motion of the major and minor plate in East Asia and provides a rigorous framework for interpreting seismicity and the kinematics, especially for seismically active Japan. However, the southern boundary of the suggested Amurian plate is poorly understood. Figure 1 shows some recent proposed boun‐ daries. The upper blue solid line is from [4], the middle black dash boundary line is from

**Figure 1.** Tectonics setting in East Asia. The un-continuous lines are the main fault lines and the dash lines are the undefined plate boundary. Circles are the earthquakes from Harvard CMT catalogue (1976-2005, Mw >5.0). The upper blue solid line is from [4], the middle black dash boundary line is from [34], the bottom red dash line is from [8], and

Therefore, the existence of the Amurian plate and its boundary geometry remain controver‐ sial. Over the last two decades a number of investigations of the micro-plates tectonics in East Asia have been conducted using geologic, seismological and geodetic data [36, 34, 8, 27, 5, 2]. However all of these studies have suffered from limited data quality and quantity, re‐ sulting in ambiguous conclusions. For instance, [34] estimated the AM motion with earth‐ quake slip vectors and found that the spreading rates of the Baikal Rift are with the order of < 1mm/yr. [8] assumed Amuria and North China as an independent Amurian plate, and es‐ timated spreading rates at the Baikal Rift were about 10 mm/yr. [29] used only four GPS sites in or around the AM plate and concluded that AM can neither be resolved nor exclud‐

[34], and the bottom red dash line is from [8].

268 Geodetic Sciences - Observations, Modeling and Applications

the upper green solid+dash lines is this study.

The Crustal Movement Observation Network of China (CMONC) was constructed since 1998, which contains a nationwide fiducial network of 25 continuous GPS sites observed since August 1998, 56 survey mode sites with yearly occupations and more than 900 region‐ al campaign stations operated in 1999, 2001, and 2004 with continuous observations for at least 4 days during each session. In addition, the Korean GPS Network (KGN) with more than 45 permanent GPS sites was established since 2000 by the Korea Astronomy and Space Science Institute (KASI), the Ministry Of Governmental Administration and Home Affairs (MOGAHA), and the National Geographic Information Institute (NGI). Japan GPS Earth Observation Network (GEONET) was established since 1996 by the Geographical Survey In‐ stitute of Japan. These denser GPS observations can investigate and study detailed crustal deformation and kinematics in East Asia. Here we collect about 1000 GPS sites (1999.1- 2004.12) in East Asia with 54 core IGS sites that were used for ITRF2000 [1] and 10 perma‐ nent IGS sites located in East Asia. These GPS sites are shown in Figure 2. The pentagon stands for the permanent GPS sites (2000-2005), the dot denotes the campaign GPS sites and the triangle is the yearly observed GPS sites (1999-2005).

**Figure 2.** GPS sites distribution in this study. The triangle is the yearly observed GPS site, the pentagon is the continu‐ ous GPS site (2000-2005), and the dot is the campaign GPS site (1999-2005). The solid lines are the known plate boun‐ daries.

All available GPS data were processed in single-day solutions using the GAMIT software [16] in a three-step approach. At the first step, loose a-priori constraints were set for all pa‐ rameters and double-differenced GPS phase observations from each day were used to esti‐ mate station coordinates and the zenith tropospheric delay (ZTD) at each station every 2 hours. The IGS final orbits, IERS Earth orientation parameters, azimuth- and elevation-de‐ pendent antenna phase center models, as recommended by the IGS were used in the data processing. The 54 global IGS stations served as ties to the ITRF2000 frame [1]. At the second step, the regional daily solutions were combined with global solutions produced by the Scripps Orbital and Position Analysis Center (SOPAC, http://sopac.ucsd.edu/) using the GLOBK software [9], and the reference frame was applied to the solution by performing a seven-parameter transformation to align it to ITRF2000 (via the global 54 core stations). At the third step, the site velocities were estimated by least square linear fitting to time varia‐ tions of the daily coordinates for each station [13]. For example, Figure 3 shows the North, East and Vertical GPS time series at SUWN site (South Korea). The ~1000 GPS site velocity field combining with recently published velocities at the Baikal Rift and Mongolia [5] are shown in Figure 4 with respect to the Eurasian plate.

**Figure 3.** GPS time series at SUWN site (South Korea).

### **3. Sub-plates and kinematics**

#### **3.1. Block modeling method**

deformation and kinematics in East Asia. Here we collect about 1000 GPS sites (1999.1- 2004.12) in East Asia with 54 core IGS sites that were used for ITRF2000 [1] and 10 perma‐ nent IGS sites located in East Asia. These GPS sites are shown in Figure 2. The pentagon stands for the permanent GPS sites (2000-2005), the dot denotes the campaign GPS sites and

**Figure 2.** GPS sites distribution in this study. The triangle is the yearly observed GPS site, the pentagon is the continu‐ ous GPS site (2000-2005), and the dot is the campaign GPS site (1999-2005). The solid lines are the known plate boun‐

All available GPS data were processed in single-day solutions using the GAMIT software [16] in a three-step approach. At the first step, loose a-priori constraints were set for all pa‐ rameters and double-differenced GPS phase observations from each day were used to esti‐ mate station coordinates and the zenith tropospheric delay (ZTD) at each station every 2 hours. The IGS final orbits, IERS Earth orientation parameters, azimuth- and elevation-de‐ pendent antenna phase center models, as recommended by the IGS were used in the data processing. The 54 global IGS stations served as ties to the ITRF2000 frame [1]. At the second step, the regional daily solutions were combined with global solutions produced by the Scripps Orbital and Position Analysis Center (SOPAC, http://sopac.ucsd.edu/) using the GLOBK software [9], and the reference frame was applied to the solution by performing a seven-parameter transformation to align it to ITRF2000 (via the global 54 core stations). At the third step, the site velocities were estimated by least square linear fitting to time varia‐ tions of the daily coordinates for each station [13]. For example, Figure 3 shows the North, East and Vertical GPS time series at SUWN site (South Korea). The ~1000 GPS site velocity field combining with recently published velocities at the Baikal Rift and Mongolia [5] are

the triangle is the yearly observed GPS sites (1999-2005).

270 Geodetic Sciences - Observations, Modeling and Applications

shown in Figure 4 with respect to the Eurasian plate.

daries.

The definition of the micro-plate geometries in East Asia is quite unclear, especially their boundaries. Denser space geodetic data can be used to define the plate boundary by testing an independent plate rotation about a best-fit Euler's rotation pole obtained by geodetic ve‐ locities. Here we assumed several plates or blocks in NE Asia whose plate boundary geome‐ tries were defined principally on the basis of seismicity and faults: EU (Eurasia), AM (Amuria), SK (South Korea), NC (North China), SC (South China), AK(AM+SK), AN (AK +NC), EA (East Asia) (see Figure 1). The boundaries are respectively the Yin Shan-Yan Shan Mountain belts for AM-NC, the Qinling-Dabie fault for NC-SC and the Tanlu fault for NC-SK blocks. As the blocks have large areas, the plates are just defined as rigid blocks motions without considering the elastic strain [32], i.e.

$$
\Delta \mathbf{v} = \boldsymbol{\Omega} \times \mathbf{r} \tag{1}
$$

where *v* is the velocity and *r* is the position vector. In order to calculate Euler's rotation pa‐ rameters, we project the fault geometry and stations positions from spherical to planar ge‐ ometry [20]. Then the rigid block motion are modeled to estimate the Euler ration parameters of each block (see Table 1). Here only continuously and yearly observed GPS sta‐ tions are used. In addition, the angular velocity of the Eurasian plate was estimated from the velocities of 22 IGS sites located on the Eurasian plate (TROM, MADR, HERS, BRUS, KOSG, 7203, ZIMM, VILL, OBER, ONSA, WETT, POTS, GOPE, GRAS, BOR1, LAMA, KIRU, JOZE, ZWEN, IRKT, KIT3, KSTU), including the core sites in the Eurasian plate employed for the orientation and maintenance of the ITRF2000 [1] (see Table 1).

**Figure 4.** Crustal deformation rates in the Eurasia plate (EU) fixed reference frame with error ellipses in 95% confi‐ dence limits.

#### **3.2. Testing results**

To test whether the microblocks are independent of the Eurasian plate, we used a *χ* <sup>2</sup> test that compares how well two different models fit a set of data. *χ* <sup>2</sup> is a sum of the squares of weighted residuals defined as:

$$\chi^2 = \sum\_{l=1}^{N} \frac{\left(\nu\_o(l) - \nu\_w(l)\right)^2}{\sigma\_o^2(l)}\tag{2}$$

where *vo*(*i*) is the observation velocity of site *i*, *vm*(*i*) is the calculated velocity of site *i* from the plate rotation model, *σ<sup>o</sup>* 2 (*i*) is the variance of observation velocity in site *i*, and *N* is the total number of observations. Table 2 shows the *χ* <sup>2</sup> for each model. The *χ* <sup>2</sup> for the model of independent AM, SK, SC, NC blocks is smaller than the Eurasian plate and a 2-block model (EU and EA), respectively. To check whether AM, SK, SC and NC are independent blocks, we perform an F-ratio statistical test.

*v r* =W´ (1)

where *v* is the velocity and *r* is the position vector. In order to calculate Euler's rotation pa‐ rameters, we project the fault geometry and stations positions from spherical to planar ge‐ ometry [20]. Then the rigid block motion are modeled to estimate the Euler ration parameters of each block (see Table 1). Here only continuously and yearly observed GPS sta‐ tions are used. In addition, the angular velocity of the Eurasian plate was estimated from the velocities of 22 IGS sites located on the Eurasian plate (TROM, MADR, HERS, BRUS, KOSG, 7203, ZIMM, VILL, OBER, ONSA, WETT, POTS, GOPE, GRAS, BOR1, LAMA, KIRU, JOZE, ZWEN, IRKT, KIT3, KSTU), including the core sites in the Eurasian plate employed for the

**Figure 4.** Crustal deformation rates in the Eurasia plate (EU) fixed reference frame with error ellipses in 95% confi‐

To test whether the microblocks are independent of the Eurasian plate, we used a *χ* <sup>2</sup> test that compares how well two different models fit a set of data. *χ* <sup>2</sup> is a sum of the squares of

2

( ( ) ( )) ( )

*vi vi i*

*o m*

s

2

c

1

=

*i o*

*N*

2


orientation and maintenance of the ITRF2000 [1] (see Table 1).

272 Geodetic Sciences - Observations, Modeling and Applications

dence limits.

**3.2. Testing results**

weighted residuals defined as:

First, we assume the one-block model in which East Asia (EA) is part of the Eurasian plate and the two-block model in which the East Asia (EA) plate rotates independently with re‐ spect to the Eurasian plate (EU). The 3-block model assumes that the EA is divided into the AN (Amuria+South Korea+North China) and South China (SC) plates, while the 4-block model contains the EU, SC (S.China), NC (N.China) and AK (Amuria+S.Korea) plates. The 5 block model is the EU, SC, NC, AM and SK plates. These blocks were defined principally on the basis of seismicity and faults (see Figure 1). We compare the misfit of each model inver‐ sion and test for significance using the F-ratio [30]:

$$F = \frac{[\mathbb{X}^2(1\ black) - \mathbb{X}^2(2\ black)]/3}{\mathbb{X}^2(2\ black)/(N-3)}\tag{3}$$

where *χ* <sup>2</sup> (1*block*) and *χ* <sup>2</sup> (2*block*) stand for the sum of the squares of weighted residuals in one-block and two-block models, respectively. In Table 2 one can clearly see that for the two-block model (EU and EA) the reduced chi-squared misfit of GPS velocity observations has been greatly reduced from 6.4 (for the one-block model (EU+EA)) to 1.4 (for the twoblock model (EU and EA)) and the F-ratio statistic (Eq.2) is 1086.5, which is well above the 99% confidence level of 3.8. The reduced chi-squared misfits of other independently rotating blocks for the 3-block, 4-block and 5-block models are also greatly degraded, and the calcu‐ lated F-statistics between the 2-3 block and 3-4 block models are well above the 99% confi‐ dence level of 3.8 as well as between the 3-4 block and 4-5 block models. These results indicate that the AM, SK, NC and NC are independent of the Eurasian plate motion. Fur‐ thermore, it shows that the South Korea block (SK) is excluded from the Amurian plate (AM), coinciding with recent test results using fewer GPS sites [12]. Figure 5 shows residual velocities (observed minus predicted) at yearly observed and continuous GPS site in East Asia with respective to the 1 block (Eurasia), 2 blocks (Eurasia and East Asia) and 4 blocks (Eurasia, Amuria+S.Korea, North China and South China). The mean residual in 4-block model is much smaller than fewer blocks.

The estimated relative motions along the block boundaries are further obtained (Fig. 6). Comparisons of spreading or converging rates and directions along these boundaries are showing in Figures 7 and 8. The GPS-derived relative motion directions are nearly the same as the earthquake slip vector directions (Fig. 6). However there are some differences at Bai‐ kal Rift. It may be due to the larger uncertainty of earthquake slip vectors ( ±15 ). Another is possibly the Euler vector problem of the large non-rigid Eurasian plate (EU). For instance, Table 1 lists different Euler vectors of the Eurasian plate, and larger discrepancy is found between [5] and other geodetic results. The tectonic boundaries between the North China and Amuria plates, the Yin Shan-Yan Shan Mountain belts, are extending at about 2.4 mm/yr. The Qinling-Dabie fault between the North China and South China plates is moving left laterally at about 3.1 mm/yr. This difference between Amuria and South China predicted rates is about 5.5 mm/yr, almost consistent with geological results by [3] and [24]. The Tanlu fault between the North China and South Korea blocks is moving right laterally at about 3.8 mm/yr. The Amuria and South Korea blocks are extending at about 1.8 mm/yr. The conver‐ gent rates at the boundaries of the AM and Okhotsk [2] are from 9 to 17 mm/yr, similar to the seismic results of [17]. The spreading rates in the Baikal Rift zone are about 3.0 ± 1.0 mm/yr, consistent with [5] at 4 ± 1 mm/yr.

**Figure 5.** Residual velocities (observed minus predicted) at yearly observed and continuous GPS site in East Asia with respective to the 1 block (Eurasia), 2 blocks (Eurasia and East Asia) and 4 blocks (Eurasia, Amuria+S.Korea, North China and South China).



between [5] and other geodetic results. The tectonic boundaries between the North China and Amuria plates, the Yin Shan-Yan Shan Mountain belts, are extending at about 2.4 mm/yr. The Qinling-Dabie fault between the North China and South China plates is moving left laterally at about 3.1 mm/yr. This difference between Amuria and South China predicted rates is about 5.5 mm/yr, almost consistent with geological results by [3] and [24]. The Tanlu fault between the North China and South Korea blocks is moving right laterally at about 3.8 mm/yr. The Amuria and South Korea blocks are extending at about 1.8 mm/yr. The conver‐ gent rates at the boundaries of the AM and Okhotsk [2] are from 9 to 17 mm/yr, similar to the seismic results of [17]. The spreading rates in the Baikal Rift zone are about 3.0 ± 1.0

**Figure 5.** Residual velocities (observed minus predicted) at yearly observed and continuous GPS site in East Asia with respective to the 1 block (Eurasia), 2 blocks (Eurasia and East Asia) and 4 blocks (Eurasia, Amuria+S.Korea, North China

**(∘ /My)**

**Eurasian plate**

This study -100.655 56.995 0.257± 0.002 0.6 0.1 49

Altamimi et al.[2002] -99.374 57.965 0.260± 0.005 - - - Calais et al. [2003] -107.022 52.266 0.245± 0.005 - - -

**Pole σ***maj*

**Error σmin**

**Ellipse Azimuth**

**Plates Longitude Latitude Angular rate**

NNR-1Aa -112.3 50.6 0.234

mm/yr, consistent with [5] at 4 ± 1 mm/yr.

274 Geodetic Sciences - Observations, Modeling and Applications

and South China).

**Table 1.** Rotation is in a clockwise direction about the pole. The error ellipses of the poles are described by the 1 σ semi-major and semi-minor axes of each error ellipse and the clockwise angle from true north of the semi-major axis. a No-Net-Rotation NUVEL-1A (NNR-1A) model [6] b Angular velocity vectors were estimated from 12 years of CGPS in ITRF2000 [25]Absolute and relative angular velocity vectors for the Eurasian, Amurian, South Korea, North China and South China plates.


**Table 2.** <sup>a</sup> 1: EU (Eurasia); 2: EU and EA (East Asia); 3: EU, AN (Amuria+South Korea+North China) and SC (South China); 4: EU, AK (Amuria+S.Korea), NC, SC; 5: EU, NC, SC, AM and SK; 6: EU, NC, SC, SK, West AM and East AM. *f* is the number of degrees of freedom and χ*<sup>f</sup>* 2 is the reduced χ <sup>2</sup> as the ratio of χ <sup>2</sup> to *f* .Statistic tests of different block models

#### **4. Crustal Strain rates and Seismic Risks**

Monitoring the pattern of crustal strain and comprehensive understanding of strain accu‐ mulation intensity are beneficial to reveal the physical process of crustal tectonic activities and to evaluate the earthquake risk. As the first step in the earthquake risk potential evalua‐ tion in East Asia, the strain parameters are estimated from the estimated GPS displacement rate field. In order to reduce the effects of abnormal site motions, the subnetwork with four GPS sites is used to estimate the strain parameters. Under the hypothesis that the velocity field *v* varies linearly inside each small sub-network covering the GPS sites, we can calculate the average horizontal velocity gradient *g* = *grad*(*v*) over each subnetwork. Because the ve‐ locity gradient generally incorporates both deformation and rotation, this 2-D tensor is asymmetric [21]. The crustal strain rate in East Asia can be derived from GPS deformation velocities by [12, 14]:

$$\begin{aligned} \boldsymbol{\nu}\_{el} &= \frac{\widehat{\boldsymbol{\mathcal{O}}} \boldsymbol{\nu}\_{el}}{\widehat{\boldsymbol{\mathcal{O}}} \boldsymbol{\mathcal{X}}\_{el}} \boldsymbol{\mathcal{X}}\_{el} + \frac{\widehat{\boldsymbol{\mathcal{O}}} \boldsymbol{\mathcal{X}}\_{el}}{\widehat{\boldsymbol{\mathcal{O}}} \boldsymbol{\mathcal{X}}\_{nl}} \boldsymbol{\mathcal{X}}\_{nl} \\ \boldsymbol{\nu}\_{nl} &= \frac{\widehat{\boldsymbol{\mathcal{O}}} \boldsymbol{\nu}\_{el}}{\widehat{\boldsymbol{\mathcal{O}}} \boldsymbol{\mathcal{X}}\_{el}} \boldsymbol{\mathcal{X}}\_{el} + \frac{\widehat{\boldsymbol{\mathcal{O}}} \boldsymbol{\mathcal{V}}\_{nl}}{\widehat{\boldsymbol{\mathcal{O}}} \boldsymbol{\mathcal{X}}\_{nl}} \boldsymbol{\mathcal{X}}\_{nl} \end{aligned} \tag{4}$$

**Figure 6.** Relative spreading or converging motions at the plate boundaries in East Asia.

**Figure 7.** Relative spreading or converging rates at the plate boundaries in East Asia.

**4. Crustal Strain rates and Seismic Risks**

276 Geodetic Sciences - Observations, Modeling and Applications

velocities by [12, 14]:

Monitoring the pattern of crustal strain and comprehensive understanding of strain accu‐ mulation intensity are beneficial to reveal the physical process of crustal tectonic activities and to evaluate the earthquake risk. As the first step in the earthquake risk potential evalua‐ tion in East Asia, the strain parameters are estimated from the estimated GPS displacement rate field. In order to reduce the effects of abnormal site motions, the subnetwork with four GPS sites is used to estimate the strain parameters. Under the hypothesis that the velocity field *v* varies linearly inside each small sub-network covering the GPS sites, we can calculate the average horizontal velocity gradient *g* = *grad*(*v*) over each subnetwork. Because the ve‐ locity gradient generally incorporates both deformation and rotation, this 2-D tensor is asymmetric [21]. The crustal strain rate in East Asia can be derived from GPS deformation

> *ei ei ei ei ni ei ni ni ni ni ei ni ei ni*

(4)

*v v vxx x x v v v xx x x*

¶ ¶ = + ¶ ¶ ¶ ¶ = + ¶ ¶

**Figure 6.** Relative spreading or converging motions at the plate boundaries in East Asia.

where *vei* and *vni* are the east and north component velocity at the site *i* located at *( xei* , *xni* ). Strain components *ε*˙*ee*, *ε*˙ *nn* and *ε*˙ en are expressed as ∂*ve* ∂ *xe* , ∂*vn* ∂ *xn* and <sup>1</sup> 2 ( ∂*ve* ∂ *xn* + ∂*vn* ∂ *xe* ) , respec‐ tively. The dilation rates i show that East Asia is under the compressional strain regime at WNW-ENE, consistent with the focal mechanism of earthquakes in Northeast Asia. The high dilation rates appear in North China, Southwest Japan and at the boundary of Philip‐ pine Sea plate. The strong compression rates are probably caused by the extrusion force due to the subduction of the Philippine Sea and Pacific plates and the expulsion of Eurasian plate with Indian plate collision, causing frequent earthquakes in these regions. Inversely, the South Korea and South China blocks have relatively lower dilation rates, indicating a lower indirect effect of push and subduction forces or as if such forces are transmitted through South Korean peninsula and South China without causing any deformation/strain, alternatively. This may be attributed to the strong rheology or/and absence of relatively weak zones in the region.

In addition, we estimate the scalar strain rate, defined as

$$
\dot{\varepsilon} = \sqrt{\dot{\varepsilon}\_{ev}^{2} + \dot{\varepsilon}\_{m}^{2} + 2\dot{\varepsilon}\_{en}^{2}}\tag{5}
$$

where *e* and *n* are longitude and latitude directions, respectively. Figure 9 shows the contour map of scalar strain rates in Northeast Asia, implying the Korean peninsula and South Chi‐ na as stable blocks with low strain rates. It once again highlights that high strain rates con‐ centrate in North China, Southwest Japan and the boundary of Philippine Sea plate with Eurasian plate, consistent with high seismicity in these areas (Figure 1).

**Figure 8.** Relative spreading or converging directions at the plate boundaries in East Asia.

In addition, the accumulated strain energy is generally released through earthquakes until the adjacent fault blocks or plates reach a new state of equilibrium [26, 33]. Therefore, the release of tectonic strain energy stored within crustal rock is the cause of major earthquakes. The strain energy per unit volume (i.e. the strain energy density) is an important index re‐ flecting the intensity of crustal activities, and its variation rate indicates the long-term trend of accumulated energy within the crust. The larger the variation rate of strain energy densi‐ ty, the higher energy accumulated in the crust, which would more probably result in earth‐ quakes. Therefore, for the earthquake risk evaluation and prediction, it is important to estimate the strain energy density from surface displacement observations and determine the state of strain energy density within the crust and its temporal variations.

**Figure 9.** Map of scalar strain rates in East Asia from GPS observations.

where *e* and *n* are longitude and latitude directions, respectively. Figure 9 shows the contour map of scalar strain rates in Northeast Asia, implying the Korean peninsula and South Chi‐ na as stable blocks with low strain rates. It once again highlights that high strain rates con‐ centrate in North China, Southwest Japan and the boundary of Philippine Sea plate with

S.China-Philippine Amuria-Philippine

Latitude (degree)

In addition, the accumulated strain energy is generally released through earthquakes until the adjacent fault blocks or plates reach a new state of equilibrium [26, 33]. Therefore, the release of tectonic strain energy stored within crustal rock is the cause of major earthquakes. The strain energy per unit volume (i.e. the strain energy density) is an important index re‐ flecting the intensity of crustal activities, and its variation rate indicates the long-term trend of accumulated energy within the crust. The larger the variation rate of strain energy densi‐ ty, the higher energy accumulated in the crust, which would more probably result in earth‐ quakes. Therefore, for the earthquake risk evaluation and prediction, it is important to estimate the strain energy density from surface displacement observations and determine

Amuria-Okhotsk Baikal Rift (Amuria-Eurasia)

30 31 32 33 34

This study

Kreemer et al

Sella et al.

Seno et al.

Sella et al. This study

Wei and Seno

Prawirodirdjo and Bock

Heki et al.

Wei and seno

52 54 56 58

Eurasian plate, consistent with high seismicity in these areas (Figure 1).

15 20 25 30 35

Sella et al.

Nuvel-1 This study Wei and seno

40 45 50 55 60

**Figure 8.** Relative spreading or converging directions at the plate boundaries in East Asia.

the state of strain energy density within the crust and its temporal variations.

110

40

60

80

100

Slip direction (degree)

120

130

Seno et al. This study Sella et al.

140

150

278 Geodetic Sciences - Observations, Modeling and Applications

For an elastic body, the strain energy equals the work done by external forces and its density is the strain energy per unit volume. The general tensor form for the strain energy density can be expressed in terms of strain and stress using Hooker's Law:

$$U = \frac{1}{2}\sigma\_y \varepsilon\_y\tag{6}$$

where *U* is the strain energy density (Unit: *J*.*m*−<sup>3</sup> ), *σij* and *εij* are the stress and strain, re‐ spectively. And the variation rate of strain energy density can be further derived from Eq. (6). The stress is obtained through the laws of elasticity theory as follows [31]:

$$
\Delta \sigma\_y = 2\,\mu\varepsilon\_y + \delta\_y \lambda\Delta \tag{7}
$$

where *μ* is the modulus of rigidity, *λ* is the Lame parameter, *δij* is Kroneckerdelta, *Δ* is the 2- D surface dilation (∑ *i*=1 2 *εii* ) . For Poisson's ratio *v* =0.25, *λ* = *μ* , and it is assumed as standard value of 3×1010 Pa [7].

The stress ( *σij* ), strain ( *εij* ) and their rates can be derived from GPS displacements (1999-2004) and velocities, respectively. Using Eq. (6), the strain energy density variation rate in East Asia can be obtained using the derived the strain, stress and their rates, which is shown in Figure 10. The distribution of strain energy density variation rates shows that the most active areas are in North China, Southwest Japan and west margin of Philippine Sea plate, respectively, again consistent with high seismic activity zones. As the GPS measure‐ ments are made after the large historic earthquakes, the strain energy density rates derived from GPS displacement rate may include contributions from postseismic relaxation. These regions with anomalous large strain energy density rates probably indicate a high earth‐ quake risk in future, and the lower strain energy density rates in the South Korean peninsu‐ la and South China imply that low seismicity may continue in the future.

**Figure 10.** Variation rate of crustal strain energy density in East Asia

### **5. Conclusions**

GPS data (1998-2005) from more than 85 continuous and about 1000 campaign stations in East Asia have been processed. The kinematics of East Asia is studied by modeling GPS-de‐ rived velocities with rigid block rotations and elastic deformation. It has been found that the deformation in East Asia can be well described by a number of rotating blocks, which are independent of the Eurasian plate motion with statistical significance above the 99% confi‐ dence level. The tectonic boundary between the North China and Amuria plates is the Yin Shan-Yan Shan Mountain belts with about 2.4 mm/yr extension. The boundary between North China and South China is the Qinling-Dabie fault, moving left laterally at about 3.1 mm/yr. The Amuria and South Korea blocks are extending at about 1.8 mm/yr. The Baikal Rift between the Amurian and Eurasian plates is spreading at about 3.0 mm/yr. The 9~17 mm/yr relative motion between the Amuria and Okhotsk blocks is accommodated at the East Sea-Japan trench zone. Furthermore, the relative motion rates and deformation types are nearly consistent with seismic and geological solutions along their boundaries. In addi‐ tion, the AM, SK and SC blocks are almost rigid with residual velocities on order of 1.0~1.2 mm/yr, while the NC block has larger residual velocities on order of 1.6 mm/yr, indicating un-modeled deformation in block boundaries. Localized deformation near the Qinling-Da‐ bie fault and Yin Shan-Yan Shan Mountain belts may be elastic strain accumulation due to interseismic locking of faults.

The strain and energy density rates in East Asia are investigated with GPS observations (1999.1-2004.12). The dilation rates show that East Asia is under the compressional strain re‐ gime at WNW-ENE, consistent with the focal mechanism of earthquakes in East Asia. The high dilation rates focus on North China, Southwest Japan and the boundary of Philippine Sea plate, probably caused by the compression force due to the subduction of the Philippine Sea and Pacific plates and the expulsion of Eurasian plate with Indian plate collision. In con‐ trast, the South Korean Peninsula and South China blocks havbe relatively lower dilation rates, indicating a possible lower effect of push and subduction forces or that such forces are transmitted through South Korean peninsula and South China without causing any defor‐ mation/strain. This may be attribute to the strong rheology or/and absence of weak zones in the region, which leads to fewer earthquakes. Moreover, the scalar strain rates and strain en‐ ergy density rates further imply the Korean peninsula and South China as a stable block with low rates, and high rates mainly concentrate on North China and Southwest Japan and the western boundary of Philippine Sea plate, consistent with highly seismic occurrences in these areas. In addition, the strain energy density rate reflects a long-term trend of strain en‐ ergy accumulation and release. Therefore, North China, Southwest Japan and western boun‐ dary of Philippine Sea plate with high strain energy density rates are still highly seismic and the low seismicity in South Korea and South China with lower strain energy density rates may continue in the future.

### **Acknowledgements**

most active areas are in North China, Southwest Japan and west margin of Philippine Sea plate, respectively, again consistent with high seismic activity zones. As the GPS measure‐ ments are made after the large historic earthquakes, the strain energy density rates derived from GPS displacement rate may include contributions from postseismic relaxation. These regions with anomalous large strain energy density rates probably indicate a high earth‐ quake risk in future, and the lower strain energy density rates in the South Korean peninsu‐

110 115 120 125 130 135

Longitude(°)

72.5

106

GPS data (1998-2005) from more than 85 continuous and about 1000 campaign stations in East Asia have been processed. The kinematics of East Asia is studied by modeling GPS-de‐ rived velocities with rigid block rotations and elastic deformation. It has been found that the deformation in East Asia can be well described by a number of rotating blocks, which are independent of the Eurasian plate motion with statistical significance above the 99% confi‐ dence level. The tectonic boundary between the North China and Amuria plates is the Yin Shan-Yan Shan Mountain belts with about 2.4 mm/yr extension. The boundary between North China and South China is the Qinling-Dabie fault, moving left laterally at about 3.1 mm/yr. The Amuria and South Korea blocks are extending at about 1.8 mm/yr. The Baikal Rift between the Amurian and Eurasian plates is spreading at about 3.0 mm/yr. The 9~17 mm/yr relative motion between the Amuria and Okhotsk blocks is accommodated at the East Sea-Japan trench zone. Furthermore, the relative motion rates and deformation types

5.71

39.139.1

5.71

5.71

39.1

5.71

39.1 39.1

72.5 72.5

106

139173 206 240

5.71

39.1

0

50

100

150

200

250 mJ.m-3/yr

26

**5. Conclusions**

28

30

32

34

Latitude(°)

5.71

5.71

**Figure 10.** Variation rate of crustal strain energy density in East Asia

5.71 5.71

5.71

5.71 5.715.71

39.1

72.5

106

5.71

5.71

5.71

5.71

5.71

Custal strain energy density rate

36

38

40

42

280 Geodetic Sciences - Observations, Modeling and Applications

la and South China imply that low seismicity may continue in the future.

Figures were made with the public domain software GMT [Wessel and Smith, 1998]. We are grateful to those who created the Crustal Motion Observation Network of China and made the observation data available, Korean GPS Network and Japan GPS Earth Observation Net‐ work (GEONET).

### **Author details**

Shuanggen Jin

Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai, China

### **References**


[15] Kato, T., Kaotake, Y., & Nakao, S. (1998). Initial results from WING, the continuous GPS network in the western Pacific area. *Geophysical Research Letters*, 125(3), 369-372.

**References**

107(A10), 2214, 10.1029/2001JB000561.

282 Geodetic Sciences - Observations, Modeling and Applications

Asia, . *Geophys. Res. Lett*, 20, 895-898.

*syst*, 4(3), 1027, 10.1029/2001GC000252.

*Res*, 108(B10), 2501, 10.1029/2002JB002373.

*Int*, 101, 425-478.

*phys. Res*, 104, 29147-29155.

2348-2350.

[1] Altamimi, Z., Sillard, P., & Boucher, C. (2002). ITRF2000: A New Release of the Inter‐ national Terrestrial Reference Frame for Earth Science Applications. *J. Geophys. Res*,

[2] Apel, E. V., Burgmann, R., Steblov, G., Vasilenko, N., King, R., & Prytkov, A. (2006). Independent active microplate tectonics of northeast Asia from GPS velocities and

[3] Avouac, J. P., & Tapponnier, P. (1993). Kinematic model of deformation in central

[4] Bird, P. (2003). An updated digital model of plate boundaries. *Geochem. Geophys. Geo‐*

[5] Calais, E., Vergnolle, M., San'kov, V., Lukhnev, A., Miroshnitchenko, A., Amarjargal, S., & Déverchère, J. (2003). GPS measurements of crustal deformation in the Baikal-Mongolia area (1994-2002): Implications for current kinematics of Asia. *J. Geophys.*

[6] De Mets, C., Argus, D. F., Gordon, R. G, et al. (1990). Current plate motions. *Geophys J*

[7] Hanks, T. C., & Kanamori, H. (1979). A moment-magnitude scale. *J. Geophys. Res.*, 84,

[8] Heki, K., Miyazaki, S., Takahashi, H., Kasahara, M., Kimata, F., Miura, S., & An, K. (1999). The Amurian plate motion and current plate kinematics in East Asia. *J. Geo‐*

[9] Herring, A. (2002). GLOBK global Kalman filter VLBI and GPS analysis program,

[10] Jin, S. G., & Zhu, W. (2002). Present-day spreading motion of the mid-Atlantic ridge.

[11] Jin, S. G., & Zhu, W. Y. (2003). Active Motion of Tectonic Blocks in Eastern Asia: Evi‐ dence from GPS Measurements. *ACTA Geological Sinica-English Edition*, 77(1), 59-63.

[12] Jin, S. G., & Park, P. (2006). Crustal stress and strain energy density rates in South Korea deduced from GPS observations. *Terr. Atmos. Ocean. Sci*, 17(1), 169-178.

[13] Jin, S. G., Park, P., & Zhu, W. (2007). Micro-plate tectonics and kinematics in North‐ east Asia inferred from a dense set of GPS observations. *Earth Planet. Sci. Lett*,

[14] Jin, S. G., Park, P., & Park, J. (2007 b). Why is the South Korean peninsula largely

version 10.0, Mass. *Inst. of Technol., Cambridge Mass, USA*.

*Chin. Sci. Bull*, 47(18), 1551-1555, 10.1360/02tb9342.

257(3-4), 486-496, 10.1016/j.epsl.2007.03.011, 2007a.

aseismic? Geodetic evidences. *Curr. Science*, 93(2), 250-253.

block modeling. *Geophys. Res*, 33, L11303, 10.1029/2006GL026077.


## **Earth Rotation – Basic Theory and Features**

### Sung-Ho Na

[30] Stein, S., & Gordon, R. (1984). Statistical tests of additional plate boundaries from

[31] Straub, C. (1996). Recent crustal defoirmation and strain accumulation in the Mar‐ mara sea region, inferred from GPS measurements. *Ist. of Geod. and Photogram. FTHZ*

[32] Thatcher, W. (2007). Microplate model for the present-day deformation of Tibet. *J.*

[33] Weber, J., Stein, S., & Engeln, J. (1998). Estimation of intraplate strain accumulation in the New Madrid seismic zone from repeat GPS surveys. *Tectonics*, 17, 250-266.

[34] Wei, D., & Seno, T. (1998). Determination of the Amurian plate motion, in Mantle dy‐ namics and plate interaction in East Asia. *edited by M. Flower, S. Chung, C. Lo and T.*

[35] Wessel, P., & Smith, W. H. F. (1998). New, improved version of Generic Mapping

[36] Zonenshain, L. P., & Savostin, L. A. (1981). Geodynamics of the Baikal rift zone and

plate motion inversions. *Earth Planet. Sci. lett*, 69, 401-412.

*Geophys. Res*, 112, B01401, 10.1029/2005JB004244.

Tools released: Eos. *Trans. Amer. Geophys. Union*, 79, 579.

plate tectonics of Asia. *Tectonophysics*, 76, 1-45.

*Mitt.*, 58.

284 Geodetic Sciences - Observations, Modeling and Applications

*Lee*, 337-346.

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/54584

### **1. Introduction**

Earth rotation is, in most case, meant to be spin rotation of the Earth. Diverse seasonal varia‐ tions are the result of Earth's obliquity to the ecliptic and orbital rotation of the Earth around the Sun. Although the linear motion of Earth's orbital rotation is faster than that of Earth's spin rotation, human beings cannot easily recognize it except aberration, which was first noticed by Bradley. In the Renaissance and also ancient times, scholars comprehended that orderly move‐ ments of stars in night sky are due to the Earth's spin rotation. To them, the planets, Moon and Sun moved with different periodicities in complicated way on the celestial sphere. As a matter of fact, the Earth's orbital/spin rotations are quite stable, and their stabilities far exceed human perception and most man-made instruments. However, the variations in both angular speed and direction of Earth spin have become detectable as technologies improved. Earth rotation is one of the most interesting scientific phenomena ever known. Moreover the importance of ac‐ curate knowledge of Earth rotation cannot be overestimated, because both the spatial and time systems of human civilization are referenced to the Earth and its spin rotational state.

Precession of equinox has been known since early day astronomy. The period of Earth's pre‐ cession is about 26 thousand years. Nutation, which is periodic perturbation of the Earth's spin axis, is much smaller in amplitude than precession and often regarded as associated motion of precession. While precession and nutation are the motions of the Earth's spin axis (angular mo‐ mentum) viewed by an observer in space outside of the Earth, the pole of Earth's spin rotation also changes with respect to observer on the Earth. Nowadays the Earth's rotational pole is usually represented by the Celestial Intermediate Pole (CIP). Slow drift of the Earth's pole in a time scale of thousand years is called polar wander. The polar wander is mostly due to the gla‐ cial isostatic adjustment and slow internal processes in the Earth. The pole of Earth's reference ellipsoid coincides with the Reference Pole, which was determined as average position of the Earth's pole between 1900 and 1905 (formerly called the Conventional International Origin). Pole position (CIP) with respect to the Reference Pole is represented by rectangular coordi‐ nates (*x*p, *y*p), which slowly draws a rough circle on the Earth's surface in a year or so. This pole offset is termed as polar motion.

© 2013 Na; licensee InTech. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2013 Na; licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Not only the orientation of the Earth's spin rotation but also its speed is slightly variable. Before the advent of quartz clocks, rotation of the Earth was regarded as the most reliable clock except planetary motion. Periodic variations in the length of day (LOD) have been found since the observations of star transit across meridian with accuracy better than 1 millisecond in the 1950s. Annual, semi-annual, and fortnightly perturbations are the first ones identified in the LOD variation. The Universal Time (UT), which is based on Earth rotation has been replaced by the Coordinated Universal Time (UTC) or the International Atomic Time (TAI), which are based on atomic clocks. Secular deceleration of the Earth due to tidal friction has been presumed since G. Darwin, and was confirmed later on.

Space geodetic technology since 1980s greatly enhanced precision and accuracy of measuring Earth rotation. Very Long Baseline Interferometry (VLBI) on the electromagnetic waves, which were emitted from quasar and detected at stations on Earth, has provided most important dataset, particularly for UT variation. VLBI and others as Global Positioning System (GPS), Satellite Laser Ranging (SLR) are capable of measuring pole offset with precision better than 1 milliarcsec as subdaily basis. Recently Ring Laser Gyroscope emerged as a unique and promising instrument to measure directly the Earth's spin rotational angular velocity with unprecedented accuracy.

Any position on the Earth's surface is usually denoted by its latitude and longitude in the Terrestrial Reference Frame (TRF), while direction to an astronomical object can be conven‐ iently represented by its right ascension and declination in the Celestial Reference Frame (CRF). Transformation between TRF and CRF is attained with the necessary information; Earth rotation angle due to time passage from a certain epoch and small changes in orientation due to precession-nutation and polar motion.

In this chapter, characteristics of aforementioned variations in the Earth's rotation are com‐ piled. Explanation of theoretical principle and observational features for each aspects are given one by one. But at times of much elaboration, more thorough treatment is avoided and proper references are recommended instead. The pre-requisite mechanics and mathematics are summarized in Appendix, where elementary vector algebra, harmonic oscillator, and basics for rotational mechanics *etc.* are covered. Old but still worth-reading monographs about Earth rotation are Munk and Macdonald [1], Lambeck [2], and Moritz and Mueller [3]. Explanations on terms and concepts of Earth rotation can be found in Seidelmann [4]. Extensive description about recent developments of Earth rotation study with emphasis on wobble and LOD was given by Gross [5]. In the same volume, Dehant and Mathews gave summary of Earth rotation theories with emphasis on nutation [6]. A technical note of International Earth Rotation and Reference Systems Service; IERS Conventions 2010 gives descriptions of new definitions in Earth rotation [7].

### **2. Precession and nutation**

Precession of equinox has been noticed by careful observers since early days. Earth's spin rotational pole and vernal/autumnal equinoxes slowly change their positions on the celestial sphere. Like any other gyroscopic motions, the equation of motion for precession of the Earth can be simply expressed as follows.

Not only the orientation of the Earth's spin rotation but also its speed is slightly variable. Before the advent of quartz clocks, rotation of the Earth was regarded as the most reliable clock except planetary motion. Periodic variations in the length of day (LOD) have been found since the observations of star transit across meridian with accuracy better than 1 millisecond in the 1950s. Annual, semi-annual, and fortnightly perturbations are the first ones identified in the LOD variation. The Universal Time (UT), which is based on Earth rotation has been replaced by the Coordinated Universal Time (UTC) or the International Atomic Time (TAI), which are based on atomic clocks. Secular deceleration of the Earth due to tidal friction has been presumed

Space geodetic technology since 1980s greatly enhanced precision and accuracy of measuring Earth rotation. Very Long Baseline Interferometry (VLBI) on the electromagnetic waves, which were emitted from quasar and detected at stations on Earth, has provided most important dataset, particularly for UT variation. VLBI and others as Global Positioning System (GPS), Satellite Laser Ranging (SLR) are capable of measuring pole offset with precision better than 1 milliarcsec as subdaily basis. Recently Ring Laser Gyroscope emerged as a unique and promising instrument to measure directly the Earth's spin rotational angular velocity with

Any position on the Earth's surface is usually denoted by its latitude and longitude in the Terrestrial Reference Frame (TRF), while direction to an astronomical object can be conven‐ iently represented by its right ascension and declination in the Celestial Reference Frame (CRF). Transformation between TRF and CRF is attained with the necessary information; Earth rotation angle due to time passage from a certain epoch and small changes in orientation due

In this chapter, characteristics of aforementioned variations in the Earth's rotation are com‐ piled. Explanation of theoretical principle and observational features for each aspects are given one by one. But at times of much elaboration, more thorough treatment is avoided and proper references are recommended instead. The pre-requisite mechanics and mathematics are summarized in Appendix, where elementary vector algebra, harmonic oscillator, and basics for rotational mechanics *etc.* are covered. Old but still worth-reading monographs about Earth rotation are Munk and Macdonald [1], Lambeck [2], and Moritz and Mueller [3]. Explanations on terms and concepts of Earth rotation can be found in Seidelmann [4]. Extensive description about recent developments of Earth rotation study with emphasis on wobble and LOD was given by Gross [5]. In the same volume, Dehant and Mathews gave summary of Earth rotation theories with emphasis on nutation [6]. A technical note of International Earth Rotation and Reference Systems Service; IERS Conventions 2010 gives descriptions of new definitions in

Precession of equinox has been noticed by careful observers since early days. Earth's spin rotational pole and vernal/autumnal equinoxes slowly change their positions on the celestial

since G. Darwin, and was confirmed later on.

286 Geodetic Sciences - Observations, Modeling and Applications

to precession-nutation and polar motion.

unprecedented accuracy.

Earth rotation [7].

**2. Precession and nutation**

$$\frac{d\vec{L}}{dt} = \vec{\tau} \tag{1}$$

Precessional torque comes mainly from the Moon and Sun; their gravitational pull to the Earth's equatorial bulge. But planets in the solar system also affect and contribute to the precession of the Earth in small amount. Earth's precession caused by all those effects together is called general precession, while luni-solar precession is referred to that caused by the Moon and Sun only. Recently it has been recognized that planetary perturbations lead to minute change in Earth's orbital plane, *i.e.* ecliptic. It is now recommended to use terms as precession of equator and precession of ecliptic [7]. Compared with the angular speed of the Earth's precession, lunar orbital rotation and apparent solar annual rotation around the Earth are much faster (similar comparison holds for the 18.6 years lunar orbital precession). Therefore, the luni-solar precessional torque can be evaluated by treating the masses of the Moon and Sun as hula-hoop shaped distribution (Figure 1).

**Figure 1.** Due to the Earth's equatorial bulge, the Moon and Sun exert torque to Earth's spin rotation. The resulting slow precession is clockwise viewed from above the pole of ecliptic.

In fact, the precessional torque is exerted by differential gravitational force (= tidal force) on the Earth's equatorial bulge. Qualitative understanding may be directly acquired from Figure 1. If one considers both directions of the Earth's spin angular momentum *L* <sup>→</sup>(solid arrow in the Figure) and the time averaged solar precessional torque *τ* <sup>→</sup> (the direction of the torque is perpendicular and into the paper), then, according to Equation; *ΔL* <sup>→</sup> <sup>=</sup> *<sup>τ</sup>* <sup>→</sup> *Δt*, the small change *ΔL* <sup>→</sup> of *<sup>L</sup>* <sup>→</sup> in time increment *Δt* should be *<sup>τ</sup>* <sup>→</sup> *Δt*. Average solar torque for the Earth's precession can be written as

$$\tau\_s = \frac{3GM\_s}{2r\_s^3} (\text{C} - A) \sin \varepsilon \cos \varepsilon \tag{2}$$

where *G*, *Ms*, *rs*, and *ε*are the constant of gravitation, mass of the Sun, mean solar distance, and the Earth's obliquity angle to the ecliptic plane. *C* and *A* are the principal moments of inertia of the Earth. Average lunar torque can be written similarly. Since the magnitude of Earth's spin angular momentum is *L* =*Cω*0, and the precession locus is shortened by a factor of sin*ε*, the total angular speed of the Earth's precession due to both of the lunar and solar torques is given as the following.

$$
\rho o\_{L.S.} = \frac{\Im G}{2o o\_E} \frac{\mathbb{C} - A}{\mathbb{C}} \left( \frac{M\_m}{r\_m^3} + \frac{M\_s}{r\_s^3} \right) \cos \varepsilon \tag{3}
$$

Derivation of Equations 2-3 can be found in [8-10], or other equivalent material. The angular velocity of general precession was recently found as *ωG*.*P*.=50.287946" / *yr*.

Due to its precession, orientation of the Earth's spin axis ceaselessly changes in the celestial sphere. This change can be represented by the three angles *ζ*, *θ*, and *z* as illustrated in Fig‐ ure 2. The transformation matrix corresponding to three successive rotations with these an‐ gles is represented as *Rprec* =*R*3(− *z*)*R*2(*θ*)*R*3(−*ζ*). Followings are the three angles of IAU2000A in terms of *T*, which is the time of given epoch in Julian century since J2000.0 - Greenwich noon on January 1st of year 2000 [11].

**Figure 2.** Representation of the Earth's precession by three angles ζ, θ, and *z*. The associated coordinate transforma‐ tion matrix is given as *R*3(− *z*)*R*2(θ)*R*3(−ζ).

```
2345
              2345
 2.5976176" + 2306.0809506" 0.3019015" 0.0179663" 0.0000327" 0.0000002" 
 2004.1917476" 0.4269353" 0.0418251" 0.0000601" 0.0000001" 
 2.5976176" + 230
            TT T T T
       TT T T T
z
z
q
 = ++--
 = ----
= - 2 345 6.0803226" 1.0947790" 0.0182273" 0.0000470" 0.0000003" T T TTT + + +-
                                           (4)
```
3

*s*

*r*

= -

*s*

t

288 Geodetic Sciences - Observations, Modeling and Applications

*L S*

w

Greenwich noon on January 1st of year 2000 [11].

tion matrix is given as *R*3(− *z*)*R*2(θ)*R*3(−ζ).

torques is given as the following.

*s*

*GM C A*

. . 3 3

w

velocity of general precession was recently found as *ωG*.*P*.=50.287946" / *yr*.

<sup>3</sup> ( )cos <sup>2</sup> *m s*

Derivation of Equations 2-3 can be found in [8-10], or other equivalent material. The angular

Due to its precession, orientation of the Earth's spin axis ceaselessly changes in the celestial sphere. This change can be represented by the three angles *ζ*, *θ*, and *z* as illustrated in Fig‐ ure 2. The transformation matrix corresponding to three successive rotations with these an‐ gles is represented as *Rprec* =*R*3(− *z*)*R*2(*θ*)*R*3(−*ζ*). Followings are the three angles of IAU2000A in terms of *T*, which is the time of given epoch in Julian century since J2000.0 -

**Figure 2.** Representation of the Earth's precession by three angles ζ, θ, and *z*. The associated coordinate transforma‐

 e


*E m s GC A M M C r r*

<sup>3</sup> ( )sin cos <sup>2</sup>

where *G*, *Ms*, *rs*, and *ε*are the constant of gravitation, mass of the Sun, mean solar distance, and the Earth's obliquity angle to the ecliptic plane. *C* and *A* are the principal moments of inertia of the Earth. Average lunar torque can be written similarly. Since the magnitude of Earth's spin angular momentum is *L* =*Cω*0, and the precession locus is shortened by a factor of sin*ε*, the total angular speed of the Earth's precession due to both of the lunar and solar

e e

(2)

Change in the right ascension and declination of any fixed position on the celestial sphere can readily be acquired with known precession matrix.

Besides precession, there exist small periodic oscillations superposed on it, and these are called nutations. The cause of nutation is the luni-solar gravitational pull to Earth's equatorial bulge, same as for precession but only with different frequency. In Figure 3, the main nutation (18.6 year period) is illustrated as a superposition on the precession. There are numerous different periodic components in the Earth's nutation. Treating the lunar and solar masses as circularly distributed rings was only an approximation to consider precession, but the true motions of the Moon and Sun relative to the Earth definitely are periodic. Accordingly effect of the periodic lunar or solar torque exist as oscillatory perturbations, *i.e*., nutations. Like the response of a simple harmonic oscillator to periodic driving force (away from resonance), which is proportional to inverse square of the driving frequency, amplitude of longer period nutation should be larger.

The nutational torque vector and resultant change in angular momentum vector can be briefly described as follows (here, precession itself is neglected for the time being. for more detailed treatment, see [3, 6]). Denote the nutational angular momentum and torque as *L* ⇀ *<sup>n</sup>* and *τ* ⇀ *n*. Write the nutational torque vector as elliptically rotating in time; *τ* ⇀ *<sup>n</sup>* = *e* ⌢ <sup>1</sup>*A*1cos*ωnt* ∓ *e* ⌢ <sup>2</sup>*A*2sin*ωnt* ('− ' for retrograde, '+ ' for prograde). Then steady state solution for *L* ⇀ *<sup>n</sup>* of the equation of motion; *d L* ⇀ *n dt* <sup>=</sup>*<sup>τ</sup>* ⇀ *<sup>n</sup>*, readily follows as *L* ⇀ *<sup>n</sup>* = *e* ⌢ 1 *A*1 *ωn* sin*ωnt* ± *e* ⌢ 2 *A*2 *ωn* cos*ωnt*

('+ ' for retrograde, '− ' for prograde). These are illustrated in Figure 4. Viewed from above, the motions of each set of retrograde/prograde nutational torque and angular momentum are clockwise/anticlockwise rotations, and nutational torque is ahead of phase *π* / 2 for both cases. The 18.6 year nutation is retrograde. Other major nutational components, such as fortnightly and semiannual, are mostly prograde, since the orbital rotations of the Moon and Sun (appa‐ rent motion relative to the Earth) are prograde. The number 1, 2, and 3 in Figure 4 represent successive times separated by a time increment *Δt*, and angular momentum change is given as *ΔL* ⇀ *<sup>n</sup>* ≃*τ* ⇀ *<sup>n</sup>Δt*. The nutational angle is given as arctan(*L <sup>n</sup>* / *Cωe*)≃ *L <sup>n</sup>* / *Cωe*. While torque itself is the same for different nutational components of same origin (lunar or solar), lower frequency nutation evidently has larger amplitude due to the inverse frequency dependence of nutation angular momentum; *L <sup>n</sup>* = | *L* ⇀ *<sup>n</sup>* | = | *τ<sup>n</sup>* | / *ωn*. This is the reason why 18.6 year nutation has much larger amplitude than other short period ones. Likewise semiannual nutation is roughly six times larger than fortnightly one (remember lunar tide is twice larger than solar tide). However, semi-annual nutation is of larger amplitude than annual one (particularly for the case of *Δε*). Similarly fortnightly nutation amplitude is larger than monthly one. These are due to symmetric nature of tidal force. If the eccentricities of lunar/solar orbits were smaller, the monthly and annual nutations would have been virtually nil.

**Figure 4.** Schematic illustrations of precessional/nutational angular momentum and nutational torque. On the left, total angular momentum is shown. In the middle, nutation angle is described by two angular momentum vectors. On the right, two cases of nutation. (above: retrograde elliptical, below: prograde circular) are shown. Nutational torque vectors are drawn in olive. The two unit vectors *e* ^ <sup>1</sup> and *e* ^ <sup>2</sup> are in the direction of longitude decrease and obliquity in‐ crease respectively. For convenience, half cycle phase is added to vectors of prograde case.

Two angles used to specify nutation are denoted as *Δε* and *Δψ*, which are illustrated in Figure 5.

**Figure 5.** Nutation described by two angles Δε and Δψ.

('+ ' for retrograde, '− ' for prograde). These are illustrated in Figure 4. Viewed from above, the motions of each set of retrograde/prograde nutational torque and angular momentum are clockwise/anticlockwise rotations, and nutational torque is ahead of phase *π* / 2 for both cases. The 18.6 year nutation is retrograde. Other major nutational components, such as fortnightly and semiannual, are mostly prograde, since the orbital rotations of the Moon and Sun (appa‐ rent motion relative to the Earth) are prograde. The number 1, 2, and 3 in Figure 4 represent successive times separated by a time increment *Δt*, and angular momentum change is given

is the same for different nutational components of same origin (lunar or solar), lower frequency nutation evidently has larger amplitude due to the inverse frequency dependence of nutation

much larger amplitude than other short period ones. Likewise semiannual nutation is roughly six times larger than fortnightly one (remember lunar tide is twice larger than solar tide). However, semi-annual nutation is of larger amplitude than annual one (particularly for the case of *Δε*). Similarly fortnightly nutation amplitude is larger than monthly one. These are due to symmetric nature of tidal force. If the eccentricities of lunar/solar orbits were smaller, the

**Figure 4.** Schematic illustrations of precessional/nutational angular momentum and nutational torque. On the left, total angular momentum is shown. In the middle, nutation angle is described by two angular momentum vectors. On the right, two cases of nutation. (above: retrograde elliptical, below: prograde circular) are shown. Nutational torque

<sup>2</sup> are in the direction of longitude decrease and obliquity in‐

^ <sup>1</sup> and *e* ^

crease respectively. For convenience, half cycle phase is added to vectors of prograde case.

*<sup>n</sup>Δt*. The nutational angle is given as arctan(*L <sup>n</sup>* / *Cωe*)≃ *L <sup>n</sup>* / *Cωe*. While torque itself

*<sup>n</sup>* | = | *τ<sup>n</sup>* | / *ωn*. This is the reason why 18.6 year nutation has

as *ΔL* ⇀ *<sup>n</sup>* ≃*τ* ⇀

angular momentum; *L <sup>n</sup>* = | *L*

290 Geodetic Sciences - Observations, Modeling and Applications

vectors are drawn in olive. The two unit vectors *e*

⇀

monthly and annual nutations would have been virtually nil.

In Table 1, several nutation components of largest amplitude are listed.


**Table 1.** Amplitudes of six largest nutation components [unit: milliarcsec]

The amplitude ellipses of four major components of Earth's nutation are shown in Figure 6.

The transformation matrix corresponding to the nutation described by the two angles *Δε* and *Δψ* is given as *Rnut* =*R*1(−*ε* −*Δε*)*R*3(−*Δψ*)*R*1(*ε*). The period of the largest nutational com‐ ponent is 18.6 years, which corresponds to the retrograde precession of the lunar orbital plane. Next largest components are semiannual, fortnightly, annual and monthly nutations. The most recent model for the precession and nutation has been reported by Mathews *et al*. [12], which is adopted by International Astronomical Union (IAU) as current standard mod‐ el [11-12]. Although out-dated, nutation calculated by using rigid Earth approximation of former studies [13-14] show close match with the most recent one [12]. This is due to the fact that the Earth behaves as an almost perfectly rigid body at such slow variations as the pre‐ cession and nutations. It is noted that the model of Mathews *et al.* was preceded by several other nutation models for nonrigid Earth. For example, see [15-18].

**Figure 6.** Comparison on size and shape of the nutation ellipses of four major components; 18.6 year, semiannual, fortnightly, and annual. A scale arrow of 5 arcsec is drawn for reference.

There are other ways to represent the coordinate transform associated with the Earth's precession and nutation. One of them is to use the position of pole projected onto the ICRF plane, of which orientation to quasars is invariable. IAU 2006 recommended to use a set of associated new terms and definitions as described below. For details, see [7] and references therein.

The transformation matrix, which incorporates both precession and nutation together, can be written as follows.

$$R\_{prec\,\ast \,mut} = \begin{pmatrix} 1 - b\,\, X \,\, ^2 & -b\,\, XY & -\,\, X \\ -b\,\, \mathbf{X}\,\, Y & \mathbf{1} - b\,\, Y \,\, ^2 & -\,\, Y \\ \mathbf{X} & \,\, Y & \mathbf{1} - b\,\, (\mathbf{X}\,\, ^2 + \,\, Y\,\, ^2) \end{pmatrix}.$$

where *b*is defined as *b* =1 / (1 + *Z*) with *Z* = 1− *X* <sup>2</sup> −*Y* <sup>2</sup> . Followings are IAU 2006 expressions for *X* and *Y* without explicitly showing all the oscillatory terms of nutation.

$$\begin{aligned} X &= -0.016617'' + 2004.191898'' \ T - 0.4297829'' \ T^2 - 0.19861834'' \ T^3 + 0.000007578'' \ T^4 \\ &+ 0.0000059285'' \ T^5 + \text{oscillator terms} \\ Y &= -0.006951'' - 0.025896'' \ T - 22.4072747'' \ T^2 + 0.00190059'' \ T^3 + 0.001112526'' \ T^4 \\ &+ 0.0000001358'' \ T^5 + \text{oscillator terms} \end{aligned}$$

**Figure 7.** Celestial Intermediate Pole referenced to the CRF equator.

#### **3. Secular deceleration**

cession and nutations. It is noted that the model of Mathews *et al.* was preceded by several

**Figure 6.** Comparison on size and shape of the nutation ellipses of four major components; 18.6 year, semiannual,

There are other ways to represent the coordinate transform associated with the Earth's precession and nutation. One of them is to use the position of pole projected onto the ICRF plane, of which orientation to quasars is invariable. IAU 2006 recommended to use a set of associated new terms and definitions as described below. For details, see [7] and references

The transformation matrix, which incorporates both precession and nutation together, can be

. Followings are IAU 2006 expressions

) )

*X* = − 0.016617" + 2004.191898" *T* − 0.4297829" *T* <sup>2</sup> − 0.19861834" *T* <sup>3</sup> + 0.000007578" *T* <sup>4</sup>

*Y* = − 0.006951" − 0.025896" *T* − 22.4072747" *T* <sup>2</sup> + 0.00190059" *T* <sup>3</sup> + 0.001112526" *T* <sup>4</sup>

for *X* and *Y* without explicitly showing all the oscillatory terms of nutation.

other nutation models for nonrigid Earth. For example, see [15-18].

292 Geodetic Sciences - Observations, Modeling and Applications

fortnightly, and annual. A scale arrow of 5 arcsec is drawn for reference.

1−*bX* <sup>2</sup> −*bXY* − *X* −*bXY* 1−*bY* <sup>2</sup> −*Y*

+ 0.0000059285" *T* <sup>5</sup> + oscillatory terms

+ 0.0000001358" *T* <sup>5</sup> + oscillatory terms

*X Y* 1−*b*(*X* <sup>2</sup> + *Y* <sup>2</sup>

where *b*is defined as *b* =1 / (1 + *Z*) with *Z* = 1− *X* <sup>2</sup> −*Y* <sup>2</sup>

therein.

*Rprec*+*nut* =(

written as follows.

Although variation in the magnitude of Earth's spin rotation had not been detected easily in their days, some investigators, such as I. Kant and G. Darwin, suspected the Earth's deceler‐ ation. With careful reasoning only, they correctly concluded that the Earth should be secularly decelerated due to tidal friction in the oceans and solid Earth and that the lunar orbit should be modified accordingly. A schematic illustration for this secular interaction is given in Figure 7. In the figure, the two identical Earth tidal bulges exist with minute phase delay due to tidal friction. Amplitude of body tide in the Earth exceeds 20 cm in most area over the world, and that of ocean tide is usually larger. Phase lag of body tide is found to be a few degrees, and differs for each tidal constituent. Phase lag of ocean tide is known to vary largely at places. The associated energy dissipation in ocean and solid Earth exceeds 3.0 Terra Watt. Tidal torque, which is due to the gravitational pull from the tide raising body (either the Moon or the Sun) to the misaligned tidal bulges, reduces spin angular momentum of the Earth. Equal and opposite torque should exist and increase the angular momentum of orbital rotation. Total angular momentum of the Earth-Moon system can be approximately expressed as *L* =(*Mere* <sup>2</sup> <sup>+</sup> *Mmrm* 2 )*ω<sup>m</sup>* <sup>=</sup> *MeMm Me* + *Mm r* 2 *ωm*, where *r* is Earth-Moon distance, *re* and *rm* are two each distances from the center of mass (approximately, <sup>1</sup> <sup>81</sup> *<sup>r</sup>* and <sup>80</sup> <sup>81</sup> *<sup>r</sup>*), and *Me* and *Mm* are masses

of the Earth and Moon respectively. Since lunar orbital parameters satisfy Kepler's third law, increase of lunar angular momentum is accompanied with increase of the Earth-Moon distance and decrease of lunar angular velocity; from *r* <sup>3</sup> *ω<sup>m</sup>* <sup>2</sup> <sup>=</sup>*G*(*Me* <sup>+</sup> *Mm*), relation between the two increments *dr* and *dωm* is acquired as 3 *dr* /*r* = −2 *dω<sup>m</sup>* / *ωm*. Tide raised in the Moon by the Earth should be larger in amplitude than Earth tide, but its effect is quite smaller in the secular deceleration due to the spin-orbit synchronous rotation of the Moon. It should be noted that the above arguments are based on the approximation neglecting the Earth's obliquity and lunar orbital inclination to the ecliptic.

**Figure 8.** Schematic illustration for tidal deceleration of the Earth and Earth-Moon distance increase due to their tidal interaction.

Solar tide in the Earth is about half of lunar tide in amplitude, and it significantly contributes to Earth's deceleration (about 20%). However, unlike the lunar orbital change, the change in the Earth-Sun distance or Earth's orbital rotation period due to the Earth-Sun tidal interaction is too small to be accurately observed.

G. Darwin carried elaborate formulation and calculation for changes in the dynamical state of the Earth Moon system, and extended further arguments [19]. Faster rotation of the Earth in the geological past was later identified from paleontological evidences. Historical eclipse records also revealed positive indications about the Earth's deceleration. Direct confirmations of the secular changes in the Earth's rotational state and lunar orbit have become available after operation of accurate atomic clocks and lunar laser ranging. Lunar laser ranging started since Apollo 11 landing in 1969, and enabled direct estimate of the present lunar recession as 3.82 cm/yr. Other related studies are satellite orbit analysis and ocean tide modeling. The fact that the Moon was closer to the Earth in the past, led G. Darwin to the fission hypothesis for lunar origin. Other hypotheses suggested for lunar origin were capture, binary accretion, and impact theories. This interesting topic has been occasionally re-visited [20].

In these days, estimates about the Earth's deceleration rate by different approaches seem to converge, so that the LOD is now generally believed to increase with a rate of 1.8 millisecond / cy. Without glacial isostatic adjustment, this rate would be 2.3 millisecond /cy, however, there is large uncertainty in the variation of LOD mainly due to slow and complicated flows in the deep interior of the Earth. Recent fast retreat of glaciers and associated changes in the Earth's moments of inertia contribute to the variation of LOD by a certain amount, which is not accurately known.

deceleration due to the spin-orbit synchronous rotation of the Moon. It should be noted that the above arguments are based on the approximation neglecting the Earth's obliquity and lunar

**Figure 8.** Schematic illustration for tidal deceleration of the Earth and Earth-Moon distance increase due to their tidal

Solar tide in the Earth is about half of lunar tide in amplitude, and it significantly contributes to Earth's deceleration (about 20%). However, unlike the lunar orbital change, the change in the Earth-Sun distance or Earth's orbital rotation period due to the Earth-Sun tidal interaction

G. Darwin carried elaborate formulation and calculation for changes in the dynamical state of the Earth Moon system, and extended further arguments [19]. Faster rotation of the Earth in the geological past was later identified from paleontological evidences. Historical eclipse records also revealed positive indications about the Earth's deceleration. Direct confirmations of the secular changes in the Earth's rotational state and lunar orbit have become available after operation of accurate atomic clocks and lunar laser ranging. Lunar laser ranging started since Apollo 11 landing in 1969, and enabled direct estimate of the present lunar recession as 3.82 cm/yr. Other related studies are satellite orbit analysis and ocean tide modeling. The fact that the Moon was closer to the Earth in the past, led G. Darwin to the fission hypothesis for lunar origin. Other hypotheses suggested for lunar origin were capture, binary accretion, and

In these days, estimates about the Earth's deceleration rate by different approaches seem to converge, so that the LOD is now generally believed to increase with a rate of 1.8 millisecond / cy. Without glacial isostatic adjustment, this rate would be 2.3 millisecond /cy, however, there is large uncertainty in the variation of LOD mainly due to slow and complicated flows in the deep interior of the Earth. Recent fast retreat of glaciers and associated changes in the Earth's moments of inertia contribute to the variation of LOD by a certain amount, which is not

impact theories. This interesting topic has been occasionally re-visited [20].

orbital inclination to the ecliptic.

294 Geodetic Sciences - Observations, Modeling and Applications

is too small to be accurately observed.

interaction.

accurately known.

**Figure 9.** Recent 40 years data of delay in the universal time (UT) and Earth's angular velocity; cumulative delay in UT1 (top), difference between UT1 and UTC with leap seconds seen as steps (middle), and the excessive Earth's angular velocity deduced from UT1 delay (bottom).

In Figure 9, recent variation of UT1 is illustrated by three different ways. Above, cumulative delay in UT1 is shown, and then, in the middle, difference between UT1 and UTC (UT1-UTC) is shown. Including the last leap second introduced at the midnight of June 30th of 2012, there were total 25 leap seconds since 1972. The bottom graph shows the excessive amount of Earth's spin angular speed *Δω* from its nominal value *ω*<sup>0</sup> =7.292115×10−<sup>5</sup> (*rad* / *s*). Recent Earth's spin angular velocity, as can be seen in Figure 9(bottom), has been increasing in minute amount. This tendency can also be noticed from reduced number of leap seconds during the last twelve years or so.

Energy dissipation and angular momentum transfer via the tidal interaction process can be expressed simply in terms of *dr*(infinitesimal increase in the Earth-Moon distance) for isolated Earth-Moon system with neglecting Earth's obliquity and lunar inclination. The work involved in the process of angular momentum transfer is given as *dW* = −*d*( 1 <sup>2</sup> *Ieω<sup>e</sup>* 2 )= − *Ieωedω<sup>e</sup>* = −*τeωedt*. After a little algebra, this can re-written as *dW* =( *MeMmr* 2(*<sup>M</sup>* <sup>+</sup> *<sup>M</sup>* ) <sup>−</sup> 3*Im* <sup>2</sup>*<sup>r</sup>* )*ωmωedr*. The mechanical energy increase *dE* is given as *dE* <sup>=</sup> <sup>−</sup> <sup>1</sup> 2 ( *MeMm Me* + *Mm* )*r* − 3*Im* 2*r ω<sup>m</sup>* 2 *dr* + *GMeMm <sup>r</sup>* <sup>2</sup> *dr*. The difference between *dW* and *dE*, *i.e.*, *dQ* =*dW* −*dE*is the dissipated energy during the process. In fact, the majority of the work (96.3%) is dissipated into heat. This argument is valid regardless of the places where the dissipation occur. The estimates for the work rate *dW dt* and the tidal torque *τ* of present days are found as 3.30×10<sup>12</sup> *Watt* and 4.52×10<sup>16</sup> *N* ⋅*m*. More elaborate calculation with including Earth's obliquity, lunar inclination, and solar tidal dissipation leads to slightly different estimates; 96.6%, 3.56×10<sup>12</sup> *Watt*, and 4.88×10<sup>16</sup> *N* ⋅*m*. The corresponding values for present Earth deceleration and lunar orbital retardation are *dω<sup>e</sup> dt* <sup>=</sup> <sup>−</sup>6.08×10−<sup>22</sup> *rad* / *<sup>s</sup>* <sup>2</sup> and

$$\frac{d\,\omega\_m}{dt} = -1.26 \times 10^{-23} \,\mathrm{rad} \,\left| \,\mathrm{s}^2 = -25.9'' \right| \,\mathrm{c} \,\mathrm{y}^2.$$

Analysis on satellite orbit can yield estimate of the tidal torque, which decelerates the Earth. If there were not any other force exerting on a satellite rather than Earth's central gravita‐ tional attraction, its orbit would be a perfect ellipse. Due to luni-solar gravitational attrac‐ tion, solar radiation pressure, and others, satellites undergo certain changes in their orbital configurations. Although the Earth's gravity field perturbation due to its tidal deformation is not quite large, satellite tracking has been precise enough to detect such effect on satellite orbits since late 70s [21-23]. By extending this kind approach to the Moon, the only natural satellite of the Earth, they estimated lunar orbital retardation as 27.4, 25.3, and 24.9 " / *cy* <sup>2</sup> [21-23]. Christodoulidis *et al*. accordingly estimated the present Earth deceleration rate as −5.98×10−<sup>22</sup> *rad* /*s* <sup>2</sup> [22].

Since the late 60s, ocean tide modeling has been attempted, and global ocean tide models were acquired as numerical solution of Laplace's tidal equation with grid spacing larger than 2 ×2 . By using acquired global ocean tide model, it is possible to calculate the tidal torque that decelerate the Earth. Estimates for Earth's secular deceleration based on the first-genera‐ tion ocean tide models were summarized by Lambeck [24]. With fast development of com‐ puting devices, much more extensive modelings have been acquired lately, such as NAOJ99, FES2004, EOT11, *etc*. In fact, better performance of recent models is feasible due also to sea surface height data derived from satellite altimetry (TOPEX/POSEIDON). Ray *et al*. combined the satellite altimetry and orbit analysis and estimated the Earth's deceleration rate as −1304" / *cy* <sup>2</sup> , which is −6.35×10−<sup>22</sup> *rad* /*s* <sup>2</sup> and corresponds to LOD increase rate of 2.37 ms/cy [25]. Table 2 is based on their results.


**Table 2.** Estimate of Earth's secular deceleration by each ocean tide constituents. [unit: arcsec/cy2] [25].

Accumulation of universal time delay *ΔT* due to a constant LOD increase of 2.0 millisecond/cy for a century corresponds to 36.525s; *Δ<sup>T</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup> ×36525×0.002*<sup>s</sup>* =36.525*s*. A thousand year accu‐ mulation (ten times longer) by the same rate would result in *ΔT* = 3652.5 s (hundred times) de‐ lay, which is about one hour. Similarly two thousand year accumulation would result in about four hours delay. In fact, a correction term linear to time span should be added unless the time span is measured from the reference epoch, which is the year AD 1820 [26]. Astronomical evi‐ dences about the Earth's deceleration and lunar orbital retardation can be categorized into two different kinds of records; (i)'telescopic observation,' which has been carried for a few hun‐ dred years, and (ii) 'solar eclipse record,' which can be found from thousands year old Chinese literature or Babylonian inscriptions on clay tablets. Figure 10 is a redrawing of *ΔT* curve based on the observations of transit of Mercury and lunar occultation [27].

*dE* <sup>=</sup> <sup>−</sup> <sup>1</sup> 2 (

*dω<sup>m</sup>*

−5.98×10−<sup>22</sup> *rad* /*s* <sup>2</sup>

−1304" / *cy* <sup>2</sup>

*MeMm Me* + *Mm* )*r* − 3*Im* 2*r ω<sup>m</sup>* 2 *dr* +

296 Geodetic Sciences - Observations, Modeling and Applications

*dt* <sup>=</sup> <sup>−</sup>1.26×10−<sup>23</sup> *rad* / *<sup>s</sup>* <sup>2</sup> <sup>=</sup> <sup>−</sup>25.9" / *cy* <sup>2</sup>

[22].

, which is −6.35×10−<sup>22</sup> *rad* /*s* <sup>2</sup>

ms/cy [25]. Table 2 is based on their results.

for a century corresponds to 36.525s; *Δ<sup>T</sup>* <sup>=</sup> <sup>1</sup>

dissipation occur. The estimates for the work rate *dW*

Earth deceleration and lunar orbital retardation are

*GMeMm*

.

*dQ* =*dW* −*dE*is the dissipated energy during the process. In fact, the majority of the work (96.3%) is dissipated into heat. This argument is valid regardless of the places where the

are found as 3.30×10<sup>12</sup> *Watt* and 4.52×10<sup>16</sup> *N* ⋅*m*. More elaborate calculation with including Earth's obliquity, lunar inclination, and solar tidal dissipation leads to slightly different estimates; 96.6%, 3.56×10<sup>12</sup> *Watt*, and 4.88×10<sup>16</sup> *N* ⋅*m*. The corresponding values for present

Analysis on satellite orbit can yield estimate of the tidal torque, which decelerates the Earth. If there were not any other force exerting on a satellite rather than Earth's central gravita‐ tional attraction, its orbit would be a perfect ellipse. Due to luni-solar gravitational attrac‐ tion, solar radiation pressure, and others, satellites undergo certain changes in their orbital configurations. Although the Earth's gravity field perturbation due to its tidal deformation is not quite large, satellite tracking has been precise enough to detect such effect on satellite orbits since late 70s [21-23]. By extending this kind approach to the Moon, the only natural satellite of the Earth, they estimated lunar orbital retardation as 27.4, 25.3, and 24.9 " / *cy* <sup>2</sup> [21-23]. Christodoulidis *et al*. accordingly estimated the present Earth deceleration rate as

Since the late 60s, ocean tide modeling has been attempted, and global ocean tide models were acquired as numerical solution of Laplace's tidal equation with grid spacing larger than 2 ×2 . By using acquired global ocean tide model, it is possible to calculate the tidal torque that decelerate the Earth. Estimates for Earth's secular deceleration based on the first-genera‐ tion ocean tide models were summarized by Lambeck [24]. With fast development of com‐ puting devices, much more extensive modelings have been acquired lately, such as NAOJ99, FES2004, EOT11, *etc*. In fact, better performance of recent models is feasible due also to sea surface height data derived from satellite altimetry (TOPEX/POSEIDON). Ray *et al*. combined the satellite altimetry and orbit analysis and estimated the Earth's deceleration rate as

**Tidal constituent** *M*<sup>2</sup> *K*<sup>1</sup> *S*<sup>2</sup> *O*<sup>1</sup> *N*<sup>2</sup> *P*<sup>1</sup> *K*<sup>2</sup> *Q*<sup>1</sup> Total **Acceleration** -919 -120 -73 -71 -40 -13 -10 -3 -1304

Accumulation of universal time delay *ΔT* due to a constant LOD increase of 2.0 millisecond/cy

**Table 2.** Estimate of Earth's secular deceleration by each ocean tide constituents. [unit: arcsec/cy2] [25].

*<sup>r</sup>* <sup>2</sup> *dr*. The difference between *dW* and *dE*, *i.e.*,

*dω<sup>e</sup>*

and corresponds to LOD increase rate of 2.37

<sup>2</sup> ×36525×0.002*<sup>s</sup>* =36.525*s*. A thousand year accu‐

*dt* and the tidal torque *τ* of present days

*dt* <sup>=</sup> <sup>−</sup>6.08×10−<sup>22</sup> *rad* / *<sup>s</sup>* <sup>2</sup> and

**Figure 10.** Values of Δ*T* between year 1627 and 1975 based on the observations of lunar occultation and transit of Mercury [27].

The parabolic trend of *ΔT* curve is due to the Earth's secular deceleration. As shown in Figure 10, the secular deceleration is superposed with large fluctuations, which are ascribed to processes in the Earth's core and mantle. From these data, Morrison determined lunar orbital retardation as *dω<sup>m</sup> dt* <sup>=</sup> <sup>−</sup>26" / *cy* <sup>2</sup> [27], which is quite close to recent and reliable estimate by lunar laser ranging and others. From vast amount of historical records of solar eclipse and others, an empirical formula was found for *ΔT* as follows;

$$
\Delta T = -20 + 31t^2 \tag{5}
$$

where *ΔT* is in second, and *t* is time from AD 1820 in century [26]. Later the quadratic coefficient was adjusted from 31 to 32 [28]. Two corresponding rates of LOD increase are 1.70 and 1.75 millisecond/cy. It is noted that there were other estimates of same sorts, which largely differ from those introduced above. For example, see [29-30]. In Figure 11, curve of *ΔT* determined from observations between 500 BC and AD 1300 [26], and an approximate curve of *ΔT* based on Morrison's table [28] are shown.

**Figure 11.** Two sets of Δ*T* data and curves: Δ*T* determined from observations between 500 BC and AD 1300 and their parabola fittings. Inset is set of data and fittings for later times [26] (left), and an approximate curve of Δ*T* based on tabulated data in [28] (right, Δ*T* + 20*s* is drawn in logarithmic scale).

Paleontological evidences exist for faster rotation of the Earth in the Mesozoic and Paleozoic era. As trees retain yearly growth rings, certain organisms, such as coral and shell, can record high and low ocean tides in their hard parts due to differential growth. Careful counting of these growth lines in fossil specimen yielded the numbers of days per month and year at geological past [31-33]. Lambeck estimated Earth's deceleration rate and lunar orbital retar‐ dation rate from those numbers [33]. His estimates based on fossil bivalve data were *dω<sup>e</sup> dt* <sup>=</sup> <sup>−</sup>5.9×10−<sup>22</sup> *rad* / *<sup>s</sup>* <sup>2</sup> and *dω<sup>m</sup> dt* <sup>=</sup> <sup>−</sup>1.3×10−<sup>23</sup> *rad* / *<sup>s</sup>* <sup>2</sup>≃27" / *cy* <sup>2</sup> . These estimates were proved to be quite reasonable by later studies. Therefore, tidal dissipation and associated lunar orbital change must have occurred rather consistently during hundreds of millions years regardless of decadal fluctuations due to core/mantle processes in the Earth or other pertur‐ bations. Followings are the formulation developed by Lambeck for the analysis of paleonto‐ logical data.

12 3 1 2 3 2 1 0 0 , , a 1 1 1 , , b ( ) 1 ( )/ 10 0 <sup>1</sup> ( ) <sup>1</sup> / <sup>010</sup> / *es es ms s m s s e e es m m s m s m s s e e s m m s e m NN N dN d dN d d dN d dt dt dt dt dt dt dt N t t t t t d dt t d dt* ww ww w w w ww w w w ww w w w w w w w w w w w w w w b b w w w -- - == = - - <sup>=</sup> - - é ù <sup>D</sup> + - é ùê ú ê ú <sup>D</sup> - ê ú -- = ê úê ú ê úê ú ë ûê ú ë û ; ; 2 0 3 0 1 ( )/ c 1 ( )/ *s e s s m s N t N t* w w w w w é ù ê ú + - + - ë û (6)

*N*1, *N*2, *N*3 are the numbers of solar days per year and synodic month, and the number of synodic month per year respectively. Time derivatives of these three numbers are expressed as in equation (6b)

differ from those introduced above. For example, see [29-30]. In Figure 11, curve of *ΔT* determined from observations between 500 BC and AD 1300 [26], and an approximate curve

**Figure 11.** Two sets of Δ*T* data and curves: Δ*T* determined from observations between 500 BC and AD 1300 and their parabola fittings. Inset is set of data and fittings for later times [26] (left), and an approximate curve of Δ*T* based on

Paleontological evidences exist for faster rotation of the Earth in the Mesozoic and Paleozoic era. As trees retain yearly growth rings, certain organisms, such as coral and shell, can record high and low ocean tides in their hard parts due to differential growth. Careful counting of these growth lines in fossil specimen yielded the numbers of days per month and year at geological past [31-33]. Lambeck estimated Earth's deceleration rate and lunar orbital retar‐ dation rate from those numbers [33]. His estimates based on fossil bivalve data were

*dt* <sup>=</sup> <sup>−</sup>1.3×10−<sup>23</sup> *rad* / *<sup>s</sup>* <sup>2</sup>≃27" / *cy* <sup>2</sup>

2

 w w

w

1 ( )/ c

1 1 1 , , b ( )

; ;

*e e es m m*

 w

> ww

w w

*e e s*

 w

*s*

, , a

1 0

3 0

*N t*

1 ( )/

w

+ - ë û

w

é ù ê ú + -

 w

*e s*

 w

2 0

w

*N t*

*m s*

 w

*N t*

0

w

*s*

*t*

*s m s m s s*

 w w

proved to be quite reasonable by later studies. Therefore, tidal dissipation and associated lunar orbital change must have occurred rather consistently during hundreds of millions years regardless of decadal fluctuations due to core/mantle processes in the Earth or other pertur‐ bations. Followings are the formulation developed by Lambeck for the analysis of paleonto‐

. These estimates were

(6)

of *ΔT* based on Morrison's table [28] are shown.

298 Geodetic Sciences - Observations, Modeling and Applications

tabulated data in [28] (right, Δ*T* + 20*s* is drawn in logarithmic scale).

and

12 3

 ww


 ww

*NN N*

<sup>1</sup> ( ) <sup>1</sup>

/ <sup>010</sup>

*t t*

*t*

 b

*s e*

*t*

w *dω<sup>m</sup>*

1 2 3

*dN d dN d d dN d dt dt dt dt dt dt dt*

w

w

*m m*

1 ( )/ 10 0

/

*m*

*d dt*

é ù <sup>D</sup> + - é ùê ú ê ú <sup>D</sup> - ê ú -- = ê úê ú ê úê ú ë ûê ú ë û

w

w

w

*d dt*

*es es ms s m s s*

> w w


w

*dω<sup>e</sup>*

logical data.

w

*dt* <sup>=</sup> <sup>−</sup>5.9×10−<sup>22</sup> *rad* / *<sup>s</sup>* <sup>2</sup>

ww

w

w

b

Write *ωe*(*t*)=*ωe*(*t*0) + *Δω<sup>e</sup>* + (*dω<sup>e</sup>* / *dt*) *t* and *ωm*(*t*)=*ωm*(*t*0) + *Δω<sup>m</sup>* + (*dω<sup>m</sup>* / *dt*) *t*, then the condition equation for the unknowns *dω<sup>e</sup> dt* , *dω<sup>m</sup> dt* , *Δωe*, and *Δωm* is written as in equation (6c) with *<sup>β</sup>* <sup>=</sup> *<sup>ω</sup>e*(*t*0)−*ω<sup>s</sup> ωm*(*t*0)−*ω<sup>s</sup>* . While the estimates of different studies based on the historical eclipse records *dω<sup>e</sup> dω<sup>m</sup>*

varied widely on two quantities *dt* and *dt* , the reliability of estimate based on paleonto‐ logical record is noticeable.

Numerical extrapolation for the past/future lunar orbit was first attempted by Darwin [19], and later by many other investigators including Goldreich, Migard, Hansen, and Webb [34-37]. Goldreich calculated the past state of the Earth-Moon system by using three time steps; year, period of lunar orbital precession, and secular [34]. Goldreich and Mignard avoided time scale problem by using the Earth-Moon distance as the independent variable. Calculations of Hansen and Webb were based on each idealized ocean tide modeling. Hansen specified two kinds of paleo-ocean configurations for his calculation; one - circular continent at polar region and the other - circular continent at equator [37]. Webb used the orientation averaged ocean with one continental cap. By extending Kaula's satellite orbit theory, Lambeck derived differential equations for lunar orbital evolution. By using those Lambeck's formulae with modification, the extrapolation was recently repeated [38], and its set of the calculated past Earth-Moon distance for different tidal friction parameter is shown in Figure 12. According to this result, the Moon could have been at close approach to the Earth at different times, from 1.6 - 4.0 billion years ago or even earlier depending on the amount of reduction in tidal friction in the past.

**Figure 12.** Calculated Earth-Moon distance in the past. According to the modified Lambeck's formulae with assuming power law behavior for tidal phase lag angle. *M*2 and *S*2 tidal constituents only are considered [38].

Calculated values of five parameters (eccentricity, obliquity, lunar orbital inclination, *ωe*, and *ωm*) in the future are shown in Figure 13. The scale factor is ratio of Earth-Moon distance to its present value, *i.e*., *r* /*r*0. As the ratio increases from 1.0 (present) to 1.2-1.3, Earth's spin angular velocity decreases almost linearly, and obliquity angle increases abruptly up to and over 60 degrees, while other parameters of lunar orbit gradually change in small amounts.

**Figure 13.** Calculated Earth-Moon system configuration in future by modified Lambeck's formulae. The horizontal axis is the ratio *r* / *r*0. Solid lines correspond to those including solar tidal dissipation, while dotted lines correspond to those not including it. The units are [108m], [none], [deg], [deg], [10-5 rad/s], and [10-6 rad/s] from left top.

There exist four hypotheses suggested for the lunar origin; (i) 'fission', (ii) 'capture', (iii) 'binary accretion', and (iv) 'impact.' It is difficult to verify/trace such a far distant past event like the Earth-Moon system formation, and one may consider only feasibilities of those hypotheses. A brief sketch of these four hypotheses is given here. Fission hypothesis, suggested by Darwin [19], is an explanation of lunar origin as a result of mechanical resonance of the early Earth, which should have been rotating very fast. Since the natural period of foot ball mode of Earth's free oscillation is 54 minutes, the mechanical resonance to eject a part of mantle might have been possible, provided the early Earth had been rotating with a period about two hours or so. If the Earth would gain the angular momentum, which was lost due to both lunar and solar tidal interactions during the whole past of billion years, this fast rotation could be possible. Fission theory is also compatible with the fact that the lunar mass density is quite close to that of the Earth's mantle. Problem with the fission theory is that direct calculation of the Earth-Moon system back to the past should lead to high inclination angle of the lunar orbit. One has to imagine a certain scenario, such as close approach of Venus to perturb lunar orbit, to reconcile the imperfection. Capture hypothesis states that the Moon was initially formed independently and later captured by the Earth's gravity field during its passage near the Earth. Suppose it approached by following a hyperbolic/parabolic orbit with respect to the Earth, there should have been a mechanism to explain the change of lunar orbit. Gerstenkorn event is one such explanation claiming the Moon had approached in a retrograde orbit and then undergone its orbital change from hyperbolic to elliptic due to severe lunar tidal dissipation. Problem with the Gerstenkorn event is that, even at closest distance, tidal interaction cannot perform with such an extreme efficiency. Could there be other alternative stopping mecha‐ nism, then the capture hypothesis may remain feasible. Binary accretion hypothesis states that the Moon was formed with the Earth together as two isolated bodies by planetesimal accu‐ mulation from the beginning stage of the Earth-Moon system. Had the Moon been formed as an isolated body from the Earth, the low density of the Moon cannot be explained. Impact hypothesis emerged in the late 70s and was numerically simulated later. Suppose an object of 0.1 Earth's mass or so had collided with the early Earth, the Moon could have been formed afterwards by accretion of the remaining fragments. Impact hypothesis has been favored, because the Moon is depleted of volatile elements.

Like other major natural satellites, such as Galilean satellites, the Moon is already in synchro‐ nous rotation so that it always shows the same face toward the Earth. This synchronous rotation of the Moon can be maintained due also to another tidal interaction - dissipation of tide in the Moon raised by the Earth. Unless other impact with large third body or comparable perturbations, the Earth-Moon system will undergo secular changes as described above; the Earth's spin will slow down with its obliquity increase, and the lunar orbit will become larger with small decrease in its inclination angle.

### **4. Liouville equation and excitation function**

Calculated values of five parameters (eccentricity, obliquity, lunar orbital inclination, *ωe*, and *ωm*) in the future are shown in Figure 13. The scale factor is ratio of Earth-Moon distance to its present value, *i.e*., *r* /*r*0. As the ratio increases from 1.0 (present) to 1.2-1.3, Earth's spin angular velocity decreases almost linearly, and obliquity angle increases abruptly up to and over 60

**Figure 13.** Calculated Earth-Moon system configuration in future by modified Lambeck's formulae. The horizontal axis is the ratio *r* / *r*0. Solid lines correspond to those including solar tidal dissipation, while dotted lines correspond to

There exist four hypotheses suggested for the lunar origin; (i) 'fission', (ii) 'capture', (iii) 'binary accretion', and (iv) 'impact.' It is difficult to verify/trace such a far distant past event like the Earth-Moon system formation, and one may consider only feasibilities of those hypotheses. A brief sketch of these four hypotheses is given here. Fission hypothesis, suggested by Darwin [19], is an explanation of lunar origin as a result of mechanical resonance of the early Earth, which should have been rotating very fast. Since the natural period of foot ball mode of Earth's free oscillation is 54 minutes, the mechanical resonance to eject a part of mantle might have

those not including it. The units are [108m], [none], [deg], [deg], [10-5 rad/s], and [10-6 rad/s] from left top.

degrees, while other parameters of lunar orbit gradually change in small amounts.

300 Geodetic Sciences - Observations, Modeling and Applications

In this section, the equation of motion for Earth's angular velocity perturbation is derived in the Earth fixed reference frame by neglecting elastic deformation and other complicated features of the real Earth. First, the general equation of motion for rigid body rotation and Euler's equation for wobble are derived. Further considerations will be followed after then.

Spin angular momentum of a rotating body can be written as *L* ⇀ <sup>=</sup>*∫<sup>r</sup>* ⇀ ×*v* ⇀ *dm*, where the center of its mass can be taken as the coordinate origin. Equation of motion states that the time rate of angular momentum equals to the exerting torque; *<sup>d</sup> <sup>L</sup>* ⇀ *dt* <sup>=</sup>*<sup>τ</sup>* ⇀ . In rotating coordinate frame, time derivative should be considered with retaining terms due to change in the direction of reference frame in space (in Appendix), so that the equation of rotational motion in body frame is rewritten as follows.

$$\frac{d\vec{L}}{dt} + \vec{o} \times \vec{L} = \vec{\tau} \tag{7}$$

Taking the principal axes of body as reference frame axes, then *L* ⇀ and *<sup>ω</sup>* ⇀ × *L* ⇀ are simply expressed as follows.

$$\begin{aligned} \mathbf{L}\_i &= \sum\_{j=1}^3 I\_{ij} \alpha\_j \implies I\_i \alpha\_i \text{ , i.e.,} \bar{\mathbf{L}} = \alpha\_1 I\_1 \hat{\mathbf{e}}\_1 + \alpha\_2 I\_2 \hat{\mathbf{e}}\_2 + \alpha\_3 I\_3 \hat{\mathbf{e}}\_3 \qquad \mathbf{a} \\\ \bar{\mathbf{a}} \times \bar{\mathbf{L}} &= (\alpha\_1 \hat{\mathbf{e}}\_1 + \alpha\_2 \hat{\mathbf{e}}\_2 + \alpha\_3 \hat{\mathbf{e}}\_3) \times (\alpha\_1 I\_1 \hat{\mathbf{e}}\_1 + \alpha\_2 I\_2 \hat{\mathbf{e}}\_2 + \alpha\_3 I\_3 \hat{\mathbf{e}}\_3) \\\ &= \alpha\_2 \alpha\_3 (I\_3 - I\_2) \hat{\mathbf{e}}\_1 + \alpha\_3 \alpha\_1 (I\_1 - I\_3) \hat{\mathbf{e}}\_2 + \alpha\_1 \alpha\_2 (I\_2 - I\_1) \hat{\mathbf{e}}\_3 \qquad \mathbf{b} \end{aligned} \tag{8}$$

where the principal moments of inertia *Ii* replace inertia tensor *Iij* (underlying basic summar‐ ized in Appendix).

Then Equation 7 is rewritten as

$$I\_1 \frac{d o\_1}{dt} + o\_2 o\_3 (I\_3 - I\_2) = \tau\_1 \qquad \mathbf{a}$$

$$I\_2 \frac{d o\_2}{dt} + o\_3 o\_1 (I\_1 - I\_3) = \tau\_2 \quad \mathbf{b} \tag{9}$$

$$I\_3 \frac{d o\_3}{dt} + o\_1 o\_2 (I\_2 - I\_1) = \tau\_3 \quad \mathbf{c}$$

In case two values *I*1 and *I*2 are same (*I*<sup>1</sup> = *I*<sup>2</sup> ⇒ *I*), then Equation 9c becomes simple; *I*3 *dω*<sup>3</sup> *dt* <sup>=</sup>*τ*3. Its solution can be expressed as *ω*3(*t*)=*ω*<sup>0</sup> <sup>+</sup> 1 *I*3 *∫ t*0 *t τ*(*t* ') *dt* '. The solution *ω*3(*t*) for sinusoidal torque of frequency *ωd* as *τ*(*t*)=*τ*0cos*ω<sup>d</sup> t* is readily acquired as *ω*3(*t*)=constant + *τ*0 *I*3*ω<sup>d</sup>* sin*ω<sup>d</sup> t*. As a further restriction, if the external torque does not exist (*τ* ⇀ =0), then *ω*3(*t*) should remain constant as *ω*<sup>3</sup> =*ω*<sup>0</sup> and the two Equations 9a-b become as follows.

$$I\frac{d\alpha\_1}{dt} + \alpha\_2 \alpha\_0 (I\_3 - I) = 0,\\ I\frac{d\alpha\_2}{dt} + \alpha\_0 \alpha\_1 (I - I\_3) = 0\tag{10}$$

The coupled solution for *ω*1 and *ω*2 of Equation 10 is found as a circular motion as follows.

time derivative should be considered with retaining terms due to change in the direction of reference frame in space (in Appendix), so that the equation of rotational motion in body frame

> *dL <sup>L</sup> dt*

Taking the principal axes of body as reference frame axes, then *L*

,.

 w

1

*dt*

w

*dt*

w

w

*dt*

*dt* <sup>=</sup>*τ*3. Its solution can be expressed as *ω*3(*t*)=*ω*<sup>0</sup> <sup>+</sup>

2

3

 w

å <sup>v</sup>

.,

= Þ =+ +

´= + + ´ + +

*L I I L Ie Ie Ie*

w w

+´= w

11 2 2 33 111 22 2 333 23 3 2 1 31 1 3 2 12 2 1 3

*I Ie I Ie I Ie*

( ˆˆ ˆ ˆ ˆ ) ( ) () () () ˆˆˆ

 w

 ww

*L e e e Ie Ie Ie*

= -+ -+ -

1 23 3 2 1

+ -=

+ -=

+ -=

*<sup>d</sup> I II*

ww

*<sup>d</sup> I II*

ww

*<sup>d</sup> I II*

1 2

*dt dt*

ww

2 31 1 3 2

3 12 2 1 3

 t

ˆˆˆ a

 w

() a

 t

 t

 t

> 1 *I*3 *∫ t*0

sin*ω<sup>d</sup> t*. As a further restriction, if the external torque does not exist

*t*

() b

() c

In case two values *I*1 and *I*2 are same (*I*<sup>1</sup> = *I*<sup>2</sup> ⇒ *I*), then Equation 9c becomes simple;

sinusoidal torque of frequency *ωd* as *τ*(*t*)=*τ*0cos*ω<sup>d</sup> t* is readily acquired as

⇀ =0), then *ω*3(*t*) should remain constant as *ω*<sup>3</sup> =*ω*<sup>0</sup> and the two Equations 9a-b become as

20 3 01 3 ( ) 0, ( ) 0 *d d I I I I II*

 w

> w w

+ -= + - =

*i e*

111 22 2 333

v v (8)

w w

replace inertia tensor *Iij*

 w

> w

<sup>v</sup> <sup>v</sup> <sup>v</sup> <sup>v</sup> (7)

⇀ and *<sup>ω</sup>*

b

⇀ × *L*

(underlying basic summar‐

*τ*(*t* ') *dt* '. The solution *ω*3(*t*) for

(10)

⇀ are simply

(9)

is rewritten as follows.

expressed as follows.

ized in Appendix).

*I*3 *dω*<sup>3</sup>

(*τ*

follows.

*ω*3(*t*)=constant +

*τ*0 *I*3*ω<sup>d</sup>*

w

w w

3

302 Geodetic Sciences - Observations, Modeling and Applications

1

 ww

w w

where the principal moments of inertia *Ii*

*j*

=

w

Then Equation 7 is rewritten as

ˆ

*i ij j i i*

w

$$
\alpha\_1 = m\alpha\_0 \cos \Omega t,\\
\alpha\_2 = m\alpha\_0 \sin \Omega t \tag{11}
$$

where *Ω* = *I*<sup>3</sup> − *I <sup>I</sup> <sup>ω</sup>*0. The argument *Ωt* may be added with a phase angle; *Ω<sup>t</sup>* <sup>⇒</sup> *<sup>Ω</sup><sup>t</sup>* <sup>+</sup> *<sup>ϕ</sup>*. Amplitude *mω*0 and phase angle *ϕ* can be specified as initial condition. This rotational motion is often called Eulerian free nutation or wobble. Suppose a body initially rotating along its principal axis is agitated by a small perturbation, then the rotating body will undergo a wobbling motion as it rotates.

If the Earth were perfectly rigid, the frequency of its wobbling motion should be specified as *<sup>Ω</sup>* <sup>=</sup> *<sup>C</sup>* <sup>−</sup> *<sup>A</sup> <sup>A</sup> <sup>ω</sup>*0, of which corresponding period is 304.5 sidereal days, *i.e*., 303.6 days. As rigid Earth rotates daily, an observer on it would see that both the angular momentum and velocity vectors slowly encircle the *x*3-axis (Figure 14). The three vectors stay on a same rotating plane (not shown in the Figure), and the angles between them remain unchanged. Due mainly to finite elasticity of the Earth's mantle, the period of real Earth's wobble is increased by roughly twenty percent. The estimate according to recent observation is about 433 ~ 434 days. To honor the one who first observed this motion, it is called Chandler wobble.

**Figure 14.** Two illustrations of wobble (free Eulerian nutation). Precession of the angular velocity and momentum vec‐ tors around the principal axis (*x*3-axis) seen by an observer on the body (left). To an observer in space, the angular momentum vector remains unchanged, while the angular velocity vector and *x*3-axis rotate around it (right).

Denote small variations in inertia tensor, angular velocity, and angular momentum as *ΔIij*, *ω*0*m* ⇀ <sup>=</sup>*ω*0(*m*1, *<sup>m</sup>*2, *<sup>m</sup>*3), and *<sup>h</sup>* ⇀ =(*h*1, *<sup>h</sup>*2, *<sup>h</sup>*3), then the terms for angular momentum and Equation 7 are rewritten as follows.

First, inertia tensor *Iij* is written as

$$
\begin{bmatrix} I\_{11} & I\_{12} & I\_{13} \\ I\_{21} & I\_{22} & I\_{23} \\ I\_{31} & I\_{32} & I\_{33} \end{bmatrix} \mathbf{e} \begin{bmatrix} A & 0 & 0 \\ 0 & A & 0 \\ 0 & 0 & C \end{bmatrix} \star \begin{bmatrix} \Delta I\_{11} & \Delta I\_{12} & \Delta I\_{13} \\ \Delta I\_{21} & \Delta I\_{22} & \Delta I\_{23} \\ \Delta I\_{31} & \Delta I\_{32} & \Delta I\_{33} \end{bmatrix}
$$

It is noted that *Iij* = *I ji* , *i.e*., inertia tensor is symmetric.

Angular velocity *ω* ⇀ is written as

$$
\begin{bmatrix} \omega\_1\\ \omega\_2\\ \omega\_3 \end{bmatrix} = \omega\_0 \begin{bmatrix} 0\\ 0\\ 1 \end{bmatrix} + \omega\_0 \begin{bmatrix} m\_1\\ m\_2\\ m\_3 \end{bmatrix} = \omega\_0 \begin{bmatrix} m\_1\\ m\_2\\ 1 + m\_3 \end{bmatrix}
$$

Angular momentum *L* ⇀ in tensor form is given as *<sup>L</sup> <sup>i</sup>* <sup>=</sup><sup>∑</sup> *j Iij ω<sup>j</sup>* + *hi* , and its matrix representation is given as

$$
\begin{bmatrix} L\_{1} \\ L\_{2} \\ L\_{3} \end{bmatrix} = \begin{bmatrix} A + \Delta I\_{11} & \Delta I\_{12} & \Delta I\_{13} \\ \Delta I\_{21} & A + \Delta I\_{22} & \Delta I\_{23} \\ \Delta I\_{31} & \Delta I\_{32} & \mathbb{C} + \Delta I\_{33} \end{bmatrix} \omega\_{0} \begin{bmatrix} m\_{1} \\ m\_{2} \\ \mathbf{1} + m\_{3} \end{bmatrix} + \begin{bmatrix} h\_{1} \\ h\_{2} \\ h\_{3} \end{bmatrix} \simeq \omega\_{0} \begin{bmatrix} A m\_{1} + \Delta I\_{13} \\ A m\_{2} + \Delta I\_{23} \\ \mathbf{C} (\mathbf{1} + m\_{3}) + \Delta I\_{33} \end{bmatrix} + \begin{bmatrix} h\_{1} \\ h\_{2} \\ h\_{3} \end{bmatrix}
$$

where first order terms only are retained with neglecting much smaller higher order terms. This equation has been called Liouville's equation.

After a little algebra, three components of the equation; *<sup>d</sup> <sup>L</sup>* ⇀ *dt* <sup>+</sup> *<sup>ω</sup>* ⇀ × *L* ⇀ =0 are found as follows.

$$\begin{aligned} \frac{d\boldsymbol{h}\_1}{dt} + (\frac{d\boldsymbol{\Delta I}\_{13}}{dt} + A\frac{d\boldsymbol{m}\_1}{dt})\omega\_0 - \boldsymbol{h}\_2\omega\_0 + \mathsf{I}(\boldsymbol{C} - A)\boldsymbol{m}\_2 - \boldsymbol{\Delta I}\_{23}\mathbf{J}\omega\_0^2 &= 0 \\\\ \frac{d\boldsymbol{h}\_2}{dt} + (\frac{d\boldsymbol{\Delta I}\_{23}}{dt} + A\frac{d\boldsymbol{m}\_2}{dt})\omega\_0 + \boldsymbol{h}\_1\omega\_0 + \mathsf{I}(A - \mathsf{C})\boldsymbol{m}\_1 + \boldsymbol{\Delta I}\_{13}\mathbf{J}\omega\_0^2 &= 0 \\\\ \frac{d\boldsymbol{h}\_3}{dt} + (\frac{d\boldsymbol{\Delta I}\_{33}}{dt} + C\frac{d\boldsymbol{m}\_3}{dt})\omega\_0 &= 0 \end{aligned}$$

Divide the first and second equations by *ω*0(*C* − *A*) and rearrange terms, then the following two equations are found.

$$\begin{aligned} \frac{1}{\Omega} \frac{dm\_1}{dt} + m\_2 &= \phi\_2 - \frac{1}{a\rho\_0} \frac{d\phi\_1}{dt} & \text{a} \\ \frac{1}{\Omega} \frac{dm\_2}{dt} - m\_1 &= -\phi\_1 - \frac{1}{a\rho\_0} \frac{d\phi\_2}{dt} & \text{b} \\ \frac{dm\_3}{dt} &= -\frac{d}{dt} (\frac{h\_3}{\mathbb{C}a\_0} + \frac{\Delta I\_{33}}{\mathbb{C}}) & \text{c} \end{aligned} \tag{12}$$

where the two excitation functions *ϕ*1 and *ϕ*2 are defined as *ϕ*<sup>1</sup> <sup>=</sup> *<sup>h</sup>*<sup>1</sup> (*C* − *A*)*ω*<sup>0</sup> + *ΔI*<sup>13</sup> *<sup>C</sup>* <sup>−</sup> *<sup>A</sup>* and

$$
\phi\_2 = \frac{h\_2}{(C-A)\omega\_0} + \frac{\Delta I\_{23}}{C-A}.
$$

First, inertia tensor *Iij*

It is noted that *Iij* = *I ji*

Angular velocity *ω*

0 0 1

+ *ω*<sup>0</sup>

Angular momentum *L*

*dΔI*<sup>13</sup> *dt* <sup>+</sup> *<sup>A</sup>*

*dΔI*<sup>23</sup> *dt* <sup>+</sup> *<sup>A</sup>*

*dΔI*<sup>33</sup> *dt* <sup>+</sup> *<sup>C</sup>*

equations are found.

*m*1 *m*2 *m*3

=*ω*<sup>0</sup>

*A* + *ΔI*<sup>11</sup> *ΔI*<sup>12</sup> *ΔI*<sup>13</sup> *ΔI*<sup>21</sup> *A* + *ΔI*<sup>22</sup> *ΔI*<sup>23</sup> *ΔI*<sup>31</sup> *ΔI*<sup>32</sup> *C* + *ΔI*<sup>33</sup>

*dm*<sup>1</sup>

*dm*<sup>2</sup>

*dm*<sup>3</sup> *dt* )*ω*<sup>0</sup> =0

This equation has been called Liouville's equation.

After a little algebra, three components of the equation; *<sup>d</sup> <sup>L</sup>*

*ω*1 *ω*2 *ω*3

=*ω*<sup>0</sup>

is given as

=

*L* <sup>1</sup> *L* <sup>2</sup> *L* <sup>3</sup>

*dh*<sup>1</sup> *dt* <sup>+</sup> (

*dh*<sup>2</sup> *dt* <sup>+</sup> (

*dh*<sup>3</sup> *dt* <sup>+</sup> ( =

304 Geodetic Sciences - Observations, Modeling and Applications

*I*<sup>11</sup> *I*<sup>12</sup> *I*<sup>13</sup> *I*<sup>21</sup> *I*<sup>22</sup> *I*<sup>23</sup> *I*<sup>31</sup> *I*<sup>32</sup> *I*<sup>33</sup> is written as

*ΔI*<sup>11</sup> *ΔI*<sup>12</sup> *ΔI*<sup>13</sup> *ΔI*<sup>21</sup> *ΔI*<sup>22</sup> *ΔI*<sup>23</sup> *ΔI*<sup>31</sup> *ΔI*<sup>32</sup> *ΔI*<sup>33</sup>

, *i.e*., inertia tensor is symmetric.

⇀ in tensor form is given as *<sup>L</sup> <sup>i</sup>* <sup>=</sup><sup>∑</sup>

*ω*0

*dt* )*ω*<sup>0</sup> <sup>−</sup>*h*2*ω*<sup>0</sup> <sup>+</sup> (*<sup>C</sup>* <sup>−</sup> *<sup>A</sup>*)*m*<sup>2</sup> <sup>−</sup>*ΔI*<sup>23</sup> *<sup>ω</sup>*<sup>0</sup>

*dt* )*ω*<sup>0</sup> <sup>+</sup> *<sup>h</sup>*1*ω*<sup>0</sup> <sup>+</sup> (*A*−*C*)*m*<sup>1</sup> <sup>+</sup> *<sup>Δ</sup>I*<sup>13</sup> *<sup>ω</sup>*<sup>0</sup>

*m*1 *m*2 1 + *m*<sup>3</sup>

+

where first order terms only are retained with neglecting much smaller higher order terms.

Divide the first and second equations by *ω*0(*C* − *A*) and rearrange terms, then the following two

0

f

> f

0

( )c

D

1 1 2 2

f w

*dm <sup>d</sup> <sup>m</sup> dt dt*

+ =-

1 1

W

2 2 1 1

f w

*dm <sup>d</sup> <sup>m</sup> dt dt*

1 1 <sup>b</sup>

3 3 33 0

w

*dm h I d dt dt C C*

=- +


*h*1 *h*2 *h*3 *j Iij ω<sup>j</sup>* + *hi*

≃*ω*<sup>0</sup>

<sup>2</sup> =0

<sup>2</sup> =0

a

⇀ *dt* <sup>+</sup> *<sup>ω</sup>* ⇀ × *L*

*Am*<sup>1</sup> + *ΔI*<sup>13</sup> *Am*<sup>2</sup> + *ΔI*<sup>23</sup> *C*(1 + *m*3) + *ΔI*<sup>33</sup>

, and its matrix representation

*h*1 *h*2 *h*3

⇀ =0 are found as follows.

(12)

+

+

⇀ is written as

*m*1 *m*2 1 + *m*<sup>3</sup> The third equation can be rewritten as in equation (12c)

In case with *h*<sup>3</sup> =0, this is rewritten as *m*<sup>3</sup> = − *ΔI*<sup>33</sup> *<sup>C</sup>* . The other case with *ΔI*<sup>33</sup> =0 is also obvious; *<sup>m</sup>*<sup>3</sup> <sup>=</sup> <sup>−</sup> *<sup>h</sup>*<sup>3</sup> *Cω*<sup>0</sup> . Suppose a periodic variation in *h*3 or *ΔI*33exists, then corresponding variation of same periodicity in *m*3 should exist in an amount divided by *Cω*0 or *C*. Equation 12c can be related with changes in LOD and UT1 as follows.

$$-\frac{d}{dt}(\text{LIT1} - \text{LITC}) = \frac{\Delta LOD}{LOD} = -m\_3 = \frac{h\_3}{\text{Co}\_0} + \frac{\Delta I\_{33}}{\text{C}}\tag{13}$$

If the two excitation functions *ϕ*1 and *ϕ*2 do not exist, Equations 12a-b are the same as Equation 10, *i.e*., the body will show pure wobble. In case there exists any periodic excitation, then a perturbation in polar motion of the same periodicity should be accompanied. This can be conveniently shown by complex notation. Write *mc* =*m*<sup>1</sup> + *im*2 and *ϕ<sup>c</sup>* =*ϕ*<sup>1</sup> + *iϕ*2, then we can rewrite Equation 12a-b simply as

$$\frac{1}{\Omega} \frac{dm\_c}{dt} - im\_c = -i\phi\_c - \frac{1}{\alpha\_0} \frac{d\phi\_c}{dt} \tag{14}$$

Write a periodic excitation as *ϕ<sup>c</sup>* ∝*e iω<sup>d</sup> t* , then its corresponding solution for *mc*follows as

$$
\delta m\_c = -\frac{a\rho\_0 + a\rho\_d}{a\rho\_0} \frac{\Omega}{a\rho\_d - \Omega} \phi\_c \tag{15}
$$

### **5. Variation in LOD and polar motion**

Length of day (LOD), which is meant by the length of a solar day, is about 86400 second and slightly variable. The length of a sidereal day is about 86164 second. Besides the secular increase described in a former section, there exist periodic variations in LOD. Even before atomic clocks, precise quartz clocks provided measurement of seasonal perturbation in star transit time; being behind in spring and ahead in late summer by 20-30 millisecond, which should be accumula‐ tion of LOD variation of the same periodicity. Along with the seasonal perturbation, fortnightly and monthly variations in LOD exist. Amplitudes of LOD variations of these different periodic components are in the order of one millisecond. Some amounts of these periodic perturbations are associated with body/ocean tides in the Earth. However, there is strong atmospheric effect on LOD variation. There also exist large quasiperiodic variations of much longer period range, called decadal fluctuation.

Body tidal variation of LOD is briefly described below. Even though Earth's angular momen‐ tum remains unchanged (ignoring the small secular deceleration), angular velocity of the Earth will change when its moment of inertia along the rotational axis changes. This situation is a slight modification from Equation 13, as follows.

$$\frac{\Delta LOD}{LOD} = -m\_3 = \frac{\Delta I\_{33}}{C}$$

where *ΔLOD* is the excessive amount of LOD, and related are *z*-components of Earth's spin angular velocity and inertia tensor; *ω*<sup>3</sup> =(1 + *m*3)*ω*0 and *I*<sup>33</sup> =*C* + *ΔI*33 respectively. Accumulation of *ΔLOD* lead to delay of UT1, and this relation can be expressed as

$$
\Delta UT \, 1 - \text{LTC} = -\sum\_{d \neq s} \Delta LOD = -\int \frac{\Delta LOD}{LOD} dt
$$

As *ΔLOD* due to tidal deformation can be expected, tidal periodicities exist in the *ΔLOD* spectrum. The *ΔLOD* associated with the change in the moment of inertia due to the zonal body tide was early studied by Woolard [39] and later by Yoder and Merriam [40-41]. Follow‐ ing them, the expression for *ΔLOD* due to elastic body tidal deformation is written as

$$
\Delta LOD = -86400 \text{s} \times \frac{k\_2 M R\_e \,^5 \text{Ar}}{\text{Cr}^4} \sqrt{\frac{5}{\pi}} (1 - 3 \sin^2 \delta) \tag{16}
$$

where *k*2 is tidal Love number, *M* is the mass of tide raising body (Moon or Sun), *Re* is the mean radius of the Earth, *r* is the distance to the body from the Earth, *δ* is the declination of the body, and *Δr* is the deviation of *r*from its average value. As will be described below, atmospheric and oceanic contributions to *ΔLOD* exist in large amount, so that calculated body tidal *ΔLOD* alone does not suffice to explain observational data.

In Figure 15, the excessive LOD time series (IERS EOP 08 C04 – simply C04) since 1981 is shown (top). The curve of variation with period longer than 2000 days is superposed with the data (top). Below three different period band components are separately shown (period between 500 and 2000 days, between 100 and 500 days, and less than 100 days). The amplitude of interannual variation is about 0.3 millisecond or less. Decomposition of the time series of Figure 15, was done by simple spectral windowing. It is interesting to note that the present overall trend of LOD shown in this figure is decreasing, which is the reverse of secular increase of LOD due to tidal dissipation. This is ascribed to certain geophysical processes in the Earth's core and mantle, such as geodynamo. Recent fast retreat of glaciers might be related as well. The *ΔLOD* time series show clear annual, semiannual, monthly signals and other minor ones distributed over wide spectral range.

and monthly variations in LOD exist. Amplitudes of LOD variations of these different periodic components are in the order of one millisecond. Some amounts of these periodic perturbations are associated with body/ocean tides in the Earth. However, there is strong atmospheric effect on LOD variation. There also exist large quasiperiodic variations of much longer period range,

Body tidal variation of LOD is briefly described below. Even though Earth's angular momen‐ tum remains unchanged (ignoring the small secular deceleration), angular velocity of the Earth will change when its moment of inertia along the rotational axis changes. This situation is a

where *ΔLOD* is the excessive amount of LOD, and related are *z*-components of Earth's spin angular velocity and inertia tensor; *ω*<sup>3</sup> =(1 + *m*3)*ω*0 and *I*<sup>33</sup> =*C* + *ΔI*33 respectively. Accumulation

As *ΔLOD* due to tidal deformation can be expected, tidal periodicities exist in the *ΔLOD* spectrum. The *ΔLOD* associated with the change in the moment of inertia due to the zonal body tide was early studied by Woolard [39] and later by Yoder and Merriam [40-41]. Follow‐

5

where *k*2 is tidal Love number, *M* is the mass of tide raising body (Moon or Sun), *Re* is the mean radius of the Earth, *r* is the distance to the body from the Earth, *δ* is the declination of the body, and *Δr* is the deviation of *r*from its average value. As will be described below, atmospheric and oceanic contributions to *ΔLOD* exist in large amount, so that calculated body

In Figure 15, the excessive LOD time series (IERS EOP 08 C04 – simply C04) since 1981 is shown (top). The curve of variation with period longer than 2000 days is superposed with the data (top). Below three different period band components are separately shown (period between 500 and 2000 days, between 100 and 500 days, and less than 100 days). The amplitude of interannual variation is about 0.3 millisecond or less. Decomposition of the time series of Figure 15, was done by simple spectral windowing. It is interesting to note that the present overall trend of LOD shown in this figure is decreasing, which is the reverse of secular increase of LOD due to tidal dissipation. This is ascribed to certain geophysical processes in the Earth's core and mantle, such as geodynamo. Recent fast retreat of glaciers might be related as well.

D

4 <sup>5</sup> <sup>86400</sup> (1 3sin ) *<sup>e</sup> k MR r LOD s Cr*

2 2

p

D =- ´ - (16)

d

ing them, the expression for *ΔLOD* due to elastic body tidal deformation is written as

called decadal fluctuation.

*ΔLOD*

*LOD* <sup>=</sup> <sup>−</sup>*m*<sup>3</sup> <sup>=</sup>

*UT* 1−*UTC* = − ∑

slight modification from Equation 13, as follows.

306 Geodetic Sciences - Observations, Modeling and Applications

*<sup>Δ</sup>LOD* <sup>=</sup> <sup>−</sup>*<sup>∫</sup> <sup>Δ</sup>LOD*

of *ΔLOD* lead to delay of UT1, and this relation can be expressed as

*LOD dt*

tidal *ΔLOD* alone does not suffice to explain observational data.

*ΔI*<sup>33</sup> *C*

*days*

**Figure 15.** Variation of LOD since 1981. From the top, excessive LOD time series with its decadal trend, components of period between 500 and 2000 days, between 100 and 500 days, and less than 100 days. The dotted line is zero level. Three lower graphs are shifted for convenience, but all are in same vertical scale.

Fourier spectra of the LOD time series are shown in Figure 16. The two graphs are equivalent but differ in representation only by the axis (frequency/period). To identify the long period component accurately, calculation for these spectra was done through integration of the time series multiplied by sine/cosine sinusoid of each frequency not by using fast Fourier transform. Four largest peaks are labeled as f1, f2, f3, and f4. The two peaks f2 and f3 are annual and semiannual ones. The period of f4 is 13.7 days. As shown in lower spectrum, the period of peak f1 is spread over a wide range. Spectral peak of f1 is split by a few minor peaks. The main period of f1 peak is about 6905 days. Including the second largest peak at 4650 days, total width of half amplitude is about 5290 days between 4150 and 9440 days. This broad peak of f1 is generally believed to be associated with geomagnetic field generation in the Earth's core and 18.6 year precession of lunar orbital plane. Peak sequence in power is follows; f1 (decadal), f4 (fortnightly), f2 (annual), f3 (semiannual). Two spectral peaks of which periods are 27.6 and 9.13 days are noticed as next large ones. There exist other smaller peaks, including 14.8 day period one.

**Figure 16.** Fourier amplitude spectra of the LOD time series shown in Figure 15.

In Figure 17, two short period band components of *ΔLOD* time series of Figure 15 are shown again for three year time span between 2000 and 2003 with the body tidal *ΔLOD* calculated by using Equation 16. Both the data and calculated body tidal *ΔLOD* in Figure 17 well show annual, semiannual, monthly and fortnightly periodic components. However, there are certain differences between the two. First, the fortnightly periodicity is not quite certain in the calculated body tidal *ΔLOD*, and neither is the semiannual one. Moreover, annual and semiannual signals are much stronger in the data than the calculated ones. Evidently this discrepancy should be ascribed to other effects than body tidal perturbation.

**Figure 17.** Observed LOD variation of period less than 500 days (above) and calculated LOD variation due to body tidal perturbation (below).

The effects of global zonal wind pattern fluctuation and ocean tidal angular momentum variation on UT and LOD had been suspected [1-2]. These effects were confirmed later [42-47], and the zonal wind pattern change was found to do the dominant role. Hydrological mass transport is also identified to have effect on annual and semiannual variation. Ocean tidal effect is found to be the main cause for diurnal and semidiurnal variation in LOD. These investiga‐ tions became possible due to recent accurate ocean tide modeling and atmospheric angular momentum data. In Figure 18, both observed data and modeled LOD time series are shown after Chen *et al.* [46].

**Figure 16.** Fourier amplitude spectra of the LOD time series shown in Figure 15.

308 Geodetic Sciences - Observations, Modeling and Applications

discrepancy should be ascribed to other effects than body tidal perturbation.

In Figure 17, two short period band components of *ΔLOD* time series of Figure 15 are shown again for three year time span between 2000 and 2003 with the body tidal *ΔLOD* calculated by using Equation 16. Both the data and calculated body tidal *ΔLOD* in Figure 17 well show annual, semiannual, monthly and fortnightly periodic components. However, there are certain differences between the two. First, the fortnightly periodicity is not quite certain in the calculated body tidal *ΔLOD*, and neither is the semiannual one. Moreover, annual and semiannual signals are much stronger in the data than the calculated ones. Evidently this

**Figure 17.** Observed LOD variation of period less than 500 days (above) and calculated LOD variation due to body

tidal perturbation (below).

**Figure 18.** Observation and model for LOD; Space95 time series and calculated LOD excitation by atmospheric flow and pressure (after Chen *et al*. [46])

Position of Earth's rotational pole on the Earth's surface is slowly but ceaselessly changing. Average position of pole between 1900 and 1905 was taken as the reference point of geodetic latitude and longitude. This position was formerly called Conventional International Origin, and now is renamed as Reference Pole. Polar motion is referred to the offset of Earth's rotational pole with respect to the Reference Pole. The main feature of polar motion during the past century can be summarized as; (i) slow drift along the direction between West 70 and 80°, (ii) Chandler wobble of amplitude about 210 (150 – 280) milliarcsec, and (iii) annual wobble of amplitude about 120 (90-150) milliarcsec. There also exist various components of smaller amplitude, among which is semiannual component.

In the former section, free and forced wobble of rotating rigid Earth was considered. It is rather elaborate to calculate Chandler wobble period for more realistic Earth model. Chandler wobble frequency of elastic and oceanless Earth was acquired as *<sup>Ω</sup>* <sup>=</sup> *<sup>C</sup>* <sup>−</sup> *<sup>A</sup>*−*<sup>D</sup> <sup>A</sup>* <sup>+</sup> *<sup>D</sup> <sup>ω</sup>*0, where *D* =*kRe* <sup>5</sup>*ω*<sup>0</sup> <sup>2</sup> / 3*G* [48]. Corresponding period is 446.2 days. Assuming equilibrium pole tide, the estimate was adjusted into 425.5 day [48]. Its recent theoretical estimate is 423.5 days [49]. Among many observational estimates for Chandler period and Q value, 433.0 days and 179 are the most representative ones [50].

Three other modes of Earth rotation exist due to the core of the Earth as followings; free core nutation, free inner core nutation, and inner core wobble. For an observer on the Earth's surface, free core nutation and free inner core nutation are retrograde motions having their periods approximately one day, therefore these two are called nearly diurnal free wobbles [50-52]. For an observer in space, the periods of these two are much longer; about 430 and 1000 days. Inner core wobble period has been estimated as 900-2500 days [53-54], however its existence has not yet been reported [55].

Formerly, perturbation in the Earth's angular velocity and pole offset were considered to be the same with only difference in their directional notation, *i.e*., (*m*1, *m*2)=(*xp*, − *yp*). In fact, the relation (*m*1, *m*2)≃(*xp*, − *yp*) approximately holds for long period components. Their exact relation is the following; *mc* <sup>=</sup> *pc* <sup>−</sup> *<sup>i</sup> ω*0 *d pc dt* , where *mc* <sup>=</sup>*m*<sup>1</sup> <sup>+</sup> *im*<sup>2</sup> and *pc* <sup>=</sup> *xp* <sup>−</sup>*<sup>i</sup> yp*. Celestial Inter‐ mediate Pole (CIP), denoted by *x*p and *y*p, is the one to be used for TRF-CRF conversion. For more detail, see [5] or related references cited therein. It is noted that, for dynamical effects such as pole tide, values of (*m*1, *m*2) should be referenced to *x*3-principal axis of the Earth. Since both the *x*p and *y*<sup>p</sup> are small, the transformation matrix for polar motion necessary for CRF to TRF transformation can be written as *Rpm*≃*R*2(− *xp*)*R*1(− *yp*)≃*R*1(− *yp*)*R*2(− *xp*) as first order approximation (see Appendix for more detail).

Recent polar motion is illustrated in Figure 19. The graphs are based on the polar motion dataset of IERS EOP 08 C04.

**Figure 19.** Polar motion; two components *x*p and *y*p of polar motion since 1962. Least square fit lines are drawn for whole time span and latter time span since 1981 (left). Two each (long/short) locus of pole position with respect to Reference Pole are drawn (right).

From the lines fitted with least square error in Figure 19, the linear trend of recent polar motion is read as 8.1 cm/yr along W 59° for short time span since 1981 and 12 cm/yr along W 64° for long time span since 1962. Some former estimates of the linear trend in polar motion were 10.3 cm/yr along W 75° [56], 10.9 cm/yr along W 79°[57], and 10.3 cm/yr along W 76°[58]. Therefore, recent pole drift is comparatively slower and tilted by several degrees to the East. This might be associated with the recent rapid glacier melting in Greenland.

periods approximately one day, therefore these two are called nearly diurnal free wobbles [50-52]. For an observer in space, the periods of these two are much longer; about 430 and 1000 days. Inner core wobble period has been estimated as 900-2500 days [53-54], however its

Formerly, perturbation in the Earth's angular velocity and pole offset were considered to be the same with only difference in their directional notation, *i.e*., (*m*1, *m*2)=(*xp*, − *yp*). In fact, the relation (*m*1, *m*2)≃(*xp*, − *yp*) approximately holds for long period components. Their exact

mediate Pole (CIP), denoted by *x*p and *y*p, is the one to be used for TRF-CRF conversion. For more detail, see [5] or related references cited therein. It is noted that, for dynamical effects such as pole tide, values of (*m*1, *m*2) should be referenced to *x*3-principal axis of the Earth. Since both the *x*p and *y*<sup>p</sup> are small, the transformation matrix for polar motion necessary for CRF to TRF transformation can be written as *Rpm*≃*R*2(− *xp*)*R*1(− *yp*)≃*R*1(− *yp*)*R*2(− *xp*) as first order

Recent polar motion is illustrated in Figure 19. The graphs are based on the polar motion

**Figure 19.** Polar motion; two components *x*p and *y*p of polar motion since 1962. Least square fit lines are drawn for whole time span and latter time span since 1981 (left). Two each (long/short) locus of pole position with respect to

*dt* , where *mc* <sup>=</sup>*m*<sup>1</sup> <sup>+</sup> *im*<sup>2</sup> and *pc* <sup>=</sup> *xp* <sup>−</sup>*<sup>i</sup> yp*. Celestial Inter‐

*ω*0

*d pc*

existence has not yet been reported [55].

310 Geodetic Sciences - Observations, Modeling and Applications

relation is the following; *mc* <sup>=</sup> *pc* <sup>−</sup> *<sup>i</sup>*

dataset of IERS EOP 08 C04.

Reference Pole are drawn (right).

approximation (see Appendix for more detail).

Two main components of polar motion are Chandler and annual wobbles. Other known minor components are semi-annual, semi-Chandler, Markowitz, *etc*. There have been investigations to identify different components of polar motion and their characteristics (see, for example, [59-61]). For the polar motion time series as shown in Figure 19, two spectra were acquired by using fast Fourier transform and maximum entropy method, and are illustrated in Figure 20. Minor peaks p1, p2, and p3 in the spectrum are previously known components; semi-annual, semi-Chandler, and 300-day period ones. The 300-day period component is regarded to be associated with atmospheric phenomena [60]. Existence of 490-day period component (peak p4) in polar motion was suggested after numerical experiment and evidence of same perio‐ dicity in other geophysical phenomena [62].

**Figure 20.** Fourier and Maximum likelihood amplitude spectra of recent polar motion shown in Figure 19.

In Figure 21, different period band polar motion components are shown. The each separate time series are: *x*p with long time trend, long period, Chandler wobble, annual wobble, and short period components, which were acquired by spectral windowings.

**Figure 21.** Polar motion and its components of different period ranges. From above, original total data with its trend, long period, Chandler, annual, short period components. (*x*p of C04 between 1981 and 2012 summer)

Variability of the Chandler wobble period has been assumed, because of observationally derived period of Chandlerian motion showed such instability. However, it was asserted that the period of Chandler wobble should be a constant, which is determined by the whole mechanical structure of the Earth [63]. The apparent variation of Chandler period is caused by variable excitation.

While Equation 14 relates excitation function and polar motion of rigid Earth, similar equation for the real Earth is acquired by replacing the Chandler frequency with the following one; *<sup>Ω</sup>* <sup>=</sup> *<sup>ω</sup>*<sup>0</sup> <sup>433</sup> (1 <sup>+</sup> *i* <sup>2</sup>*<sup>Q</sup>* ). Considering the relation *mc* <sup>=</sup> *pc* <sup>−</sup> *<sup>i</sup> ω*0 *d pc dt* together, the excitation function for the real Earth can be conveniently expressed in frequency domain as follows.

$$\Phi(\alpha) = \frac{\alpha\_0}{\alpha\_0 + \alpha} \frac{\Omega - \alpha}{\Omega} M(\alpha) = \frac{\Omega - \alpha}{\Omega} P(\alpha) \tag{17}$$

where *Φ*(*ω*), *M* (*ω*), and *P*(*ω*) are the Fourier transform pairs of *ϕc*(*t*), *mc*(*t*) and *pc*(*t*) respec‐ tively. In Figure 22, both real and imaginary parts of the calculated excitation function *ϕc*(*t*) are illustrated (Q value is taken as 180) together with polar motion.

In Figure 21, different period band polar motion components are shown. The each separate time series are: *x*p with long time trend, long period, Chandler wobble, annual wobble, and

**Figure 21.** Polar motion and its components of different period ranges. From above, original total data with its trend,

Variability of the Chandler wobble period has been assumed, because of observationally derived period of Chandlerian motion showed such instability. However, it was asserted that the period of Chandler wobble should be a constant, which is determined by the whole mechanical structure of the Earth [63]. The apparent variation of Chandler period is caused by

While Equation 14 relates excitation function and polar motion of rigid Earth, similar equation for the real Earth is acquired by replacing the Chandler frequency with the following one;

*ω*0

*d pc*

 w

W- W- F = <sup>=</sup> +W W (17)

ww

*dt* together, the excitation function for

long period, Chandler, annual, short period components. (*x*p of C04 between 1981 and 2012 summer)

<sup>2</sup>*<sup>Q</sup>* ). Considering the relation *mc* <sup>=</sup> *pc* <sup>−</sup> *<sup>i</sup>*

0 0

w

are illustrated (Q value is taken as 180) together with polar motion.

w w

w

the real Earth can be conveniently expressed in frequency domain as follows.

( ) *M P* () ()

where *Φ*(*ω*), *M* (*ω*), and *P*(*ω*) are the Fourier transform pairs of *ϕc*(*t*), *mc*(*t*) and *pc*(*t*) respec‐ tively. In Figure 22, both real and imaginary parts of the calculated excitation function *ϕc*(*t*)

w

variable excitation.

*i*

*<sup>Ω</sup>* <sup>=</sup> *<sup>ω</sup>*<sup>0</sup> <sup>433</sup> (1 <sup>+</sup>

short period components, which were acquired by spectral windowings.

312 Geodetic Sciences - Observations, Modeling and Applications

**Figure 22.** Polar motion and excitation function between year 2000 and 2003 [unit: arcsec]. Polar motion compo‐ nents are delineated, and excitation function is calculated by using Equation 17. Two components of excitation func‐ tion are shifted by ±0.1 arcsec for convenience.

Excitation mechanism of Chandler wobble has been investigated for a long time, and, nowa‐ days, fluid sphere forcing at the Earth's surface is regarded to be the main source [64-71]. Annual wobble exists with amplitude slightly smaller than Chandler wobble. Much smaller semi-annual wobble also exists. Major part of annual and semiannual wobble should be due to atmospheric excitation. Eurasian continent, North Atlantic, North Pacific, and southern oceans were found to be large sources of the atmospheric excitation for annual wobble. Explanation for polar motion excitation is sought by simultaneously considering effects from wind, atmospheric pressure, ocean current, ocean bottom pressure. In Table 3, annual and semiannual excitation components are listed after Gross *et al*. [67]. Table 3 is consisted of prograde components only, however, retrograde ones were reported of the same order of magnitude.


**Table 3.** Calculated and observed polar motion excitation of annual, semiannual, and terannual components (prograde components only) after Gross *et al.* [67] (m.a.s. = milliarcsec, Atmo. = Atmospheric, O.B. = Ocean Bottom). Oceanic excitation of periods between daily to seasonal has been recognized. Atmospheric excitation is found more important for LOD, and both oceanic/atmospheric effects are found important for polar motion [72-73]. For subdaily polar motion, ocean effect is known to dominate [72]. Zhou *et al*. found better assessment of the atmospheric excitation by considering the Earth's surface topography [74]. Nowadays, daily and sub-daily variations of LOD and polar motion are observed and modeled in submicrosecond and microarcsec levels; for example, see [75-76]. However, observation and model do not match completely in all spectral range both for LOD and polar motion. This is due to insufficient coverage of observational data for atmospheric/oceanic/hydrologic excitations [77-79]. Jin *et al*. analyzed hydrologic/ oceanic excitation to polar motion by analyzing GRACE data [77, 79]. Better explanation for annual and semiannual LOD variation was found by GRACE+SLR analysis for the Earth's principal moment of inertia [78].

### **6. Time system and coordinate transformation between TRF and CRF**

Universal Time (UT) is the hour angle of the apparent Sun from the Greenwich meridian. Greenwich Sidereal Time is the hour angle of the vernal equinox from the Greenwich meridian. UT1 is corrected of tiny daily oscillation of UT, which is due to the pole offset. Terrestrial Time (TT), also called as Terrestrial Dynamical Time (TDT), is uniform time in Earth based coordi‐ nate frame. International Atomic Time (TAI) is attained by atomic clocks, and has a constant time difference with TT. Coordinated Universal Time (UTC) is based on TAI and maintained close to UT1 within 0.9 second difference by assigning leap seconds. The difference between UT1 from UTC is dUT1.

$$\begin{aligned} \text{TT = TA1 + 32.184s} \\ \text{UTC = TA1 - \text{(sum of leap seconds)}} \\ \text{UT1 = UTC + dUT1} \end{aligned} \tag{18}$$

Two Earth centered coordinate frames in common use are celestial reference frame (CRF) and terrestrial reference frame (TRF). For most civilian purposes, TRF is the one to use, while CRF is convenient in astronomy. Due to the spin rotation of the Earth, transformation between TRF and CRF is needed. As summarized in former sections, not only the simple rotation angle specified as Greenwich sidereal time, but also corrections due to the precession, nutation and polar motion are necessary.

Coordinate transformation from CRF to TRF can be expressed as follows (IAU 2000A).

$$
\begin{bmatrix} \mathbf{x}\_1 \\ \mathbf{x}\_2 \\ \mathbf{x}\_3 \end{bmatrix}\_{Terential} = \mathbf{R}\_{pm} \mathbf{R}\_{spin} \mathbf{R}\_{unit} \mathbf{R}\_{prec} \begin{bmatrix} \mathbf{x}\_1 \\ \mathbf{x}\_2 \\ \mathbf{x}\_3 \end{bmatrix}\_{Colstein} \tag{19}
$$

where *Rprec*, *Rnut*, *Rspin*, and *Rpm* represent the rotation matrices corresponding to the preces‐ sion, nutation, and spin with time passage, and polar motion.

$$\begin{aligned} R\_{proc} &= R\_3(-\varepsilon)R\_2(\Theta)R\_3(-\zeta) \\ R\_{mut} &= R\_1(-\varepsilon - \Delta\varepsilon)R\_3(-\Delta\psi)R\_1(\varepsilon) \\ R\_{spin} &= R\_3(GAS^-) \\ R\_{pm} &= R\_2(-\chi\_p)R\_1(-\chi\_p) \end{aligned}$$

Oceanic excitation of periods between daily to seasonal has been recognized. Atmospheric excitation is found more important for LOD, and both oceanic/atmospheric effects are found important for polar motion [72-73]. For subdaily polar motion, ocean effect is known to dominate [72]. Zhou *et al*. found better assessment of the atmospheric excitation by considering the Earth's surface topography [74]. Nowadays, daily and sub-daily variations of LOD and polar motion are observed and modeled in submicrosecond and microarcsec levels; for example, see [75-76]. However, observation and model do not match completely in all spectral range both for LOD and polar motion. This is due to insufficient coverage of observational data for atmospheric/oceanic/hydrologic excitations [77-79]. Jin *et al*. analyzed hydrologic/ oceanic excitation to polar motion by analyzing GRACE data [77, 79]. Better explanation for annual and semiannual LOD variation was found by GRACE+SLR analysis for the Earth's

**6. Time system and coordinate transformation between TRF and CRF**

Universal Time (UT) is the hour angle of the apparent Sun from the Greenwich meridian. Greenwich Sidereal Time is the hour angle of the vernal equinox from the Greenwich meridian. UT1 is corrected of tiny daily oscillation of UT, which is due to the pole offset. Terrestrial Time (TT), also called as Terrestrial Dynamical Time (TDT), is uniform time in Earth based coordi‐ nate frame. International Atomic Time (TAI) is attained by atomic clocks, and has a constant time difference with TT. Coordinated Universal Time (UTC) is based on TAI and maintained close to UT1 within 0.9 second difference by assigning leap seconds. The difference between

( )

Two Earth centered coordinate frames in common use are celestial reference frame (CRF) and terrestrial reference frame (TRF). For most civilian purposes, TRF is the one to use, while CRF is convenient in astronomy. Due to the spin rotation of the Earth, transformation between TRF and CRF is needed. As summarized in former sections, not only the simple rotation angle specified as Greenwich sidereal time, but also corrections due to the precession, nutation and

(18)

(19)

UTC = TAI – sum of leap seconds

Coordinate transformation from CRF to TRF can be expressed as follows (IAU 2000A).

*pm spin nut prec Terrestrial Celestial*

1 1 2 2 3 3

*x x x RR RR x x x* é ù é ù ê ú ê ú <sup>=</sup> ê ú ê ú ê ú ê ú ë û ë û

TT = TAI + 32.184s

UT1 = UTC + dUT1

principal moment of inertia [78].

314 Geodetic Sciences - Observations, Modeling and Applications

UT1 from UTC is dUT1.

polar motion are necessary.

Explanations for GAST (Greenwich Apparent Sidereal Time) and Time System are given below.

*GAST* =*GMST* + *Δψ*cos*ε*(in radian)

where GMST is the Greenwich Mean Sidereal Time in radian, and *Δψ*cos*ε* is the nutation of right ascension.

*GMST* =*GMST*<sup>0</sup> + *αUT* 1(in second)

GMST0 is given as

*GMST*<sup>0</sup> =24110.54841 <sup>+</sup> 8640184.812866*<sup>T</sup>* <sup>+</sup> 0.093104*<sup>T</sup>* <sup>2</sup> <sup>−</sup>6.2×10−<sup>6</sup> *T* 3(in second)

where *T* is time in Julian century (36525 days) from J2000.0 to 0h UT1 of the day. Julian day conversion from calendar date is preferred to calculate *T* (for the conversion, see [4]). *α*is the factor for sidereal time conversion.

*α* =1.002737909350795 + 5.9006×10−11*T* −5.9×10−15*T* <sup>2</sup>

Brief description for the transformation between TRF and CRF adopted by the IAU 2006 is given here.

$$
\begin{bmatrix} \mathbf{x}\_1 \\ \mathbf{x}\_2 \\ \mathbf{x}\_3 \end{bmatrix}\_{Tereristal} = \mathbf{R}\_{pm} \mathbf{R}\_{spin} \mathbf{R}\_{prec+mut} \begin{bmatrix} \mathbf{x}\_1 \\ \mathbf{x}\_2 \\ \mathbf{x}\_3 \end{bmatrix}\_{Colstein} \tag{20}
$$

where *Rprec*+*nut*, *Rspin*, and *Rpm* are defined as follows.

$$\mathcal{R}\_{proc\text{\textquotedblleft}unut} = \begin{pmatrix} 1 - b\,\text{X} \,\,^2 & -b\,\text{XY} & -\text{X} \\ -b\,\text{XY} & \mathbf{1} - b\,\text{Y} \,\,^2 & -\text{Y} \\ \mathbf{X} & \mathbf{Y} & \mathbf{1} - b\,\text{(}\text{X}\,\,^2 + \text{Y}\,\,^2\text{)} \end{pmatrix}$$

*Rspin* =*R*3(*θ* −*s* + *s*')

$$R\_{pm} = R\_1(-y\_p)R\_2(-x\_p)R\_3(\frac{\mathcal{X}\_p y\_p}{2})$$

For definitions of *θ*, *s*, and *s*', see [7] or other references cited therein.

### **7. Conclusion**

Different kinds of variations in the Earth's spin rotation are classified and explained. With emphasis on basic principles, each aspects of Earth rotation; precession, nutation, secular deceleration, LOD variation, and polar motion are described concisely. Euler equation and its formal solution for rigid Earth are included with their modifications for real Earth. Transfor‐ mation between CRF and TRF is summarized. For convenience, underlying mechanics and mathematics are briefly summarized in Appendix.

### **Appendix: Summary of pre-requisite mechanics and mathematics**

In this appendix, certain basic physical and mathematical concepts needed for understanding the content of this chapter are explained.

#### **A1. Vector algebra**

A vector in three dimension can be written as *A* <sup>→</sup> <sup>=</sup> *<sup>i</sup>* ^ *A*<sup>1</sup> + *j* ^ *A*<sup>2</sup> + *k* ^ *A*3, where *i* ^ , *j* ^ , and *k* ^ are the unit vectors in each directions. The unit vectors are taken to be cyclic in right handed coordinate system. Magnitude of a vector is defined as *A*= | *A* <sup>→</sup> <sup>|</sup> <sup>=</sup> *<sup>A</sup>*<sup>1</sup> <sup>2</sup> <sup>+</sup> *<sup>A</sup>*<sup>2</sup> <sup>2</sup> <sup>+</sup> *<sup>A</sup>*<sup>3</sup> 2 . Inner product of two vectors is defined as *A* <sup>→</sup> · *<sup>B</sup>* <sup>→</sup> <sup>=</sup> *<sup>A</sup>*1*B*<sup>1</sup> <sup>+</sup> *<sup>A</sup>*2*B*<sup>2</sup> <sup>+</sup> *<sup>A</sup>*3*B*3. Cross product of two vectors is defined as *A* <sup>→</sup> <sup>×</sup> *<sup>B</sup>* <sup>→</sup> <sup>=</sup> *<sup>i</sup>* ^ (*A*2*B*<sup>3</sup> − *A*3*B*2) + *j* ^ (*A*3*B*<sup>1</sup> − *A*1*B*3) + *k* ^ (*A*1*B*<sup>2</sup> − *A*2*B*1). The magnitudes of the two products are | *A* <sup>→</sup> · *<sup>B</sup>* <sup>→</sup> <sup>|</sup> <sup>=</sup> <sup>|</sup> *<sup>A</sup>* ⇀ | | *<sup>B</sup>* ⇀ |cos*θ* and | *<sup>A</sup>* <sup>→</sup> <sup>×</sup> *<sup>B</sup>* <sup>→</sup> <sup>|</sup> <sup>=</sup> <sup>|</sup> *<sup>A</sup>* ⇀ | | *<sup>B</sup>* ⇀ |sin*θ*, where *θ* is the angle between *<sup>A</sup>* ⇀ and *B* ⇀ .

One immediately following fact is that cross product of a vector with itself or any other parallel ones vanish; *A* <sup>→</sup> <sup>×</sup> *<sup>A</sup>* ⇀ =0. A triple cross product can be composed as; *A* <sup>→</sup> ×(*<sup>B</sup>* <sup>→</sup> <sup>×</sup>*<sup>C</sup>* <sup>→</sup> )=(*<sup>A</sup>* <sup>→</sup> · *<sup>C</sup>* <sup>→</sup> )*<sup>B</sup>* <sup>→</sup> <sup>−</sup>(*<sup>A</sup>* <sup>→</sup> · *<sup>B</sup>* <sup>→</sup> )*<sup>C</sup>* → .

### **A2. Harmonic oscillator**

**7. Conclusion**

**A1. Vector algebra**

vectors is defined as *A*

<sup>→</sup> <sup>|</sup> <sup>=</sup> <sup>|</sup> *<sup>A</sup>*

(*A*2*B*<sup>3</sup> − *A*3*B*2) + *j*

⇀ | | *<sup>B</sup>*

*A* <sup>→</sup> <sup>×</sup> *<sup>B</sup>* <sup>→</sup> <sup>=</sup> *<sup>i</sup>* ^

*B* ⇀ .


mathematics are briefly summarized in Appendix.

316 Geodetic Sciences - Observations, Modeling and Applications

the content of this chapter are explained.

A vector in three dimension can be written as *A*

system. Magnitude of a vector is defined as *A*= | *A*

(*A*3*B*<sup>1</sup> − *A*1*B*3) + *k*

<sup>→</sup> <sup>×</sup> *<sup>B</sup>*

^

<sup>→</sup> <sup>|</sup> <sup>=</sup> <sup>|</sup> *<sup>A</sup>*

<sup>→</sup> · *<sup>B</sup>*

⇀ |cos*θ* and | *<sup>A</sup>*

^

Different kinds of variations in the Earth's spin rotation are classified and explained. With emphasis on basic principles, each aspects of Earth rotation; precession, nutation, secular deceleration, LOD variation, and polar motion are described concisely. Euler equation and its formal solution for rigid Earth are included with their modifications for real Earth. Transfor‐ mation between CRF and TRF is summarized. For convenience, underlying mechanics and

**Appendix: Summary of pre-requisite mechanics and mathematics**

In this appendix, certain basic physical and mathematical concepts needed for understanding

<sup>→</sup> <sup>=</sup> *<sup>i</sup>* ^ *A*<sup>1</sup> + *j* ^ *A*<sup>2</sup> + *k* ^

<sup>→</sup> <sup>|</sup> <sup>=</sup> *<sup>A</sup>*<sup>1</sup>

vectors in each directions. The unit vectors are taken to be cyclic in right handed coordinate

⇀ | | *<sup>B</sup>*

*A*3, where *i*

<sup>2</sup> <sup>+</sup> *<sup>A</sup>*<sup>3</sup> 2

(*A*1*B*<sup>2</sup> − *A*2*B*1). The magnitudes of the two products are

⇀ |sin*θ*, where *θ* is the angle between *<sup>A</sup>*

<sup>2</sup> <sup>+</sup> *<sup>A</sup>*<sup>2</sup>

<sup>→</sup> <sup>=</sup> *<sup>A</sup>*1*B*<sup>1</sup> <sup>+</sup> *<sup>A</sup>*2*B*<sup>2</sup> <sup>+</sup> *<sup>A</sup>*3*B*3. Cross product of two vectors is defined as

^ , *j* ^ , and *k* ^

are the unit

⇀ and

. Inner product of two

Oscillation and wave are two phenomena of fundamental importance in science and technol‐ ogy. Wave can be regarded as succession of harmonic oscillation in space. Mechanical oscillation exists in a system, where restoring force is accompanied to a deformation from equilibrium. Theory of one-dimensional harmonic oscillator is summarized below.

#### **Linear restoring force and harmonic motion**

For a displacement *x*of a spring from its equilibrium position, restoring force is exerted by the spring and can be written as *f* = −*kx*. Equation of motion for such system can be written as *m d* 2 *x dt* <sup>2</sup> <sup>=</sup> *<sup>f</sup>* <sup>=</sup> <sup>−</sup>*kx*. The solution for this equation is found as *<sup>x</sup>* <sup>=</sup> *<sup>A</sup>*cos*ω*1*<sup>t</sup>* <sup>+</sup> *<sup>B</sup>*sin*ω*1*t*, where the

angular frequency and period are given as *ω*<sup>1</sup> <sup>=</sup> *<sup>k</sup> m* and *T*<sup>1</sup> =2*π* / *ω*1.

#### **Slightly damped harmonic motion**

Assuming existence of a small viscous dragging force, which varies linearly with the velocity, the equation of motion is given as *m d* 2 *x dt* <sup>2</sup> <sup>=</sup> <sup>−</sup>*kx* <sup>−</sup>*<sup>c</sup> dx dt* . The solution is found as exponentially decaying oscillation.

$$\text{ax}(t) = Ae^{-\gamma t}\cos(\omega\_2 t - \theta), \text{ where } \gamma = \frac{c}{2m} \text{ and } \omega\_2 = \sqrt{\omega\_1^2 - \gamma^2}.$$

#### **Forced harmonic motion**

When a damped harmonic oscillator is driven by external periodic force of angular frequency *ω<sup>d</sup>* , its equation of motion can be written as *m d* 2 *x dt* <sup>2</sup> <sup>=</sup> <sup>−</sup>*kx* <sup>−</sup>*<sup>c</sup> dx dt* <sup>+</sup> *Beiω<sup>d</sup> <sup>t</sup>* , and the solution is found as *x*(*t*)= *Aei*(*ω<sup>d</sup> <sup>t</sup>*−*θ*) , where the amplitude and phase angle are given as *<sup>A</sup>*<sup>=</sup> *<sup>B</sup>* / *<sup>m</sup>* (*ω*<sup>1</sup> <sup>2</sup> <sup>−</sup>*ω<sup>d</sup>* 2 ) <sup>2</sup> + *c* <sup>2</sup> *ωd* <sup>2</sup> / *<sup>m</sup>*<sup>2</sup> 1/2 and tan*<sup>θ</sup>* <sup>=</sup> *cω<sup>d</sup>* / *m ω*1 <sup>2</sup> <sup>−</sup>*ω<sup>d</sup>* 2 .

#### **Euler formula**

The convenient complex notation *ei<sup>α</sup>* =cos*α* + *i*sin*α*, so called Euler formula, can be verified by using three series expansions as follows.

$$\cos\alpha = 1 - \frac{\alpha^2}{2!} + \frac{\alpha^4}{4!} - \frac{\alpha^6}{6!} + \dots \\ \sin\alpha = \alpha - \frac{\alpha^3}{3!} + \frac{\alpha^5}{5!} - \frac{\alpha^7}{7!} + \dots \\ \dots \\ \text{ } e^\beta = 1 + \beta + \frac{\beta^2}{2!} + \frac{\beta^3}{3!} + \frac{\beta^4}{4!} + \dots \\ \dots$$

Substitution *β* =*iα*into the third expression leads to Euler formula.

$$e^{ia} = 1 + i\alpha + \frac{(ia)^2}{2!} + \frac{(ia)^3}{3!} + \frac{(ia)^4}{4!} + \dots \\
= 1 - \frac{\alpha^2}{2!} + \frac{\alpha^4}{4!} + \dots \\
\dots \\
+ i(\alpha - \frac{\alpha^3}{3!} + \frac{\alpha^5}{5!} - \dots) \\
= \cos\alpha + i\sin\alpha$$

#### **A3. Rotational mechanics**

#### **Torque and angular momentum**

Torque *τ* ⇀ exerted by a force *f* ⇀ acting at position *<sup>r</sup>* ⇀ is defined as *τ* ⇀ =*r* ⇀ × *f* ⇀ . Angular momentum *<sup>l</sup>* ⇀ of a single moving body of mass *m* and velocity *v* ⇀ at position *r* ⇀ is defined as *l* ⇀ =*r* ⇀ ×*mv*⇀ . Torque applied to a system results in change in its angular momentum. This relation is stated as *τ* ⇀ <sup>=</sup> *dl*⇀ *dt* , which can be readily verified by time differentiation on *<sup>l</sup>* ⇀ ; *dl*⇀ *dt* <sup>=</sup> *<sup>d</sup> dt* (*<sup>r</sup>* ⇀ ×*mv*⇀ )=*v* ⇀ ×*mv*⇀ + *r* ⇀ ×*m dv*⇀ *dt* <sup>=</sup>*<sup>r</sup>* ⇀ ×*ma*⇀ =*r* ⇀ × *f* ⇀ <sup>=</sup>*<sup>τ</sup>* ⇀ .

#### **Angular velocity**

For a simple rotation occurring in a plane, angular velocity is defined as time rate of rotation angle *<sup>ω</sup>* <sup>=</sup> *<sup>d</sup><sup>θ</sup> dt* . The speed due to rotation is given as *<sup>v</sup>* <sup>=</sup>*ωr*, where *r* is the radius from the rotation axis. It should be kept in mind that rotational motion (velocity) is perpendicular to the rotating axis (angular velocity vector). For a three dimensional rotation, angular velocity vector *ω* ⇀ is defined by its three components; *ω* ⇀ <sup>=</sup>*ω*1*<sup>e</sup>* ^ <sup>1</sup> + *ω*2*e* ^ <sup>2</sup> + *ω*3*e* ^ <sup>3</sup>. Velocity *v* ⇀ of any point *r* ⇀ in the rotating body is given as *v* ⇀ =*ω* ⇀ ×*r* ⇀ .

#### **Inertia tensor**

Consider the angular momentum of a rotating rigid body. Write mass element *dm*at position *r* ⇀ , then its linear momentum is given as *v* ⇀ *dm*=*ω* ⇀ ×*r* ⇀ *dm*. Total angular momentum is expressed as *L* ⇀ <sup>=</sup>*∫<sup>r</sup>* ⇀ ×*v* ⇀ *dm*=*∫r* ⇀ ×(*ω* ⇀ ×*r* ⇀ )*dm*=*∫*(*r* <sup>2</sup> *ω* ⇀ −(*r* ⇀ · *ω* ⇀ )*r* ⇀ )*dm*. Three components of the total angular momentum are given as follows.

$$L\_{-1} = \left[ (r^2 \omega\_1 - (\mathbf{x}\_1 \omega\_1 + \mathbf{x}\_2 \omega\_2 + \mathbf{x}\_3 \omega\_3) \mathbf{x}\_1) dm \right]$$

$$L\_{-2} = \left[ (r^2 \omega\_2 - (\mathbf{x}\_1 \omega\_1 + \mathbf{x}\_2 \omega\_2 + \mathbf{x}\_3 \omega\_3) \mathbf{x}\_2) dm \right]$$

$$L\_{-3} = \left[ (r^2 \omega\_3 - (\mathbf{x}\_1 \omega\_1 + \mathbf{x}\_2 \omega\_2 + \mathbf{x}\_3 \omega\_3) \mathbf{x}\_3) dm \right]$$

These can be rewritten as follows.

$$L\_{1} = I\_{11}\omega\_{1} + I\_{12}\omega\_{2} + I\_{13}\omega\_{3} \quad L\_{2} = I\_{21}\omega\_{1} + I\_{22}\omega\_{2} + I\_{23}\omega\_{3} \quad L\_{3} = I\_{31}\omega\_{1} + I\_{32}\omega\_{2} + I\_{33}\omega\_{3}$$

where components of inertia tensor *Iij* are defined as elements of the following matrix.

$$\begin{bmatrix} I\_{ij} \end{bmatrix} = \begin{bmatrix} I\_{11} & I\_{12} & I\_{13} \\ I\_{21} & I\_{22} & I\_{23} \\ I\_{31} & I\_{32} & I\_{33} \end{bmatrix} = \begin{bmatrix} \int (\mathbf{x}\_{2}^{2} + \mathbf{x}\_{3}^{2}) dm & -\int \mathbf{x}\_{1}\mathbf{x}\_{2} dm & -\int \mathbf{x}\_{1}\mathbf{x}\_{3} dm \\\\ -\int \mathbf{x}\_{2}\mathbf{x}\_{1} dm & \int (\mathbf{x}\_{3}^{2} + \mathbf{x}\_{1}^{2}) dm & -\int \mathbf{x}\_{2}\mathbf{x}\_{3} dm \\\\ -\int \mathbf{x}\_{3}\mathbf{x}\_{1} dm & -\int \mathbf{x}\_{3}\mathbf{x}\_{2} dm & \int (\mathbf{x}\_{1}^{2} + \mathbf{x}\_{2}^{2}) dm \end{bmatrix}$$

#### **Rotational kinetic energy**

Kinetic energy of any moving body is defined as *<sup>T</sup>* <sup>=</sup>*<sup>∫</sup>* <sup>1</sup> <sup>2</sup> *<sup>v</sup>* <sup>2</sup> *dm*. The kinetic energy *T* of a rotating body is given as; *<sup>∫</sup>* <sup>1</sup> <sup>2</sup> *<sup>v</sup>* <sup>2</sup> *dm*<sup>=</sup> <sup>1</sup> 2 *∫v* ⇀ · (*ω* ⇀ ×*r* ⇀ )*dm*<sup>=</sup> <sup>1</sup> 2 *ω* ⇀ ·*∫*(*<sup>r</sup>* ⇀ ×*v* ⇀ )*dm*<sup>=</sup> <sup>1</sup> 2 *ω* ⇀ · *L* ⇀. Therefore, *<sup>T</sup>* <sup>=</sup>*<sup>∫</sup>* <sup>1</sup> <sup>2</sup> *<sup>v</sup>* <sup>2</sup> *dm* = 1 2 *ω* ⇀ · *L* ⇀ or *<sup>T</sup>* <sup>=</sup> <sup>1</sup> 2 ∑ *i*=1 3 ∑ *j*=1 3 *Iijω<sup>i</sup> ωj* .

### **Principal axis**

*<sup>e</sup>i<sup>α</sup>* = 1+ *<sup>i</sup><sup>α</sup>* <sup>+</sup> (*iα*)

Torque *τ*

*τ* ⇀ <sup>=</sup> *dl*⇀

*dl*⇀ *dt* <sup>=</sup> *<sup>d</sup> dt* (*<sup>r</sup>*

2 2! <sup>+</sup> (*iα*)

**Torque and angular momentum**

⇀ exerted by a force *f*

⇀ ×*mv*⇀ + *r*

**A3. Rotational mechanics**

⇀ ×*mv*⇀ )=*v*

defined by its three components; *ω*

⇀ =*ω* ⇀ ×*r* ⇀ .

⇀ , then its linear momentum is given as *v*

⇀ ×(*ω* ⇀ ×*r*

*ω*<sup>1</sup> −(*x*1*ω*<sup>1</sup> + *x*2*ω*<sup>2</sup> + *x*3*ω*3)*x*1)*dm*

*ω*<sup>2</sup> −(*x*1*ω*<sup>1</sup> + *x*2*ω*<sup>2</sup> + *x*3*ω*3)*x*2)*dm*

*ω*<sup>3</sup> −(*x*1*ω*<sup>1</sup> + *x*2*ω*<sup>2</sup> + *x*3*ω*3)*x*3)*dm*

momentum are given as follows.

These can be rewritten as follows.

⇀ )*dm*=*∫*(*r* <sup>2</sup>

*ω* ⇀ −(*r* ⇀ · *ω* ⇀ )*r*

*L* <sup>1</sup> = *I*11*ω*<sup>1</sup> + *I*12*ω*<sup>2</sup> + *I*13*ω*3, *L* <sup>2</sup> = *I*21*ω*<sup>1</sup> + *I*22*ω*<sup>2</sup> + *I*23*ω*3, *L* <sup>3</sup> = *I*31*ω*<sup>1</sup> + *I*32*ω*<sup>2</sup> + *I*33*ω*<sup>3</sup>

where components of inertia tensor *Iij* are defined as elements of the following matrix.

**Angular velocity**

body is given as *v*

**Inertia tensor**

*r*

as *L* ⇀ <sup>=</sup>*∫<sup>r</sup>* ⇀ ×*v* ⇀ *dm*=*∫r*

*<sup>L</sup>* <sup>1</sup> <sup>=</sup>*∫*(*<sup>r</sup>* <sup>2</sup>

*<sup>L</sup>* <sup>2</sup> <sup>=</sup>*∫*(*<sup>r</sup>* <sup>2</sup>

*<sup>L</sup>* <sup>3</sup> <sup>=</sup>*∫*(*<sup>r</sup>* <sup>2</sup>

angle *<sup>ω</sup>* <sup>=</sup> *<sup>d</sup><sup>θ</sup>*

3 3! <sup>+</sup> (*iα*)

318 Geodetic Sciences - Observations, Modeling and Applications

of a single moving body of mass *m* and velocity *v*

⇀ ×*m dv*⇀ *dt* <sup>=</sup>*<sup>r</sup>*

4

4! + . . .=1<sup>−</sup> *<sup>α</sup>* <sup>2</sup>

⇀ acting at position *<sup>r</sup>*

⇀ ×*ma*⇀ =*r*

⇀ <sup>=</sup>*ω*1*<sup>e</sup>* ^ <sup>1</sup> + *ω*2*e* ^ <sup>2</sup> + *ω*3*e* ^

2! <sup>+</sup>

applied to a system results in change in its angular momentum. This relation is stated as

For a simple rotation occurring in a plane, angular velocity is defined as time rate of rotation

axis. It should be kept in mind that rotational motion (velocity) is perpendicular to the rotating axis (angular velocity vector). For a three dimensional rotation, angular velocity vector *ω*

Consider the angular momentum of a rotating rigid body. Write mass element *dm*at position

⇀ *dm*=*ω* ⇀ ×*r*

*dt* . The speed due to rotation is given as *<sup>v</sup>* <sup>=</sup>*ωr*, where *r* is the radius from the rotation

<sup>3</sup>. Velocity *v*

⇀ × *f* ⇀ <sup>=</sup>*<sup>τ</sup>* ⇀ .

*dt* , which can be readily verified by time differentiation on *<sup>l</sup>*

*α* 4

4! + . . . + *<sup>i</sup>*(*<sup>α</sup>* <sup>−</sup> *<sup>α</sup>* <sup>3</sup>

⇀ is defined as *τ*

⇀ at position *r*

3! <sup>+</sup>

⇀ =*r* ⇀ × *f*

⇀ is defined as *l*

⇀ of any point *r*

⇀ *dm*. Total angular momentum is expressed

⇀ )*dm*. Three components of the total angular

*α* 5

5! <sup>−</sup> . . .)=cos*<sup>α</sup>* <sup>+</sup> *<sup>i</sup>*sin*<sup>α</sup>*

⇀ . Angular momentum *<sup>l</sup>*

⇀ ×*mv*⇀ . Torque

⇀ =*r*

⇀

⇀ ;

⇀ is

⇀ in the rotating

For any rigid body, set of three particular orthogonal axes (called principal axes) exists. If body frame of reference coincides with those principal axes, only diagonal component of inertia tensor remain as nonzero.

$$\begin{bmatrix} I\_{i\bar{1}} \end{bmatrix} = \begin{bmatrix} I\_{11} & 0 & 0 \\ 0 & I\_{22} & 0 \\ 0 & 0 & I\_{33} \end{bmatrix}$$

When a body rotates along its principal axis, its motion may be regarded as a simple planar rotation. As an example, suppose a rotation of a body along its principal axis *x*3, then, the angular momentum and rotational kinetic energy are simply given as *L* = *L* <sup>3</sup> = *I*3*ω*3and *<sup>T</sup>* <sup>=</sup> <sup>1</sup> <sup>2</sup> *<sup>I</sup>*3*ω*<sup>3</sup> 2. It is noted that set of principal axis and principal moment of inertia can be deter‐ 3

mined from the following condition; *<sup>L</sup> <sup>i</sup>* <sup>=</sup>∑ *j*=1 *Iijω<sup>j</sup>* = *Iλω<sup>i</sup>* .

#### **Rotating coordinate frame**

Two coordinate frames of common origin, one stationary in space and the other rotating with angular velocity *ω* ⇀ , are described in the figure (left). On the right figure, the unit vector *e* ^ <sup>2</sup> ' at time *t*and *t* + *Δt* are shown. The change *Δe* ^ <sup>2</sup> ' in the unit vector *e* ^ <sup>2</sup> ' is *Δe* ^ 2 '= *e* ^ <sup>2</sup> '(*t* + *Δt*)− *e* ^ <sup>2</sup> '(*t*) ≃ sin*ϕ Δθ n* ^ =sin*<sup>ϕ</sup> ωΔ<sup>t</sup> <sup>n</sup>* ^. Here the direction of *Δ<sup>e</sup>* ^ <sup>2</sup> ' is denoted by a unit vector *n* ^, which is parallel to *<sup>ω</sup>* ⇀ × *e* ^ <sup>2</sup> '. Therefore the time derivative of *e* ^ <sup>2</sup> ' is expressed as *d e*^ 2 ' *dt* <sup>=</sup>*<sup>ω</sup>* ⇀ × *e* ^ <sup>2</sup> '. Likewise time derivatives of *e* ^ <sup>3</sup> ' and *e* ^ <sup>1</sup> ' are found as *d e*^ 3 ' *dt* <sup>=</sup>*<sup>ω</sup>* ⇀ × *e* ^ <sup>3</sup> ' and *d e*^ 1 ' *dt* <sup>=</sup>*<sup>ω</sup>* ⇀ × *e* ^ <sup>1</sup> '. Write the coordinates of a position vector *r* ⇀ in the two reference frame as (*x*1, *x*2, *x*3) and (*x*<sup>1</sup> ', *x*<sup>2</sup> ', *x*<sup>3</sup> '). Then we have *r* ⇀ = *x*1*e* ^ <sup>1</sup> + *x*2*e* ^ <sup>2</sup> + *x*3*e* ^ <sup>3</sup> = *x*<sup>1</sup> '*e* ^ <sup>1</sup> ' + *x*<sup>2</sup> '*e* ^ <sup>2</sup> ' + *x*<sup>3</sup> '*e* ^ 3 '=*r* ⇀ '.

The velocity vector can be written as *d x*<sup>1</sup> *dt <sup>e</sup>* ^ 1 + *d x*<sup>2</sup> *dt <sup>e</sup>* ^ 2 + *d x*<sup>3</sup> *dt <sup>e</sup>* ^ <sup>3</sup> or equivalently as *d x*<sup>1</sup> *dt <sup>e</sup>* ^ <sup>1</sup> ' + *d x*<sup>2</sup> *dt <sup>e</sup>* ^ <sup>2</sup> ' + *d x*<sup>3</sup> *dt <sup>e</sup>* ^ <sup>3</sup> ' + *x*<sup>1</sup> ' *d e*^ 1 ' *dt* <sup>+</sup> *<sup>x</sup>*<sup>2</sup> ' *d e*^ 2 ' *dt* <sup>+</sup> *<sup>x</sup>*<sup>3</sup> ' *d e*^ 3 ' *dt* . Since the time derivatives of unit vectors are *d e*^ *i* ' *dt* <sup>=</sup>*<sup>ω</sup>* ⇀ × *e* ^ *<sup>i</sup>* ', the three terms *x*<sup>1</sup> ' *d e*^ 1 ' *dt* <sup>+</sup> *<sup>x</sup>*<sup>2</sup> ' *d e*^ 2 ' *dt* <sup>+</sup> *<sup>x</sup>*<sup>3</sup> ' *d e*^ 3 ' *dt* can be rewritten as *ω* ⇀ ×(*x*<sup>1</sup> '*<sup>e</sup>* ^ <sup>1</sup> ' + *x*<sup>2</sup> '*e* ^ <sup>2</sup> ' + *x*3*e* ^ <sup>3</sup> ')=*ω* ⇀ ×*r* ⇀ '. Then the relation is expressed as *v* ⇀ =*v* ⇀ ' + *ω* ⇀ ×*r* ⇀ '

This time derivative expression derived for position vector in rotating frame can be extended to other vectors.

*d dt* <sup>|</sup> *inertial* <sup>=</sup> *<sup>d</sup> dt* | *body* + *ω* ⇀ ×

#### **A4. Fourier series and fourier transform**

A periodic function *f* (*t*) of period *T* can be represented by its Fourier series.

$$f'(t) = \frac{a\_0}{2} + \sum\_{n=1}^{\circ} \left( a\_n \cos \frac{2\pi nt}{T} + b\_n \sin \frac{2\pi nt}{T} \right)$$

where the coefficients are defined as follows.

$$a\_n = \frac{2}{T} \int\_{-T/2}^{T/2} f\left(t\right) \cos\frac{2\pi nt}{T} dt \text{ ' and } b\_n = \frac{2}{T} \int\_{-T/2}^{T/2} f\left(t\right) \sin\frac{2\pi nt}{T} dt \text{ ', }$$

Equivalent series for *f* (*t*) can be written as

$$f(t) = \sum\_{n=-\alpha}^{\omega} c\_n \exp(i\frac{2\pi n}{T}t) \text{ with coefficient } c\_n = \frac{1}{T} \int\_{-T/2}^{T/2} f(t') \exp(-i\frac{2\pi n}{T}t')dt'.$$

As can be derived from above expression by taking period *T* as infinity, Fourier transform *F* (*ω*) of arbitrary well-behaved function *f* (*t*) is defined as follows.

$$F(\omega) = \frac{1}{\sqrt{2\pi}} \stackrel{\rightharpoonup}{\int} f(t) \exp(-i\omega t) dt$$

The inverse transform is given as

$$f'(t) = \frac{1}{\sqrt{2\pi}} \int\_{-\omega}^{\omega} F(\omega) \exp(i\omega t) d\omega$$

The velocity vector can be written as

320 Geodetic Sciences - Observations, Modeling and Applications

*d e*^ 1 ' *dt* <sup>+</sup> *<sup>x</sup>*<sup>2</sup> '

*<sup>i</sup>* ', the three terms *x*<sup>1</sup> '

A periodic function *f* (*t*) of period *T* can be represented by its Fourier series.

*<sup>T</sup> ∫* −*T* /2

*T* /2

*f* (*t* ')sin

2*πnt* ' *<sup>T</sup> dt* '

2*πnt <sup>T</sup>* )

*d x*<sup>1</sup> *dt <sup>e</sup>* ^ <sup>1</sup> ' + *d x*<sup>2</sup> *dt <sup>e</sup>* ^ <sup>2</sup> ' + *d x*<sup>3</sup> *dt <sup>e</sup>* ^ <sup>3</sup> ' + *x*<sup>1</sup> '

*ω* ⇀ ×(*x*<sup>1</sup> '*<sup>e</sup>* ^ <sup>1</sup> ' + *x*<sup>2</sup> '*e* ^ <sup>2</sup> ' + *x*3*e* ^ <sup>3</sup> ')=*ω* ⇀ ×*r*

*v* ⇀ =*v* ⇀ ' + *ω* ⇀ ×*r* ⇀ '

*d*

*f* (*t*)= *a*0 <sup>2</sup> <sup>+</sup> ∑ *n*=1 *∞*

*an* <sup>=</sup> <sup>2</sup> *<sup>T</sup> ∫* −*T* /2

*T* /2

vectors are

to other vectors.

*dt* <sup>|</sup> *inertial* <sup>=</sup> *<sup>d</sup>*

*d e*^ *i* ' *dt* <sup>=</sup>*<sup>ω</sup>* ⇀ × *e* ^

*dt* | *body*

(*an*cos

*f* (*t* ')cos

+ *ω* ⇀ ×

**A4. Fourier series and fourier transform**

2*πnt*

where the coefficients are defined as follows.

2*πnt* '

Equivalent series for *f* (*t*) can be written as

*<sup>T</sup>* <sup>+</sup> *bn*sin

*<sup>T</sup> dt* ' and *bn* <sup>=</sup> <sup>2</sup>

*d x*<sup>1</sup> *dt <sup>e</sup>* ^ 1 + *d x*<sup>2</sup> *dt <sup>e</sup>* ^ 2 + *d x*<sup>3</sup> *dt <sup>e</sup>* ^

*d e*^ 1 ' *dt* <sup>+</sup> *<sup>x</sup>*<sup>2</sup> '

⇀ '. Then the relation is expressed as

This time derivative expression derived for position vector in rotating frame can be extended

*d e*^ 3 '

> *d e*^ 2 ' *dt* <sup>+</sup> *<sup>x</sup>*<sup>3</sup> '

*d e*^ 2 ' *dt* <sup>+</sup> *<sup>x</sup>*<sup>3</sup> ' <sup>3</sup> or equivalently as

*dt* can be rewritten as

*dt* . Since the time derivatives of unit

*d e*^ 3 ' Both Fourier series and Fourier transform are widely used. While Fourier transform is more convenient in theoretical development, quite often discrete Fourier transform is used for calculation in practice. For a given sequence *f n* , its discrete Fourier transform is defined as follows.

$$\,\_1F\mathbf{I}\mathbf{J} = \frac{1}{\sqrt{N}} \sum\_{n=0}^{N-1} f\mathbf{I} n \mathbf{J} \exp(-i\frac{2\pi}{N}nk)$$

The inverse transform is given as

$$f[\mathbf{L}n] = \frac{1}{\sqrt{N}} \sum\_{k=0}^{N-1} F[k] \mathbf{exp}(i\frac{2\pi}{N}nk)$$

Both *f n* and *F k* are periodic sequences such that *f n* = *f N* + *n* and *F k* = *F N* + *k* . Fast Fourier transform is an algorithm to evaluate discrete Fourier transform by minimum number of calculations.

#### **A5. Rotation matrix**

Coordinate transformation due to a rotation along one of reference axis by an angle *θ* is represented as follows.

$$R\_1(\theta) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \theta & \sin \theta \\ 0 & -\sin \theta & \cos \theta \end{bmatrix} \\ R\_2(\theta) = \begin{bmatrix} \cos \theta & 0 & -\sin \theta \\ 0 & 1 & 0 \\ \sin \theta & 0 & \cos \theta \end{bmatrix} \\ R\_3(\theta) = \begin{bmatrix} \cos \theta & \sin \theta & 0 \\ -\sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

Any two different axis rotations of finite angles do not commute, unless the angles are infinitesimal. For polar motion, the transformation matrix is given as *Rpm*≃*R*2(− *xp*)*R*1(− *yp*)≃*R*1(− *yp*)*R*2(− *xp*) for first order approximation. Retaining second order terms, the transformation matrix is given as *Rpm*≃*R*2(− *xp*)*R*1(− *yp*)*R*3(− *xp yp* <sup>2</sup> )≃*R*1(<sup>−</sup> *yp*)*R*2(<sup>−</sup> *xp*)*R*3( *xp yp* <sup>2</sup> ).

### **Acknowledgements**

This chapter was written by the author during his stay at Korea Astronomy and Space Science Institute (KASI) under support of KASI Basic Core Technology Development Programs. He thanks Dr. Gross R. and Dr. Jin S. for their kind influence on his works about Earth rotation.

### **Author details**

Sung-Ho Na

Korea Astronomy and Space Science Institute, Korea

### **References**


[12] Mathews PM., Herring TA., Buffett BA. Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth's interior, Journal of Ge‐ ophysical Research 2002; 107, B4, 2068.

**Acknowledgements**

322 Geodetic Sciences - Observations, Modeling and Applications

**Author details**

Korea Astronomy and Space Science Institute, Korea

Cambridge University Press; 1980.

versity Science Books; 1992.

Sung-Ho Na

**References**

1960.

This chapter was written by the author during his stay at Korea Astronomy and Space Science Institute (KASI) under support of KASI Basic Core Technology Development Programs. He thanks Dr. Gross R. and Dr. Jin S. for their kind influence on his works about Earth rotation.

[1] Munk WH., Macdonald GJF. The Rotation of the Earth. Cambridge University Press;

[2] Lambeck K. The Earth's Variable Rotation: Geophysical Causes and Consequences.

[4] Seidelmann PK., editor. Explanatory Supplement to the Astronomical Almanac. Uni‐

[5] Gross RS. Earth Rotation Variations – Long Period. In: Herring T., editor. Geodesy,

[6] Dehant V., Mathews PM. Earth Rotation Variations. In: Herring T., editor. Geodesy,

[7] Petit G., Luzum B., editor. IERS Conventions (2010), IERS Technical Note No. 36, In‐ ternational Earth Rotation and Reference Systems Service (IERS). Verlag des Bunde‐

[8] Danby JMA. Fundamentals of Celestial Mechanics 2nd ed. Willman-Bell Inc; 1988.

[10] Stacey FD., Davis PM. Physics of the Earth 4th ed., Cambridge University Press; 2008.

[11] Capitaine N., Wallace PT., Chapront J. Expressions for IAU 2000 precession quanti‐

[3] Moritz H., Mueller II. Earth Rotation Theory and Observation. Ungar; 1986.

Treatise on Geophysics vol. 3. Elsevier; 2009. p239-294.

Treatise on Geophysics vol. 3. Elsevier; 2009. p295-349.

[9] Melchior P. The Tides of the Planet Earth. Pergamon Press; 1978.

ties, Astronomy and Astrophysics 2003; 412, 567-586.

samts für Kartographie und Geodäsie; 2010.


data and illustration exist in the following article: Morrison LV., Ward CG. An Anal‐ ysis of the Transits of Mercury: 1677-1973, Monthly Notice of Royal Astronomical So‐ ciety 1975; 173, 183-206.


[43] McCarthy DD., Luzum BJ. An analysis of tidal variations in the length of day, Geo‐ physical Journal International 1993; 114, 341-346.

data and illustration exist in the following article: Morrison LV., Ward CG. An Anal‐ ysis of the Transits of Mercury: 1677-1973, Monthly Notice of Royal Astronomical So‐

[28] Morrison LV., Stephenson LV. Historical Values of the Earth's Clock Error DT and the Calculation of Eclipses, Journal for the History of Astronomy 2004; 25, 327-336.

[29] Spencer Jones H. The Rotation of the Earth and the Secular Acceleration of the Sun, Moon and Planets, Monthly Notice of Royal Astronomical Society 1939; 99, 541-548.

[30] Newton R. Secular Accelerations of the Earth and Moon, Science 1969; 166, 825-831.

[31] Rosenberg GD., Runcorn SK., editor. Growth Rhythms and the History of the Earth's

[32] Scrutton CT. Periodic Growth Features in Fossil Organisms and the Length of the Day and Month. In Brosche P, Sündermann J, editors. Tidal Friction and the Earth's

[33] Lambeck K. The Earth's Paleorotation. In Brosche P., Sündermann J., editors. Tidal

[35] Mignard F. Long Time Integration of the Moon's Orbit and references therein. In Bro‐ sche P, Sündermann J, editors. Tidal Friction and the Earth's Rotation II. Springer-

[36] Hansen KS. Secular Effects of Oceanic Tidal Dissipation on the Moon's Orbit and the Earth's Rotation, Review of Geophysics and Space Physics 1982; 20, 457-480.

[37] Webb DJ. On the Reduction in Tidal Dissipation Produced by Increases in the Earth's Rotation Rate and Its Effect on the Long-Term History of the Moon's Orbit. In Bro‐ sche P., Sündermann J., editors. Tidal Friction and the Earth's Rotation II. Springer-

[38] Na S. Tidal Evolution of Lunar Orbit and Earth Rotation, Journal of the Korean As‐

[39] Woolard EW. Inequalities in mean Solar Time from Tidal Variations in the Rotation

[40] Yoder CF., Williams JG., Parke ME. Tidal Variations of Earth Rotation, Journal of Ge‐

[41] Merriam JB. Tidal Terms in Universal Time: Effects of Zonal Winds and Mantle Q,

[34] Goldreich P. History of the lunar orbit, Reviews of Geophysics 1966; 4, 411-439.

Friction and the Earth's Rotation. Springer-Verlag; 1978. p145-153.

ciety 1975; 173, 183-206.

324 Geodetic Sciences - Observations, Modeling and Applications

Rotation. Wiley; 1975.

Verlag; 1982. p67-91.

Verlag; 1982. p210-221.

tronomical Society 2012; 45, 49-57.

ophysical Research 1981; 86, B2, 881-891.

of the Earth, Astronomical Journal 1959; 64, 140-142.

Journal of Geophysical Research 1984; 89, B12, 10109-10114.

[42] Hide R., Dickey JO. Earth's Variable Rotation, Science 1991; 253, 629-637.

Rotation. Springer-Verlag; 1978. p154-196.


[73] Chao BF., Ray RD., Gipson JM., Egbert GD., Ma C. Diurnal/semidiurnal polar motion excited by oceanic tidal angular momentum, Journal of Geophysical research 1996; 101, B9, 20151-20163.

[57] McCarthy DD., Luzum B. Path of the mean rotational pole from 1899 to 1994, Geo‐

[58] Gross RS., Vondrák J. Astrometric and space-geodetic observations of polar wander,

[59] Schuh H., Nagel S., Seitz T. Linear drift and periodic variations observed in long time

[60] Höpfner J. Low-frequency variations, Chandler and annual wobbles of polar motion

[61] Höpfner J. Parameter variability of the observed periodic oscillations of polar motion

[62] Na, S., Cho, J., Baek J., Kwak Y., Yoo SM. Spectral Analysis on Earth's Spin Rotation for the Recent 30 years, Journal of the Korean Physical Society 2012; 61, 152-157.

[63] Jochmann H. Period variations of the Chandler wobble, Journal of Geodesy 2003; 77,

[64] Gross RS. The excitation of the Chandler wobble, Geophysical Research Letter 2000;

[65] Aoyama Y., Naito I. Atmospheric excitation of the Chandler wobble, 1983-1998, Jour‐

[66] Brzezinski A., Nastula J. Oceanic excitation of Chandler wobble, Advances in Space

[67] Gross RS., Fukumori I., Menemenlis D. Atmospheric and oceanic excitation of the Earth's wobbles during 1980-2000, Journal of Geophysical Research 2003; 108, B8,

[68] Seitz F., Schmidt M. Atmospheric and oceanic contributions to Chandler wobble ex‐ citation determined by wavelet filtering, Journal of Geophysical Research 2005; 110,

[69] Salstein DA., Rosen RD. Regional Contributions to the Atmospheric Excitation of Rapid Polar Motions, Journal of Geophysical Research 1989; 94, D7, 9971-9978.

[70] Chao BF., Au AY. Atmospheric Excitation of the Earth's Annual Wobble: 1980-1988,

[71] Nastula J., Kolaczek B. Seasonal Oscillations in Regional and Global Atmospheric Ex‐

[72] Ponte RM., Oceanic excitation of daily to seasonal signals in Earth rotation: results from a constant-density numerical model, Geophysical Journal International 1997;

citation of Polar Motion, Advances in Space Research 2002; 30, 381-386.

physical Journal International 1996; 125, 623-629.

326 Geodetic Sciences - Observations, Modeling and Applications

Geophysical Research Letter 1999; 26, 2085-2088.

454-458.

2370.

B11406.

130, 469-474.

27, 2329-2332.

Research 2002; 30, 195-200.

series of polar motion, Journal of Geodesy 2001; 74, 701-710.

as observed over one century, Surveys in Geophysics 2004; 25, 1-54.

with smaller amplitudes, Journal of Geodesy 2003; 77, 388-401.

nal of Geophysical Research 2001; 106, B5, 8941-8954.

Journal of Geophysical Research 1991; 96, B4, 6577-6582.

