Inertia Tensor for MORVEL Tectonic Plates

Li Chunxiao

(1. Yunnan Observatories, Chinese Academy of Sciences, Kunming 650011, China, Email: lcx@ynao.ac.cn;2. University of Chinese Academy of Scinces, Beijing 100049, China)

?

Inertia Tensor for MORVEL Tectonic Plates

Li Chunxiao1,2

(1. Yunnan Observatories, Chinese Academy of Sciences, Kunming 650011, China, Email: lcx@ynao.ac.cn;2. University of Chinese Academy of Scinces, Beijing 100049, China)

Abstract:The NNR (No-Net-Rotation)-MORVEL (Mid-Ocean Ridge VELocity) 56 is a set of angular velocities describing the motions of 56 plates relative to a No-Net-Rotation reference frame. These plates can be adjusted in terms of non-overlapping polygonal regions, separated by plate boundaries on a unit sphere. During the calculation on the kinematic parameters for these 56 plates in a NNR reference frame using the International Terrestrial Reference Frame (ITRF) velocity field, the geometric parameters of tectonic plates play a significant role in establishing an absolute plate motion model based on space geodesy results. The computational method for these geometric parameters implemented as a FORTRAN90 program is described in this paper, allowing an evaluation of the area and the inertia tensor of a polygonal region on a unit sphere. This program is mainly built on a triangulation algorithm and the adaptive Simpson’s double integral method for spherical polygons, which produces highly reliable results for all 56 modern plates.

Key words:Tectonic plate; Spherical polygon; Inertia tensor; NNR-MORVEL56

1Introduction

Most of Earth’s major features can be understood from the interactions between tectonic plates, which move independently, separating from, colliding with, and sliding against one another. Until the middle 1960s an unifying theory was developed to explain Earth’s dynamics[1]. Several decades after the inception of the theory on plate tectonics, the plate dynamic models constructed using the geological and geophysical data have been dominant, until long time-span geodetic observations were gathered to estimate contemporary plate kinematic parameters[2-5]. As one of the most representative geological plate motion models, NUVEL-1A is one of the mainstream models regarding the plate dynamics and kinematics. With the setting up of the enhanced amount and quality of the geologic or geodetic data during the last few years, the MORVEL refined the precision and accuracy of the geometric and kinematic parameters for 56 plates that are partly taken from an updated digital model of pate boundaries by Bird[6]. Relative to the NUVEL-1A, the MORVEL incorporates more than twice as many plates and covers more of Earth’s surface, and nearly all the NUVEL-1A angular velocities differ significantly from its MORVEL counterparts[7].

To derive an absolute motion model in a NNR reference frame, however, the inertia tensors are always considered as indispensable attribute of all these plates. Despite various established methods for calculating plate inertia tensors corresponding to the NUVEL-1A model presented in many papers[8-9], however it is necessary to recalculate a new set for the NNR-MORVEL56 model[10], given the considerable discrepancy between the NUVEL-1A and the MORVEL.

When a polygon on the unit sphere is employed for the representation of a tectonic plate, a simplified analysis of the plate inertia tensor can be performed through a numerical method, which is carried out over all 56 plates. The method for calculating all 9 components of the inertia tensor is illustrated in this paper and this method requires the precise knowledge of the plate boundaries. The boundary file contains a 2-column sequence of the latitude-longitude plate boundary coordinates that fully enclose the plate in the counterclockwise direction.

In the first section, we introduce some concepts regarding the no-net-rotation conditions and indicate the calculation of the Euler vector in an absolute motion model. The following section describes the detailed mathematical models to estimate the area and the inertia tensor of the spherical polygons. The last section of this paper is dedicated to manifest the results for 56 modern plates and the appendix gives the original FORTRAN90 program for obtaining the aforementioned results.

2Net Lithosphere Rotation

A no-net-rotation model for the lithosphere assumes that the integral ofv×rover the Earth’s surface equals zero[11], i.e.

(1)

where,ris the radial vector of the surface element on a unit sphere, andvcorresponds to the horizontal velocity at that position. The angular velocity of net rotationωnetwas computed as the total angular momentum of all plates divided by the moment of inertia of the entire lithosphere, using the equation[12-13]

(2)

Then it is convenient to convert the Eq.(2) into the following form[8]:

(3)

whereωiis the Euler vector describing the motion of plateirelative to an inertial reference frame, such as ITRF2008, andQiis the inertia tensor of platei, whereigoes from 1 ton. The angular velocities of the plates relative to the NNR reference frame were then found by vector subtraction, namely,

(4)

For those tectonic plates where angular velocities are not available in the geodetic model such as the ITRF2000-PMM, due to a lack of sufficient data, Altamimi[14]tested four cases to perfect the incomplete geodetic model. The fourth case described a method for estimating the missing angular velocity. Here we employ it in this paper by using a simple equation written as:

(5)

3Area and inertia tensor of plates

3.1Evaluation of the plate area

The spherical polygons are defined by great circle arcs connecting points on the sphere, the positions of which are given by latitudes and longitudes. In fact, an algorithm for determining the area of a spherical polygon of arbitrary shape has been presented by Bevis and Cambareri[15], where the kernel idea is to compute the interior angle at each vertex of the spherical polygon. In this paper, however, we employ a somewhat similar method to that of Miller[16], trying to determine the area of spherical polygon by summing the signed areas of component triangles.

For a spherical polygon ofnsides, the spherical excessEis generalized as

(6)

whereαi(i=1…n) are the interior angles of the polygon. Considering a spherical polygon ABCD as shown in Fig.1(a), the north-pole combined with any two adjacent vertices of the polygon can constitute a spherical triangle, such as NAB. The two sides of the triangle are known from the latitudes of their vertices, i.e.an=π/2-latitude(A) andbn=π/2-latitude(B). Taking the previously obtained two sides and the included angle specified by the difference between longitudes ofAandB, the opposite sideabcan be calculated via the formula:

(7)

where the haversine function is defined ashavx=(1-cosx)/2 withxin radians. Having obtained all three sides of the spherical polygon, we can use the formula to obtain its excess:

(8)

To find the area of a spherical polygon, first, one may use the successive vertices in pair to form a spherical triangle. Each spherical triangle employs the north pole as a common vertex to make the calculations convenient. When calculating the areas of the individual triangles, we adopt a convention that the sign of the triangle area(which has a same value as the spherical excess for a polygon on a unit sphere) is identical to the sign of the difference between the longitudes of a pair of adjacent vertices. If the longitude of the first vertex is less than the second one, then the sign of this triangle area is defined as a positive value for a set of points arranged in the anticlockwise way and vice versa. Therefore the area of the spherical polygon isregarded as the absolute value of the sum of the signed spherical excesses for each of the spherical triangles. Taking the facility of the calculation for the upcoming inertia tensor, a provision is crucial that the vertex point traversing the polygon prefersto be enumerated in the counterclockwise direction.

Fig.1Spherical polygon and triangles illustrating calculation method for area and inertia tensor discussed in text.

(a) N triangles constructed from counterclockwise spherical polygon of n sides;

(b) Spherical triangle encompassing Sorth Pole and the sides of polygon traversing the 180th meridian

3.2Estimation of plate inertia tensor

The components of the symmetric inertia tensorQcan be calculated for a regionPusing the following formula:

(9)

wherexμ(μ=1,2,3) are the Cartesian coordinates,δμν(μ,ν=1,2,3) are the elements of the identity matrix, and the integration is carried out over the surface of a plateP. These inertia tensors are based on the hypothesis that the surface density of the plate is unit one, and entirely describe the plate geometry. For instance, the plate areaAis easily calculated by taking the trace ofQ:

(10)

This implies invariance of the trace under coordinate rotations and the sum of the diagonal components is always double the area of the polygon. Generally speaking, non-diagonal components indicate the asymmetry of the polygon with respect to the Cartesian axes, and all of the diagonal components have a positive value, which is useful, together with the Eq.(10), as the verification test for the calculation results.

In this paper we propose a somewhat different method from Schettino[17]for constructing the spherical triangle. In fact, we will see that the integral at the right-hand side of (Eq.(9)) is easily calculated for spherical triangles. The components of the total tensor are therefore given by:

(11)

whereAis the total polygon area, which has been illustrated in the last section. Let a point on the sphere be given in spherical coordinates (θ,λ), whereθis the latitude andλis the longitude, so that its Cartesian coordinates are given by

(12)

then the area element dAat this position is equal to cosθdλdθ. The components of the inertia tensor for a triangleNABare therefore written, in spherical coordinates, as:

(13)

whereλ1,λ2are the longitudes of verticesA,Band the function f is given by

(14)

Next, as a result of the symmetry of the inertia tensor, the 6 independent components of f are expressed in the following way:

(15)

(16)

(17)

f22(θ,λ)=cos3θsin2λ

(18)

(19)

(20)

The upper limit of the inner integral about the latitude is always set to π/2, because the North Pole is considered as the common vertex of each of the spherical triangles. In contrast, with the simple upper limit, the lower limit functionθ(λ) can be obtained from a serious of derivations, whose concrete form is formulated as

(21)

whereC1,C2,C3represent three constants. Once the integrated triangles are determined by one side of the polygon, such asAB, we can write their expression in the following form:

(22)

where (θ1,λ1), (θ2,λ2) are the latitude and the longitude of verticesAandB, respectively.

4Results and analysis

The NNR-NUVEL56 model contains 56 tectonic plates around the earth, and the software OSXStereonet developed by Cardozo and Allmendinger[18]was applied to plot the global plate distribution map, as illustrated in Fig.2. Utilizing the Fortran program, we estimated geometry parameters of all 56 modern plates with the accuracy the inertia tensor better than 10-6. Table 1 lists the area and the inertia tensor of all MORVEL plates, which provide essential material for calculating plate kinematic parameters in the NNR reference frame. The sum of the area of all plates equals 12.566340 steradian, which is slightly less than the surface area of the whole unit sphere, namely 4π(12.566371) with the relative errorη=0.00023%. Our results indicate that the relative errors for the six components of the total tensor are 0.00017%, 0.00029%, 0.00026%, 0.0010%, 0.00016%, and 0.00010% respectively, It is shown that the inertia tensor of the entire spherical surface is 8π/3E, whereEtakes the identity matrix. The discrepancy probably arises from either the imperfection of plates over the entire sphere or the unavoidable rounding errors in floating point arithmetic.

5Conclusion

The method for computing the areas and inertial tensors of tectonic plates has been presented. This method is based upon the triangulation algorithm and the adaptive Simpson’s double integral procedure, which can be applied to the spherical polygons representing such tectonic plates. Results for the NNR-MORVEL56 tectonic plates show that highly reliable data can be produced, as long as starting from the precise definition of the plate boundaries. In addition, a FORTRAN90 program has been attached to the end of thispaper, which is expected to be valuable to the future studies of the kinematics and dynamics associated to the motions of tectonic plates.

Fig.2Plate boundaries and geometries employed for MORVEL

Table 1 continued from previous page.

aPlate name abbreviations are as follows: am, Amur; an, Antarctic; AP, Altiplano; ar, Arabia; AS, Aegean Sea; au, Australia; BH, Birds Head; BR, Balmoral Reef; BS, Banda Sea; BU, Burma; ca, Caribbean; CL, Caroline; cp, Cocos; cp, Capricorn; CR, Caroline; EA, Easter; eu, Eurasia; FT, Futuna; GP, Galapagos; in, India; jf, Juan de Fuca; JZ, Juan Fernandez; KE, Kermadec; lw, Lwandle; MA, Mariana; MN, Manus; MO, Maoke; mq, Macquarie; MS, Molucca Sea; na, North America; NB, North Bismarck; ND, North Andes; NH, New Hebrides; NI, Niuafoou; nb, Nubia; nz, Nazca; OK, Okhotsk; ON, Okinawa; pa, Pacific; PM, Panama; ps, Philippine Sea; ri, Rivera; sa, South America; SB, South Bismarck; sc, Scotia; SL, Shetland; sm, Somalia; sr, Sur; SS, Solomon Sea; su, Sundaland; sw, Sandwich; TI, Timor; TO, Tonga; WL, Woodlark; yz, Yangtze. Plate abbreviations given in lower case are for plates included in the MORVEL. Plate abbreviations given in the upper case are for plates from Bird (2003); b. Plate areas are in steradians for a unit sphere, and the sum of which totals 4π.

References：

[1]Hamblin W K, Christiansen E H. Earth′s dynamic systems[M]. 10th ed. State of New Jersey: Prentice Hall, 2003.

[2]Larson K M, Freymueller J T, Philipsen S. Global plate velocities from the global positioning system[J]. Journal of Geophysical Research Solid Earth, 1997, 102(B5): 9961-9981.

[3]Sella G F, Dixon T H, Mao A. Revel: a model for recent plate velocities from space geodesy[J]. Journal of Geophysical Research Solid Earth, 2002, 107(B4): ETG 11-1- ETG 11-30.

[4]Corné K, Holt W E, John H A. An integrated global model of present-day plate motions and plate boundary deformation[J]. Geophysical Journal International, 2003, 154(1):8-34.

[5]Altamimi Z, Métivier L, Collilieux X. ITRF 2008 plate motion model[J]. Journal of Geophysical Research Solid Earth, 2012, 117(B7): 47-56.

[6]Peter B. An updated digital model of plate boundaries[J]. Geochemistry Geophysics Geosystems, 2003, 4(3):101-112.

[7]Demets C, Gordon R G, Argus D F. Geologically current plate motions[J]. Geophysical Journal International, 2010, 181(1): 1-80.

[8]Argus D F, Gordon R G. No-net-rotation model of current plate velocities incorporating plate motion model nuvel-1[J]. Geophysical Research Letters, 1991, 18(11): 2039-2042.

[9]Corné K, Holt W E. A no-net-rotation model of present-day surface motions[J]. Geophysical Research Letters, 2001, 28(23): 4407-4410.

[10]Argus D F, Gordon R G, Demets C. Geologically current motion of 56 plates relative to the no-net-rotation reference frame[J]. Geochemistry Geophysics Geosystems, 2011, 12(11): 75-87.

[11]Klemann V, Martinec Z, Ivins E R. Glacial isostasy and plate motion[J]. Journal of Geodynamics, 2008, 46(3-5): 95-103.

[12]Solomon S C, Sleep N H. Some simple physical models for absolute plate motions[J]. Journal of Geophysical Research, 1974, 79(17): 2557-2567.

[13]Torsvik T H, Steinbergerd B, Gurnise M, et al. Plate tectonics and net lithosphere rotation over the past 150 my[J]. Earth & Planetary Science Letters, 2010, 291(1): 106-112.

[14]Zuheir A, Patrick S, Claude B. The impact of a no-net-rotation condition on ITRF2000[J]. Geophysical Research Letters, 2003, 30(2): 36-1-36-4.

[15]Bevis M, Cambareri G. Computing the area of a spherical polygon of arbitrary shape[J]. Mathematical Geology, 1987, 19(4): 335-346.

[16]Heckbert P S. Graphics Gems IV[M]. Boston: Academic Press Professional, 1994: 132-137.

[17]Schettino A. Computational methods for calculating geometric parameters of tectonic plates[J]. Computers & Geosciences, 1999, 25(8): 897-907.

[18]Cardozo N, Allmendinger R W. Spherical projections with OSXStereonet[J]. Computers & Geosciences, 2013, 51: 193-205.

CN 53-1189/PISSN 1672-7673

MORVEL構(gòu)造板塊的轉(zhuǎn)動張量

李春曉1,2

(1. 中國科學(xué)院云南天文臺，云南 昆明650011；2. 中國科學(xué)院大學(xué)，北京100049)

關(guān)鍵詞:構(gòu)造板塊；球面多邊形；轉(zhuǎn)動慣量張量；NNR-MORVEL56

摘要:NNR-MORVEL56板塊運動模型描述了全球56個構(gòu)造板塊在無整體旋轉(zhuǎn)參考架下的角速度運動參數(shù)。這些板塊可以近似描述為單位球上的無重疊球面多邊形區(qū)域。用ITRF速度場計算這56個板塊相對于無整體旋轉(zhuǎn)參考架下的絕對運動時，板塊的幾何參數(shù)起著至關(guān)重要的作用。詳細給出了計算板塊幾何參數(shù)的方法并且編寫了FORTRAN90程序以供參考，使得計算單位球上板塊的面積和轉(zhuǎn)動慣性張量得以實現(xiàn)。文中的計算方法和程序主要采用球面三角算法和自適應(yīng)辛普森雙積分算法，并對全球56個板塊的幾何參數(shù)進行了計算，得到了較為可靠的計算結(jié)果。

基金項目：中國科學(xué)院精密導(dǎo)航定位與定時技術(shù)重點實驗室青年基金 (2014PNTT09) 資助.

收稿日期：2015-03-26；修定日期：2015-04-22

作者簡介：李春曉，男，碩士. 研究方向：地球動力學(xué). Email: lcx@ynao.ac.cn

中圖分類號：P183.2

文獻標識碼：A

文章編號:1672-7673(2016)01-0058-12