The exact implementation of a spherical harmonic model for gravimetric quantities

1 Introduction

Gravimetric quantities are defined as physical quantities related to Earth’s gravitational field that can be measured directly. To illustrate their corresponding evaluation using a spherical harmonic expansion (SHE) of the field, consider the geodetically relevant quantities, the magnitude of gravity, the components of the deflection of the vertical (DOV), and the normal height. Gravity magnitude is measured by a gravimeter; deflection components are determined by measuring astronomic coordinates at a point; and the normal height follows immediately from measurements of geopotential differences (formulas are provided in Section 2). In all cases, the global three-dimensional coordinates of the measurement point must be known, which today is easily achieved by tracking satellites of a Global Navigation Satellite System, such as GPS (Global Positioning System).

Several approximations are usually introduced in order to make the development of an SHE from gravimetric data reasonably tractable. Although the development of the best such models incorporates appropriate corrections to these approximations, one often reverts to these approximations when using the models to compute the measured quantities. Usually, such computations are based on their relationships to the disturbing potential (e.g., Ivanov et al. 2018) and thus are corrupted by spherical and/or linear approximations. Generally, the spherical approximation error is of the order of Earth’s flattening times the disturbing (or anomalous) quantity being computed, while the linear approximation error is usually less and of the order of the square of the disturbing (or anomalous) quantity. The present discussion concerns evaluating the actual measured quantity using an SHE without these usual approximations. This permits direct comparison of the model with measured data, thus enabling a more direct evaluation of the accuracy of the model, assuming that the formulation of the model is exact. That is, the only assumed model error is due to errors in the estimation of the model parameters, which includes the fact that only a finite number of parameters, from the theoretically infinite set, can be estimated. The effect of these model errors on the computed quantities is outside the scope of this report, but it is worth knowing that these are the only errors that could affect the computed quantities and that spherical and linear approximations are absent.

The basic SHE for Earth’s gravitational field as a finite series of spherical harmonic functions for the gravitational potential, V , is given in spherical polar coordinates, θ , λ , r , by

where θ is the co-latitude of the evaluation point, λ is its longitude, and r is its radius; GM is the product of Newton’s constant of gravitation and Earth’s total mass (including the atmosphere); R is a specified constant radius; and the constants, C n , m , are real coefficients associated with the spherical harmonic functions, Y ¯ n , m ( θ , λ ) , defined by

where P ¯ n , m is a fully normalized associated Legendre function of the first kind. The normalization is such that the spherical harmonic functions are orthonormal on the unit sphere. The gravity potential, W , is the gravitational potential plus the centrifugal potential due to Earth’s rate of rotation, ω E ,

The disturbing potential is then defined with respect to a normal gravity potential, U , by T = W − U . The normal gravity potential is generated by a co-rotating normal spheroid that is also an equipotential surface of the normal gravity field ( U = U 0 ). Details of these definitions may be found in the textbook by Hofmann-Wellenhof and Moritz (2005).

2 Gravimetric quantities

The geopotential number is defined by

where W 0 is the gravity potential that (regionally or globally) defines the geoid (a vertical datum). The geopotential number, C P ( W 0 ) , at a point, P , can be measured using its equivalent definition as a line integral in a conservative field,

where g = − ∂ W / ∂ ν is the component of gravity in the direction of the gradient of W .

Having measured the geopotential number at a point, the dynamic height follows immediately,

where γ 0 is a constant value of gravity that transforms the units of the geopotential number from (m²/s²) to (m). A more geometrical height, the normal height, associated with P and relative to the vertical datum defined by W 0 , is given by

where γ ¯ P ′ is the average value of normal gravity along the normal plumb line between the telluroid point, P ′ , and the normal spheroid. A formula that requires only the normal gravity, γ P 0 , on the normal spheroid is given by a truncated series (Hofmann-Wellenhof and Moritz 2005, p. 168).

where m = ω E 2 a 2 b / GM , and the parameters, a , b , f , GM , belong to the normal spheroid. The last term in the outside parentheses is less than ( H P ⁎ / a ) 2 < 2 × 10 − 6 for all topographic elevations; and thus, the relative accuracy of equation (8) is of the order of ( H P ⁎ / a ) 3 ≃ 3 × 10 − 9 , which is two orders better than millimeter accuracy in the normal height. This is better than the linear approximation inherent in the relationship of the normal height to the disturbing potential.

where h P is the geodetic height above a defined geodetic ellipsoid, and the second term on the right side is also known as the height anomaly.

The gravity vector in the spherical coordinate system, θ , λ , r , is the gradient of the gravity potential, W ,

where the gradient operator for spherical coordinates is ∇ ( θ λ r ) = 1 r ∂ ∂ θ 1 r sin θ ∂ ∂ λ ∂ ∂ r T ; hence, from equation (3),

where the components of the centrifugal acceleration are in the coordinate directions of θ , λ , r , respectively. The partial derivatives of equation (1), needed in equation (11), are readily expressed with respect to the longitude and the radius. For the co-latitude, a formula for the derivative of the associated Legendre function is given for 0 ≤ m ≤ n and n ≥ 1 by Holmes and Featherstone (2002).

The magnitude of gravity at a point (the quantity that is measured) is simply

The gravity anomaly and gravity disturbance then follow immediately with an appropriate subtraction by the normal gravity magnitude (Hofmann-Wellenhof and Moritz 2005). For the free-air gravity anomaly, the standard formula used for evaluating the corresponding SHE is

which includes both the spherical and the linear approximations.

The deflection of the vertical (DOV) has various definitions (Jekeli 1999); consider the Helmert definition, which is the angle between the direction of gravity at a point and the ellipsoid normal through that point for a given geodetic ellipsoid (Figure 1). The north and east components, ξ and η , can be determined with measurements of astronomic latitude and longitude, Φ and Λ , respectively (Pick et al. 1973),

where ϕ P , λ P are the geodetic coordinates with respect to the defined geodetic ellipsoid. The third-order terms may be ignored; even the second-order term in equation (15) for maximum deflection angles of 1 ′ is well below the observation accuracy of about 0 . 1 ″ (Hirt et al. 2004) and is usually neglected. In the local, south-east-up, geodetic coordinate system defined with the up direction along the ellipsoid normal (Figure 1), the gravity vector has components along mutually orthogonal unit vectors, ψ ˆ P , λ ˆ P , h ˆ P , where ψ = 90 ° − ϕ is the geodetic co-latitude,

where δ s ψ and δ s λ are differential elements in the south and east directions, respectively, and δ h is a differential element in the up direction. Transforming from the spherical coordinate directions to the south-east-up directions involves a rotation about the east axis by the angle, ψ P − θ P ; therefore,

Figure 1

The components of the gradient of the gravity potential, drawn such that the deflections of the vertical, ξ , η , are positive according to equations (15) and (16).

From Figure 1, the components of the Helmert DOV thus satisfy

These are exact formulas, but the differences with respect to the measurable quantities equations (15) and (16) are negligible.

These formulas differ from the standard formulas used to evaluate the SHE for the DOV components,

which contain spherical and linear approximations (Jekeli 1999).

Given the geodetic coordinates, ϕ P , λ P , h P , of a point, P , its spherical coordinates are

where e 2 = 2 f − f 2 is the square of the first eccentricity of the geodetic ellipsoid and N P = a / 1 − e 2 sin 2 ϕ P , where a , f are the semi-major axis and flattening of the ellipsoid. Thus, the gravimetric quantities, geopotential number, gravity magnitude, and components of the DOV can be readily computed from the SHE equation (1) using, respectively, equation (4), equations (11) and (13), and equations (19) and (20).

3 Numerical precision of evaluating the SHE

While the SHE-derived formulas for the gravimetric quantities are exact, there remains the question of the numerical precision of computing the SHE. This question arises particularly as the model parameter, n max , increases and the recursion formulas for the associated Legendre function and its derivative tend to deteriorate in precision. Using Fukushima’s extended-exponent algorithm (Fukushima 2012), and a comparison of double and quadruple precision in FORTRAN, the number of significant digits in the double precision evaluation of the model, EGM2008 ( n max = 2 , 190 ) (Pavlis et al. 2012), over the full range of latitudes is established. The reliability of FORTRAN quadruple precision as a standard for comparison is verified by coding the same evaluation in the programming language, Python, using 32-digit arithmetic. Figure 2 shows the absolute relative differences between the double and quadruple precision values of the gravity magnitude and the DOV components for geodetic latitudes from equator to pole. Better than 12 significant digits of accuracy are obtained for the gravity and east deflection components, while generally 10 or more significant digits are obtained for the north component (likely because derivatives of the associated Legendre functions are involved exclusively). In all cases, the numerical accuracy of evaluating gravimetric quantities using EGM2008 as described above is many orders of magnitude better than the corresponding measurement accuracy.

Figure 2

Absolute relative numerical accuracy in gravimetric quantities evaluated by EGM2008 ( n max = 2 , 190 ) .

Finally, it should be noted that the provided FORTRAN code developed for these computations is designed for a relatively small number of data points, commensurate with the amount of data gathered with typical ground surveys. It is not designed for efficiency in generating very large grids of data. Of course, efficiency is a fluid concept in view of ever-increasing computational capabilities. The double precision values associated with Figure 2 were achieved at an average rate of 0.25 s per point running in a Microsoft Windows environment using the Microsoft Fortran 77 compiler on a laptop with an Intel 8th generation i5 core processor running at 1.6 GHz and with 16 Gbytes of random access memory (actual central processing unit (CPU) time is vastly less).

4 Summary

A spherical harmonic expansion (SHE) of Earth’s gravitational potential is used to formulate exact expressions for gravimetric quantities, such as the geopotential number (which leads immediately to the dynamic height and also the normal height), the magnitude of gravity, and the components of the deflection of the vertical. This contrasts with the usual corresponding expressions derived from the SHE that are corrupted by spherical and linear approximations. Other quantities of interest, such as the gravity anomaly and geoid undulation, can be derived from these expressions, but only with additional assumptions. For example, the geoid undulation requires an assumption on the mass-density of the crust; and the free-air gravity anomaly typically requires a model for the vertical gradient of gravity. Other quantities, such as the Molodensky gravity anomaly, require only the three-dimensional definition of a normal field. Further derivatives of the potential, such as the elements of the gradient tensor, follow similarly as derived here. A FORTRAN program and sample output based on the full EGM2008 model are given as supplemental material that illustrates the computations for the basic quantities discussed in this report. It is also demonstrated that such computations (using the Fukushima extended-exponent algorithm) yield a minimum of 9 significant digits, and typically 10 to 13, when using the full EGM2008 expansion.