Theoretical gravity
In geodesy and geophysics, theoretical gravity (also known as normal gravity) is a mathematical model approximating the gravity field of the Earth, computed for points on or above a reference ellipsoid of revolution.[1] This ideal gravity includes the gravitational attraction of a uniform ellipsoid plus the centrifugal effect due to Earth's rotation, providing a standard against which observed gravity measurements are compared to determine gravity anomalies.[2] The concept originated from Isaac Newton's theory of universal gravitation in the 17th century, with significant developments in the 18th and 19th centuries through efforts to model Earth's oblate shape, such as Clairaut's theorem on the figure of the Earth.[3] International standardization began in the 20th century with formulas like the International Gravity Formula of 1930, evolving to modern systems such as the Geodetic Reference System 1980 (GRS80) and World Geodetic System 1984 (WGS84).[4] These models are essential for applications in surveying, navigation, and understanding Earth's mass distribution.Introduction
Definition and Scope
Theoretical gravity represents the idealized model of the gravitational acceleration on Earth's surface, computed as the normal gravity field arising from a reference ellipsoid that approximates the planet as a rotating oblate spheroid, excluding local geological anomalies or irregularities.[1] This model provides a standardized value for gravitational acceleration, with the internationally adopted standard of 9.80665 m/s² defined at sea level on the reference ellipsoid.[5] In essence, theoretical gravity serves as the baseline expectation for g, derived purely from the Earth's overall mass distribution, rotational dynamics, and ellipsoidal shape, without perturbations from subsurface density variations or topography.[2] A key distinction lies in how theoretical gravity differs from actual measurements: while measured gravity incorporates local anomalies caused by mass excesses or deficits (such as those from mountains, ocean trenches, or mineral deposits), theoretical gravity deliberately omits these to represent a smooth, idealized field.[6] Similarly, free-air gravity refers to observed values corrected only for elevation above the reference surface but still retaining anomaly effects, whereas theoretical gravity remains tied to the ellipsoid's predicted norm.[7] This separation enables the isolation and study of gravitational deviations for geophysical analysis. The scope of theoretical gravity extends beyond Earth, adapting the reference model to other celestial bodies through advanced representations. For planets like Mars or Jupiter, spherical harmonic expansions model the gravity field by decomposing it into Legendre polynomials that capture global mass distributions and oblateness due to rotation.[8] For irregular bodies such as the Moon or asteroids, mascon (mass concentration) models parameterize localized gravity variations as discrete mass elements, offering efficient approximations where uniform ellipsoids fail.[9] These extensions maintain the core principle of an anomaly-free baseline, tailored to each body's shape and dynamics. Theoretical gravity thus underpins geodesy by establishing reference frameworks for precise positioning and Earth orientation studies.Historical Development
The development of theoretical gravity models began in the early 20th century with efforts to model Earth's gravitational field using reference ellipsoids derived from ground-based measurements. Friedrich Robert Helmert's 1906 ellipsoid, with a semi-major axis of 6,378,200 m and flattening of 1/298.3, was based on data from over 1,600 pendulum stations worldwide, correcting for altitude and uneven distribution of observations to estimate Earth's ellipticity.[10] This model served as a key precursor in European geodesy, providing an empirical foundation for normal gravity calculations despite limitations from sparse global coverage. Similarly, John Fillmore Hayford's 1924 ellipsoid, featuring a semi-major axis of 6,378,388 m and flattening of 1/297, utilized deflections of the vertical from North American triangulation networks, achieving an ellipticity estimate of 1/297 ± 0.5. Adopted by the International Association of Geodesy (IAG) in 1924 as the International Reference Ellipsoid, it marked a milestone in standardizing global models by converging results from arc measurements, pendulum gravity, and astronomical perturbations.[10][11] Pre-1930 approximations, including those tied to Helmert and Hayford ellipsoids, relied on incomplete datasets and were gradually replaced as improved gravimetric and deflection measurements revealed inconsistencies in Earth's oblateness and equatorial radius estimates. By the late 1920s, advancements in instrumentation and data analysis, such as better control of topographic and density perturbations, necessitated a unified gravity formula. In 1930, the IAG adopted the International Gravity Formula (IGF) at its Stockholm General Assembly, proposed by Gino Cassinis and based on the 1924 Hayford ellipsoid (also known as the Cassinis ellipsoid in this context). This formula provided a standardized expression for normal gravity on the ellipsoid, incorporating the Potsdam gravity datum and Clairaut's spheroid theory, to facilitate consistent anomaly computations worldwide.[12][13] Key revisions to the IGF occurred in the mid- to late 20th century, driven by satellite geodesy and refined Earth models. The 1967 revision, adopted by the IAG at Lucerne, tied the IGF to the Geodetic Reference System 1967 (GRS67), with parameters including a semi-major axis of 6,378,160 m and flattening of 1/298.247, reflecting enhanced gravity data from global networks and early satellite observations that improved accuracy over the 1930 version by accounting for better-determined dynamical form factors.[14] Further refinement came in 1980 with the GRS80-based IGF, featuring a semi-major axis of 6,378,137 m and flattening of 1/298.257, which incorporated Doppler satellite data and higher-precision gravimetry to minimize discrepancies in the normal gravity field.[13] Following 1980, the World Geodetic System 1984 (WGS84), developed by the U.S. Department of Defense and released in 1984, refined GRS80 parameters slightly—adjusting the inverse flattening to 1/298.257223563—for compatibility with GPS navigation, while adopting its gravity formula for global applications. No new standardized IGFs have been adopted since 1980, as WGS84 and GRS80 continue to underpin modern geodesy, supported by ongoing satellite missions like GRACE that validate their adequacy without necessitating wholesale replacement.[15][16]Physical Principles
Gravitational Components
In theoretical gravity models, the pure gravitational attraction arises from the Earth's mass distribution, approximated initially by Newton's law of universal gravitation. For a spherical Earth of mass M and radius r, the gravitational acceleration g_{\text{grav}} at the surface is given by g_{\text{grav}} = \frac{GM}{r^2}, where G is the gravitational constant, yielding a value of approximately 9.80665 m/s² for the standard Earth parameters GM = 3.986004418 \times 10^{14} m³/s² and mean radius r \approx 6371 km.[17] This expression assumes uniform density and sphericity, but the Earth's oblateness requires adjustments to account for the non-uniform mass distribution. To incorporate the Earth's oblate shape, the gravitational potential is expanded in spherical harmonics, with the dominant correction from the second-degree zonal harmonic coefficient J_2, representing the equatorial bulge. The adjusted gravitational acceleration becomes g_{\text{grav}} \approx \frac{GM}{r^2} \left[ 1 - 3 J_2 \left( \frac{a}{r} \right)^2 P_2 (\sin \phi) + \ higher\ order\ terms \right], where a is the equatorial radius, \phi is the geocentric latitude, P_2 (x) = \frac{1}{2} (3 x^2 - 1) is the Legendre polynomial of degree 2, and J_2 \approx 1.08263 \times 10^{-3} quantifies the quadrupole moment due to rotational flattening.[17] This J_2 term increases gravity at the poles and decreases it at the equator, reflecting the closer proximity to the denser core at higher latitudes. The oblate spheroid configuration couples the Earth's shape to its gravity field through hydrostatic equilibrium, as described by Clairaut's theorem, which relates the flattening f of the ellipsoid to the gravitational acceleration's latitudinal variation. Specifically, Clairaut's theorem states that the ratio of the centrifugal to gravitational acceleration at the equator drives the polar excess in gravity, approximately \frac{\Delta g}{g} \approx \frac{5}{2} m - f, where m = \frac{\omega^2 a}{g} is the centrifugal ratio and \omega is Earth's angular velocity; this predicts a 0.5% stronger gravity at the poles due to both geometric flattening and mass redistribution.[18] The theorem, derived under the assumption of a self-gravitating, rotating fluid body, ensures the equipotential surface aligns with the ellipsoid, providing a foundational link between form and force in theoretical models. The normal gravity field in theoretical models is computed as the radial derivative of the gravitational potential U evaluated at the reference ellipsoid surface, yielding \gamma = -\frac{\partial U}{\partial r}, where U is the sum of the central Newtonian term and the oblateness corrections up to J_2.[19] This potential U is normalized such that U = \frac{GM}{r} \left(1 - J_2 \left(\frac{a}{r}\right)^2 [P_2](/page/polynomial)(\sin \phi)\right), with P_2 the Legendre polynomial, ensuring the field represents the idealized, rotationally symmetric attraction without local anomalies.[20] The resulting \gamma(\phi) varies smoothly from equatorial to polar values, establishing the baseline for geodetic computations.Rotational Effects
The rotation of Earth introduces a centrifugal acceleration that acts outward perpendicular to the axis of rotation, modifying the effective gravitational field experienced at the surface. This acceleration arises in the non-inertial rotating frame of reference and has magnitude a_c = \omega^2 \rho, where \omega is the angular velocity of Earth and \boldsymbol{\rho} is the perpendicular distance vector from the rotation axis to the point of interest. For a point at latitude \phi, \rho \approx R \cos \phi, where R is the Earth's radius, so a_c = \omega^2 R \cos \phi. The component along the local radial direction (opposing gravity) is a_c \cos \phi = \omega^2 R \cos^2 \phi.[21] Earth's angular velocity \omega is precisely $7.2921151467 \times 10^{-5} rad/s, corresponding to one rotation every sidereal day.[22] The centrifugal acceleration reaches its maximum at the equator (\phi = 0), where it equals approximately $0.034 m/s², directed outward along the radial direction and thus directly opposing the local gravitational attraction.[23] This results in an effective gravity g_{\text{eff}} that is the vector sum of the true gravitational acceleration g_{\text{grav}} (due to Earth's mass distribution) and the centrifugal term, approximated as g_{\text{eff}} \approx g_{\text{grav}} - \omega^2 R \cos^2 \phi for the vertical component, with the centrifugal effect reducing g_{\text{eff}} by up to about 0.3% globally (neglecting the small horizontal component).[21] At higher latitudes, the effect diminishes as \cos^2 \phi, vanishing entirely at the poles. In theoretical models, the gravitational attraction provides the baseline inward pull, but rotation introduces this dynamic centrifugal reduction that must be subtracted to obtain the observed effective gravity.[21] Additionally, Earth's rotation drives the formation of an equatorial bulge through hydrostatic equilibrium, where centrifugal forces cause the planet to assume an oblate shape, indirectly influencing the gravitational field by altering the mass distribution and thus g_{\text{grav}} itself.[24] This oblateness amplifies the latitudinal variation in effective gravity beyond the direct centrifugal component alone.Reference Ellipsoid Models
Reference ellipsoid models provide a geometric approximation of Earth's shape as an oblate spheroid, which is essential for computing theoretical gravity by defining a smooth, equipotential surface that closely matches the geoid. An oblate spheroid is characterized by its equatorial semi-major axis a and flattening f, where f = (a - b)/a and b is the polar semi-minor axis. For the WGS84 reference ellipsoid, a = 6378137 m and f = 1/298.257223563.[25] These parameters ensure the ellipsoid's surface is an equipotential, allowing the normal gravitational potential to be calculated analytically.[26] Several historical and modern reference ellipsoids have been developed to refine this approximation, each tailored to available geodetic data and intended for defining the normal gravity potential in theoretical models. The Hayford ellipsoid of 1909, with a = 6378388 m and f = 1/297, was derived from deflection of the vertical measurements and served as a foundational model for early 20th-century gravity computations.[27] Building directly on Hayford's work, the International Ellipsoid of 1924, adopted by the International Union of Geodesy and Geophysics (IUGG) in Madrid, retained the same parameters (a = 6378388 m, f = 1/297) and became the standard for global gravity reference until the mid-20th century, enabling consistent normal potential calculations across international surveys.[28][29] Later refinements incorporated satellite and more precise ground measurements. The Geodetic Reference System 1967 (GRS67), with a = 6378160 m and f = 1/298.247167427, was established by the IUGG to better align with emerging global data, providing an improved basis for normal gravity potentials in geophysical applications.[25] The Geodetic Reference System 1980 (GRS80), featuring a = 6378137 m and f = 1/298.257222101, further enhanced accuracy by integrating Doppler satellite observations and astro-geodetic data, serving as a reference for deriving closed-form normal potentials that approximate Earth's irregular gravity field.[25] These ellipsoids define the normal potential U such that its level surface coincides with the ellipsoid, allowing theoretical gravity to be computed as the gradient of U on that surface.[30] The following table summarizes the key parameters of these reference ellipsoids:| Ellipsoid | Semi-major Axis a (m) | Flattening f |
|---|---|---|
| Hayford 1909 | 6378388 | 1/297 |
| International 1924 | 6378388 | 1/297 |
| GRS67 | 6378160 | 1/298.247167427 |
| GRS80 | 6378137 | 1/298.257222101 |
Core Formulas
Basic Gravity Expression
The foundational mathematical expression for theoretical gravity provides a simplified estimate of the effective gravitational acceleration g(\phi) at latitude \phi on Earth's surface, prior to more refined models. This basic form is given by g(\phi) \approx g_e \left(1 + \beta \sin^2 \phi \right), where g_e represents the equatorial gravity (approximately 9.780 m/s²), and \beta is the latitude coefficient, empirically and theoretically determined to be about 0.0053 for Earth.[31][32] This approximation captures the primary latitudinal variation, with gravity increasing toward the poles due to Earth's oblate shape and rotation.[33] The derivation of this expression stems from classical potential theory, as developed in Clairaut's theorem (1743), which models Earth as a self-gravitating, rotating fluid body in hydrostatic equilibrium. The total potential \Psi is the sum of the gravitational potential V (from mass distribution) and the centrifugal potential \Phi_c = -\frac{1}{2} \omega^2 s^2, where \omega is Earth's angular velocity and s = r \cos \phi is the distance from the rotation axis.[31] For a basic outline without higher-order series expansions, the gravitational potential is approximated for a nearly spherical body using low-degree spherical harmonics (primarily the J₂ term for oblateness), assuming the surface is an equipotential. The effective gravity g(\phi) is then the magnitude of the gradient of \Psi, projected normal to the surface, yielding the linear dependence on \sin^2 \phi after combining terms and evaluating at the reference radius.[32] This approach integrates the gravitational attraction, which strengthens at higher latitudes due to closer proximity to the center along the polar axis, with the outward centrifugal effect that is maximal at the equator.[33] This simplified expression presupposes a uniform density approximation for the reference Earth model, treating the planet as a homogeneous ellipsoid to facilitate the initial potential calculation and ellipticity estimation via Clairaut's relation between dynamical flattening and the J₂ gravitational coefficient.[31] Such assumptions enable a closed-form estimate suitable for preliminary geophysical analyses, though real Earth deviations (e.g., density contrasts) necessitate corrections in advanced models.[32] The centrifugal reduction at the equator, contributing about 0.3% to the total variation, underscores the rotational influence embedded in \beta.[33]Somigliana Equation
The Somigliana equation provides a closed-form expression for the magnitude of normal gravity on the surface of a reference ellipsoid, accounting for both the ellipsoidal shape and Earth's rotation. It refines simpler latitude-dependent models by incorporating the exact geometry of the equipotential surface. This formula is fundamental in geodesy for defining theoretical gravity in systems like GRS80.[30] The equation for normal gravity g(\phi) at geodetic latitude \phi is given by g(\phi) = g_e \frac{1 + k \sin^2 \phi}{\sqrt{1 - e^2 \sin^2 \phi}}, where g_e is the normal gravity at the equator, g_p is the normal gravity at the poles, k = \frac{b g_p - a g_e}{a g_e}, e^2 = 2f - f^2 is the squared eccentricity, a is the semi-major axis, b is the semi-minor axis, and f is the flattening. For the GRS80 ellipsoid, the parameters are g_e = 9.780327 m/s², g_p = 9.832185 m/s², a = 6378137 m, f = 1/298.257222101, yielding b \approx 6356752.3142 m, e^2 \approx 0.00669438002290, and k \approx 0.001931852653.[30][34] This formula derives from the normal gravitational potential U of a rotating, homogeneous ellipsoid, where normal gravity is the norm of the gradient |\nabla U| evaluated on the ellipsoidal surface. The derivation integrates Clairaut's theorem, which relates the flattening f to the centrifugal potential and mass distribution, with the boundary condition that the ellipsoid is an equipotential surface. Parameter computation involves solving for g_e and g_p from the total potential constants, ensuring consistency with observed gravity values.[35] The Somigliana equation is valid precisely on the ellipsoid surface (height h = 0) and achieves high accuracy, with differences from series expansions or numerical integrations below $10^{-6} m/s² (approximately 0.1 mGal), sufficient for most geodetic applications.[30]Approximation Methods
Series Expansions Overview
Series expansions provide a practical method for approximating theoretical gravity by expanding the closed-form Somigliana equation into a power series in terms of latitude-dependent terms, facilitating easier numerical evaluation.[36] This approach treats the Somigliana formula as the base for deriving normal gravity on the reference ellipsoid.[30] The Taylor series expansion is typically performed around the equator (latitude φ = 0), yielding the form g(\phi) = g_e \left[1 + \alpha \sin^2 \phi + \beta \sin^4 \phi + \cdots \right], where g_e is the equatorial gravity, and the coefficients α, β, etc., are derived from the ellipsoid's geometric parameters such as flattening and rotational effects.[36] Common truncations occur at second order (up to \sin^2 \phi) for basic applications or fourth order (up to \sin^4 \phi) for improved accuracy, balancing computational simplicity with precision.[30] These expansions offer advantages in computational efficiency, particularly for manual calculations or early electronic devices, as they replace complex closed-form evaluations with straightforward polynomial arithmetic.[36] The series converges for latitudes |φ| < 90°, enabling reliable approximations across most of the ellipsoid surface without requiring iterative methods.[30] However, truncating the series introduces limitations, with reduced precision near the poles where higher-order terms become significant, potentially leading to errors larger than those from the full closed-form expression.[36] Historically, such expansions were prevalent before widespread computer availability, serving as essential tools for geodetic computations in the mid-20th century.[30]International Gravity Formulas
The International Gravity Formulas (IGFs) represent a series of empirical approximations derived from series expansions of the Somigliana equation, providing standardized expressions for normal gravity as a function of geodetic latitude φ on specific reference ellipsoids. These formulas were developed by the International Association of Geodesy (IAG) to facilitate consistent gravity anomaly computations in geodesy and geophysics, evolving through refinements in Earth models and measurement data. The initial formula, adopted in 1930, marked the first global standardization effort. The 1930 IGF, based on the International Reference Ellipsoid (also known as the Cassinis or Hayford ellipsoid of 1909), is given by: g(\phi) = 9.78049 \left(1 + 0.0052884 \sin^2 \phi - 0.0000059 \sin^2 2\phi \right) \, \mathrm{m/s^2} This formula incorporates the Potsdam gravity datum and Clairaut's spheroid model, achieving an accuracy of less than 0.1 mGal (where 1 mGal = 10^{-5} m/s²) for latitude-dependent variations. It was designed primarily for unifying gravity measurements from disparate national surveys, though its reliance on pre-satellite data limited precision in equatorial and polar regions.[37] Subsequent updates addressed these limitations with improved ellipsoidal parameters and observational data. The 1967 IGF, tied to the Geodetic Reference System 1967 (GRS67), refines the expression as: g(\phi) = 9.780318 \left(1 + 0.0053024 \sin^2 \phi - 0.0000058 \sin^2 2\phi \right) \, \mathrm{m/s^2} This version, approved by the IAG, enhances accuracy to a maximum error of 0.004 mGal in its precise form, or 0.1 mGal in the conventional approximation, by integrating early satellite observations and adjusting for better rotational and centrifugal effects. The shift from the 1930 formula introduces systematic differences of up to 17 mGal, primarily due to updated equatorial gravity values and flattening parameters in GRS67.[37][30] The 1980 IGF, associated with the Geodetic Reference System 1980 (GRS80), further optimizes the series for modern applications: g(\phi) = 9.780327 \left(1 + 0.0053024 \sin^2 \phi - 0.0000058 \sin^2 2\phi \right) \, \mathrm{m/s^2} Retaining the core structure of the 1967 version but with a slightly adjusted equatorial constant, it achieves a relative accuracy of 10^{-4} mGal (0.1 μGal), representing an order-of-magnitude improvement over prior formulas through incorporation of high-precision satellite altimetry and gravimetry data. Comparisons reveal differences from the 1967 IGF on the order of 0.8 mGal at the equator, tapering to smaller values at higher latitudes, enabling more reliable global gravity field modeling. This evolution reflects progressive alignment with observed Earth oblateness and mass distribution, establishing GRS80 as the basis for contemporary geodetic standards.[37][30]Variations and Corrections
Latitude Variations
Theoretical gravity on Earth varies systematically with latitude due to the planet's rotation and shape, resulting in a minimum value of approximately 9.780 m/s² at the equator and a maximum of 9.832 m/s² at the poles, yielding a total range of about 0.052 m/s².[34] This pattern reflects the combined influence of the Earth's oblateness and rotational effects, where gravity is weakest at the equator because of the greater distance from the planet's center of mass and the outward centrifugal acceleration, while it strengthens toward the poles where these factors diminish. The latitudinal variation arises from the Earth's oblateness, which alters the gravitational attraction through changes in distance to the center and mass distribution, and the direct centrifugal force from rotation, which reduces effective gravity most prominently at low latitudes. The functional dependence on latitude φ arises from terms proportional to sin²φ in established theoretical models, such as those in the International Gravity Formula, capturing the smooth increase in gravity from equatorial to polar regions.[30] Visual representations, including global gravity maps and latitudinal profiles, illustrate this polar-equatorial gradient as a monotonic rise in theoretical gravity values, peaking symmetrically at both poles and dipping at the equator, without incorporating local geological anomalies for pure theoretical assessment.[38] This gradient underscores the rotational dynamics of Earth, where the sin²φ term models the transition effectively over the full range of latitudes.Height and Altitude Adjustments
In theoretical gravity, adjustments for height and altitude account for the decrease in gravitational acceleration as measurements are taken above the reference ellipsoid surface, treating the intervening space as a vacuum without additional mass effects. This free-air correction isolates the geometric and potential changes due to elevation, essential for reducing gravity observations to a common datum. The baseline surface gravity g(\phi, 0) at latitude \phi, as established in latitude variation analyses, serves as the starting point for these vertical adjustments. The free-air correction formula, derived in classical physical geodesy, expresses the gravity at height h using the linear approximationg(\phi, h) = g(\phi, 0) - 3.086 \times 10^{-6} h \, \mathrm{m/s^2},
where h is in meters. This approximation is valid for h < 10 km, capturing the primary radial dilution of the gravitational field; the coefficient is an average value and varies slightly with latitude (by about 0.7%).[39] For small heights, \Delta g \approx -0.3086 \, \mathrm{mGal/m} suffices and is commonly used in practice.[30] These adjustments find key applications in aviation, where commercial flights reach altitudes of 10 km, resulting in a gravity reduction of approximately 0.031 m/s² (about 0.32% of sea-level value) that influences inertial navigation and altimetry systems. In satellite geodesy for low Earth orbit (typically 200–800 km), extensions of the free-air correction are incorporated into global models to handle potential variations, enabling precise orbit determination and Earth observation.[40][41]