Fact-checked by Grok 2 weeks ago

Geopotential spherical harmonic model

A geopotential spherical harmonic model is a mathematical framework used to represent the Earth's external gravitational potential as an infinite series expansion in terms of spherical harmonics, which are orthogonal functions that solve Laplace's equation in spherical coordinates and capture the spatial variations in the gravity field due to the planet's non-uniform mass distribution.^[1] These models express the potential V(r, \phi, \lambda) outside the Earth as V = \frac{[GM](/page/GM)}{r} \sum_{l=0}^{\infty} \sum_{m=0}^{l} \left( \frac{a_e}{r} \right)^l P_{lm}(\sin \phi) (C_{lm} \cos m\lambda + S_{lm} \sin m\lambda), where [GM](/page/GM) is the Earth's gravitational constant times mass, a_e is the equatorial radius, r, \phi, \lambda are radial distance, latitude, and longitude, P_{lm} are associated Legendre functions, and C_{lm}, S_{lm} are degree l and order m harmonic coefficients determined from observations.^[1]^[2] The expansion decomposes the geopotential into zonal harmonics (m=0), which primarily model latitudinal asymmetries like Earth's oblateness from rotation, and tesseral or sectorial terms (m > 0) that account for longitudinal and finer-scale irregularities arising from density heterogeneities in the crust and mantle.^[1] The degree l corresponds to the spatial resolution, with higher degrees resolving shorter wavelength features; for instance, the degree-2 zonal coefficient C_{20} quantifies the planet's flattening, approximately -0.000484, linking to geophysical parameters via Clairaut's theorem.^[2] Coefficients are typically normalized for numerical stability and estimated using data from satellite orbits, surface gravimetry, and altimetry, with truncation at a finite degree and order (e.g., up to 70 or higher) based on data availability and computational limits.^[1] Prominent examples include the Earth Gravitational Model 2008 (EGM2008), a high-resolution model complete to degree and order 2159 (with extensions to 2190), developed through least-squares collocation of satellite-derived gravity fields like ITG-GRACE03S and a global 5-arc-minute gravity anomaly grid incorporating terrestrial, airborne, and satellite altimetry data.^[3] EGM2008 achieves geoid height accuracies of ±5 to ±10 cm over well-surveyed regions such as the United States and Europe, surpassing predecessors like EGM96 by factors of three to six in precision and six in resolution, while integrating topographic corrections for unsampled areas.^[3] Earlier models, such as JGM-3 to degree 70, relied on satellite laser ranging and Doppler tracking for trajectory predictions.^[1] These models are fundamental in geodesy for defining the geoid—an equipotential surface approximating mean sea level—and computing gravity anomalies, vertical deflections, and orbital perturbations essential for satellite navigation, precise positioning, and geophysical interpretations of Earth's interior structure.^[2] In geophysics, they enable studies of mantle convection and isostatic compensation by isolating long-wavelength signals from short-scale crustal effects, while in oceanography, they support sea surface modeling when combined with mean sea surface heights.^[3] The orthogonal nature of spherical harmonics facilitates efficient computation of derived quantities like accelerations via the gradient of the potential, making them indispensable for both numerical simulations and analytical perturbation theory in space missions.^[1]

Fundamentals

Definition and Importance

The geopotential, denoted as W, represents the sum of the gravitational attraction potential V due to Earth's mass distribution and the centrifugal potential \Phi arising from its rotation, though in many geophysical contexts, the geopotential specifically refers to the gravitational component V excluding rotation effects.^[4] This potential satisfies Laplace's equation \nabla^2 V = 0 in the exterior region beyond the Earth's masses, allowing it to be expressed mathematically as a solution to this boundary value problem. The geopotential spherical harmonic model approximates V through an infinite series expansion in spherical harmonics, which captures the non-uniformities and irregularities in the gravity field as functions of geocentric latitude \phi and longitude \lambda. The basic form of this expansion is

V(r, \phi, \lambda) = \frac{GM}{r} \sum_{n=0}^{\infty} \sum_{m=0}^{n} \left( \frac{a}{r} \right)^n \overline{P}_{nm} (\sin \phi) ( \overline{C}_{nm} \cos m\lambda + \overline{S}_{nm} \sin m\lambda ),

where GM is the Earth's gravitational parameter, r is the radial distance from the Earth's center, a is a reference radius (typically the equatorial radius), \overline{P}_{nm} are normalized associated Legendre functions, and \overline{C}_{nm}, \overline{S}_{nm} are the normalized spherical harmonic coefficients describing the field's deviations from spherical symmetry (with the n=1 terms set to zero by convention, assuming the origin is at the center of mass).^[4] In practice, the series is truncated at a finite degree N and order M based on data resolution, enabling computational modeling of the field at various altitudes. The concept of expanding the Earth's gravitational potential in spherical harmonics originated in the early 19th century, with foundational theoretical work by figures like Carl Friedrich Gauss and Franz Ernst Neumann applying such series to geomagnetic fields in the 1830s–1840s, and George Gabriel Stokes extending similar harmonic techniques to gravity variations and the geoid in 1849.^[5]^[6] These early efforts laid the groundwork for 19th-century geodesy, where harmonic expansions were used to analyze surface gravity measurements and infer the planet's figure from limited observations. Spherical harmonics serve as the orthogonal basis functions for this representation, allowing decomposition of the potential into zonal (latitudinal), tesseral (mixed), and sectorial (longitudinal) components. Geopotential spherical harmonic models are critically important in modern geodesy and space science, as they enable precise modeling of gravitational perturbations that dominate satellite dynamics. Accurate orbit determination for missions like GRACE or GOCE relies on these models to predict trajectories with centimeter-level precision, mitigating errors from unmodeled mass distributions.^[1]^[7]^[8] In navigation, they underpin GPS accuracy by refining satellite ephemerides and transforming ellipsoidal heights to orthometric heights via the geoid, which is an equipotential surface derived from the model.^[1] Additionally, the coefficients reveal insights into Earth's interior, such as density heterogeneities from mantle convection or post-glacial rebound, aiding geophysical interpretations of planetary structure.^[1]

Spherical Harmonics Overview

Spherical harmonics Y_n^m(\theta, \phi) constitute a complete, orthogonal basis for square-integrable functions on the unit sphere, enabling the decomposition of angular dependencies in spherical coordinates. They are defined for nonnegative integers n (degree) and m (order) with |m| \leq n, as

Y_n^m(\theta, \phi) = (-1)^m \sqrt{\frac{(2n+1)(n-m)!}{4\pi (n+m)!}} P_n^m(\cos \theta) e^{i m \phi},

where P_n^m are the associated Legendre functions, and the functions for negative m follow from Y_n^{-m} = (-1)^m \overline{Y_n^m}. This normalization ensures orthonormality over the sphere:

\int_0^{2\pi} \int_0^\pi Y_n^m(\theta, \phi) \overline{Y_{n'}^{m'}(\theta, \phi)} \sin \theta \, d\theta \, d\phi = \delta_{nn'} \delta_{mm'}.

The spherical harmonics diagonalize the angular part of the Laplacian operator, with eigenvalues -n(n+1), making them ideal for solving boundary value problems in potential theory.^[9]^[10] The associated Legendre functions P_n^m(x), with x = \cos \theta and -1 \leq x \leq 1, are given by Rodrigues' formula:

P_n^m(x) = \frac{(-1)^m}{2^n n!} (1 - x^2)^{m/2} \frac{d^{n+m}}{dx^{n+m}} (x^2 - 1)^n,

and satisfy the associated Legendre differential equation

(1 - x^2) \frac{d^2 P_n^m}{dx^2} - 2x \frac{d P_n^m}{dx} + \left[ n(n+1) - \frac{m^2}{1 - x^2} \right] P_n^m = 0.

For m = 0, they reduce to the Legendre polynomials P_n(x). Key properties include orthogonality for fixed m:

\int_{-1}^1 P_n^m(x) P_{n'}^m(x) \, dx = \frac{2 (n + m)!}{(2n + 1) (n - m)!} \delta_{nn'},

and recursion relations such as

(n - m + 1) P_{n+1}^m(x) = (2n + 1) x P_n^m(x) - (n + m) P_{n-1}^m(x),

which facilitate efficient computation and relations between different degrees and orders. Additionally, P_n^{-m}(x) = (-1)^m \frac{(n - m)!}{(n + m)!} P_n^m(x), ensuring symmetry for negative orders.^[9] Spherical harmonics emerge as solutions to Laplace's equation \nabla^2 \Psi = 0 in spherical coordinates (r, \theta, \phi) through separation of variables. Assuming \Psi(r, \theta, \phi) = R(r) \Theta(\theta) \Phi(\phi), the azimuthal part yields \Phi(\phi) \propto e^{i m \phi}, the polar part leads to associated Legendre functions, and the radial part satisfies R(r) \propto r^n or r^{-(n+1)}. Thus, the general harmonic solution is \Psi_{n m}(r, \theta, \phi) = \left( A r^n + B r^{-(n+1)} \right) Y_n^m(\theta, \phi), where the degree n governs the radial decay (e.g., $1/r^{n+1} for exterior potentials) and the order m captures the longitudinal variation. For exterior fields like the geopotential, the r^{-(n+1)} term dominates beyond the source region.^[9]^[10]^[11] In geopotential modeling, fully normalized spherical harmonics are standard, where the associated Legendre functions \bar{P}_n^m incorporate factors such that

\int_0^{2\pi} \int_0^\pi \left[ \bar{P}_n^m(\cos \theta) \cos(m \phi) \right]^2 \sin \theta \, d\theta \, d\phi + \int_0^{2\pi} \int_0^\pi \left[ \bar{P}_n^m(\cos \theta) \sin(m \phi) \right]^2 \sin \theta \, d\theta \, d\phi = 4\pi,

with normalization \bar{P}_n^m(x) = \sqrt{ \frac{(2n+1) (n-m)! (2 - \delta_{0m}) }{2 (n+m)! } } P_n^m(x) for the real-valued forms used in expansions. This convention ensures unit mean-square value over the sphere, simplifying the Parseval relation and coefficient interpretation in gravity field representations.^[12]

Mathematical Formulation

Potential Expansion

The geopotential, representing the gravitational potential due to Earth's mass distribution, is commonly expanded in spherical harmonics to model its spatial variations outside the Earth's surface. This expansion provides a global representation suitable for geophysical and satellite applications, capturing the non-spherical nature of the planet through a series of harmonic terms.^[13] The standard form of this expansion for the gravitational potential V at a point with geocentric radius r, colatitude \theta, and longitude \phi is given by

V(r, \theta, \phi) = \frac{GM}{r} \sum_{n=0}^{\infty} \sum_{m=0}^{n} \left( \frac{R}{r} \right)^n \left[ C_{nm} \cos(m\phi) + S_{nm} \sin(m\phi) \right] \overline{P}_{nm}(\cos \theta),

where GM is the gravitational parameter (product of the gravitational constant and Earth's mass), R is the reference equatorial radius (typically 6378.137 km), C_{nm} and S_{nm} are the fully normalized spherical harmonic coefficients, and \overline{P}_{nm} are the fully normalized associated Legendre functions of degree n and order m. This formulation assumes a body-fixed spherical coordinate system, with \theta measured from the north pole (colatitude) and \phi as east longitude from the Greenwich meridian.^[13]^[14] In practice, the infinite series is truncated at a maximum degree N to form a finite model, balancing computational feasibility with accuracy; for instance, the Earth Gravitational Model 1996 (EGM96) uses N = 360, resolving features down to approximately 30 arcminutes. The fully normalized coefficients C_{nm} and S_{nm} ensure orthogonality and unit-mean-square amplitude for the harmonics, facilitating efficient numerical evaluation and error analysis.^[13]^[14] While this expansion describes the pure gravitational geopotential V, the full gravity field relevant to many applications includes the centrifugal potential due to Earth's rotation, yielding the total potential W = V + \Omega, where \Omega accounts for the oblate equipotential surface; however, \Omega is often treated separately in spherical harmonic models to focus on mass-induced anomalies.^[13]

Derivation from Laplace's Equation

The gravitational potential V inside the Earth satisfies Poisson's equation, \nabla^2 V = 4\pi G \rho , where G is the gravitational constant and \rho is the mass density, reflecting the presence of mass distribution.^[2] Outside the Earth, in mass-free regions, the equation simplifies to Laplace's equation, \nabla^2 V = 0 , as the density term vanishes.^[15] This distinction arises from the fundamental principles of Newtonian gravity, where the potential is sourced by mass concentrations within the planet.^[2] Boundary conditions are essential for uniqueness: the potential must approach zero as the radial distance r \to \infty, ensuring it diminishes appropriately at large distances from the Earth, and it must be continuous and match smoothly at the Earth's surface, connecting the interior and exterior solutions.^[15] These conditions stem from the physical requirement that the potential remains finite and well-behaved across the domain, with the surface serving as the interface between regions of differing equations.^[2] To solve Laplace's equation in the exterior domain, spherical coordinates (r, \theta, \phi) are employed, where r is the radial distance from the origin, \theta is the colatitude, and \phi is the longitude.^[16] The method of separation of variables assumes a product solution of the form V(r, \theta, \phi) = R(r) \Theta(\theta) \Phi(\phi), decoupling the radial dependence from the angular parts.^[15] Substituting into Laplace's equation yields three ordinary differential equations: one for the azimuthal part \Phi, which implies periodicity in \phi with separation constant m^2, leading to solutions \Phi(\phi) \propto e^{\pm i m \phi} for integer m \geq 0; one for the polar part \Theta, resulting in associated Legendre functions P_n^m(\cos \theta) after referencing their orthogonality properties; and one for the radial part.^[16] The radial equation takes the form r^2 R''(r) + 2r R'(r) - n(n+1) R(r) = 0, where n is the degree of separation, an integer ensuring angular finiteness.^[17] The general solution is R(r) = A r^n + B r^{-(n+1)}. For the exterior domain, the condition V \to 0 as r \to \infty requires discarding the growing r^n term, leaving R(r) \propto r^{-(n+1)}.^[15] This decaying radial dependence ensures the potential's physical behavior far from the source.^[2] The derivation assumes a reference frame centered at the Earth's center of mass, with the z-axis aligned along the rotation axis to incorporate axisymmetry, and approximates the Earth using a spherical reference surface for the boundary, though real applications often adjust for the ellipsoidal shape via perturbations.^[2] Superposition of solutions for all n and m then yields the general exterior form, matching boundary values at the surface to determine coefficients.^[15]

Harmonic Coefficients

Types and Normalization

In geopotential spherical harmonic models, the coefficients are classified into three main types based on the relationship between the degree n and the order m: zonal, tesseral, and sectorial harmonics. Zonal harmonics occur when m = 0, describing axisymmetric variations that depend only on latitude and are independent of longitude, capturing global features like the Earth's oblateness. Tesseral harmonics arise for $0 < m < n, representing general perturbations that vary in both latitude and longitude, often associated with non-axisymmetric mass distributions such as continental structures. Sectorial harmonics correspond to m = n, exhibiting strong longitude dependence with minimal latitudinal variation, typically modeling east-west gravitational anomalies.^[18]^[19] The zonal coefficients are denoted as C_{n0}, which are related to the conventional zonal harmonics J_n by J_n = -C_{n0} for n \geq 2, with J_2 being the dominant term representing equatorial flattening. For non-zonal terms, the coefficients include C_{nm} for the cosine components and S_{nm} for the sine components, reflecting the longitudinal variations in the gravitational potential. These classifications ensure physical interpretability, as zonal terms dominate low-degree models while tesseral and sectorial terms account for higher-order regional effects.^[19]^[20] To promote orthogonality, numerical stability, and consistent energy distribution across the sphere, geopotential models employ specific normalization schemes for the associated Legendre functions and spherical harmonics. Fully normalized (or 4π-normalized) coefficients, standard in modern global models like EGM2008, scale the associated Legendre functions as

\bar{P}_{nm}(\cos \theta) = \sqrt{(2 - \delta_{m0})(2n + 1) \frac{(n - m)!}{(n + m)!}} \, P_{nm}(\cos \theta),

where \delta_{m0} = 1 if m=0 and 0 otherwise, and P_{nm} are the unnormalized functions; this ensures the spherical harmonics satisfy \int_{S^2} |\bar{Y}_{nm}|^2 \, d\Omega = 4\pi. Schmidt semi-normalized coefficients, historically used in some gravity and magnetism applications, omit the (2n + 1) factor:

\hat{P}_{nm}(\cos \theta) = \sqrt{(2 - \delta_{m0}) \frac{(n - m)!}{(n + m)!}} \, P_{nm}(\cos \theta),

resulting in integrals of $4\pi for m > 0 and $2\pi for m = 0. Conversion between fully normalized \bar{C}_{nm}, \bar{S}_{nm} and Schmidt semi-normalized \hat{C}_{nm}, \hat{S}_{nm} coefficients involves scaling by \sqrt{2n + 1}, while unnormalized to fully normalized requires division by the full factor \sqrt{(2 - \delta_{m0})(2n + 1) \frac{(n - m)!}{(n + m)!}}. For zonal terms, the relation holds similarly, with J_n typically expressed in unnormalized form but convertible accordingly.^[19]^[21]^[19] Normalization choices impact error propagation in computations, particularly for high-degree models where unnormalized functions can lead to overflow or underflow in associated Legendre evaluations, amplifying numerical instabilities by orders of magnitude (e.g., up to 200 for degree 100). Fully normalized schemes mitigate this by bounding coefficient magnitudes, reducing relative errors in potential evaluations by factors related to the normalization constant, and are preferred for global models to maintain orthogonality and facilitate power spectral analysis without additional scaling. Schmidt semi-normalization, while computationally simpler for certain recursions, can propagate larger errors in tesseral and sectorial terms due to uneven energy distribution, making fully normalized the default for high-resolution geopotential applications.^[22]^[19]

Determination from Observations

The determination of spherical harmonic coefficients for geopotential models relies primarily on observations from satellite gravimetry missions, such as the Gravity Recovery and Climate Experiment (GRACE) launched in 2002 and the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) operational from 2009 to 2013. GRACE utilized inter-satellite ranging to measure time-variable gravity changes, enabling the recovery of monthly spherical harmonic coefficients up to degree and order 60 or higher, which capture mass redistributions like those from ice melt and groundwater depletion.^[23] GOCE, on the other hand, employed gravity gradiometry to sense the spatial derivatives of the gravitational field, providing enhanced resolution for static high-degree coefficients up to degree 250, particularly sensitive to the Earth's crust and upper mantle structure.^[24] Complementary data sources include satellite altimetry missions, which derive sea surface heights convertible to gravity anomalies via the remove-restore technique, and terrestrial surface gravity measurements from gravimeters, which offer dense coverage over land but require integration to mitigate regional biases.^[25] The core method for estimating these coefficients is least-squares adjustment, where observed data—such as orbital perturbations, gravity gradients, or potential differences—are fitted to the spherical harmonic expansion to solve for the coefficients C_{nm} and S_{nm}. This involves minimizing the residuals between modeled and observed potentials, often incorporating energy conservation principles for satellite orbits or direct gradient observations, with constraints like Kaula's rule to stabilize higher-degree terms.^[26] For GRACE-derived models, this approach produced monthly solutions from 2002 to 2017, revealing time-variable gravity signals with accuracies improving to about 1 cm equivalent water height at large scales, though limited by the mission's spatial resolution of approximately 300-400 km.^[23] These adjustments typically process data in short orbital arcs (days to months) to isolate gravitational signals from non-gravitational forces. As of October 2025, GRACE-FO continues to provide monthly gravity solutions, with its science team meeting held in Boulder, Colorado.^[27] Key challenges in this determination include temporal aliasing from unmodeled short-period signals in short arcs, which can introduce north-south striping artifacts in the recovered fields, and correlations between geopotential signals and atmospheric or oceanic mass variations that require precise de-aliasing models.^[28] For instance, incomplete modeling of atmospheric tides in GRACE data can alias into low-degree harmonics, reducing signal fidelity by up to 20% in mid-latitudes.^[29] Post-2020 advances have extended these methods through the Swarm mission's magnetometry-derived accelerations, which, when combined via least-squares with GRACE-FO data, bridge gaps in time-variable gravity monitoring and enhance low-degree coefficient stability up to degree 20.^[30] GRACE-FO, launched in 2018, continues the inter-satellite ranging legacy with improved microwave instruments, enabling monthly solutions that integrate with Swarm for reduced aliasing errors and better separation of geophysical signals, achieving noise levels below 0.5 cm water height equivalent in recent combined products.^[31] Future integrations with planned missions, such as NASA's GRACE-C (scheduled for launch in 2028) and ESA's MAGIC (Mass-change And Geosciences International Constellation) constellation (targeted for the early 2030s), promise even higher temporal resolution for dynamic geopotential modeling.^[32]^[33]

Dominant Terms and Effects

Zonal Harmonics

Zonal harmonics represent the axisymmetric components of the Earth's geopotential, arising from mass distributions that are symmetric about the rotation axis and thus independent of longitude (m=0 terms in the spherical harmonic expansion). These are characterized by the coefficients C_{n0}, with the zonal harmonics defined as J_n = -C_{n0} for even degrees n = 2, 4, \dots, while odd-degree coefficients vanish due to the Earth's equatorial symmetry.^[20] The dominant zonal harmonic is the degree-2 term J_2, with a normalized value of approximately +4.84 \times 10^{-4}, which primarily accounts for the Earth's oblateness and constitutes the leading contribution to the non-spherical geopotential.^[20] Physically, J_2 originates from the equatorial bulge induced by the planet's rotation, which flattens the Earth at the poles; higher-order even zonal harmonics J_4, J_6, \dots reflect additional axisymmetric mass anomalies from internal density variations and surface topography.^[34] This J_2 term exerts significant dynamical effects, including the precession of the Earth's equator under lunisolar torques enhanced by oblateness and the regression of the nodes in satellite orbits, where it induces a westward drift in the ascending node longitude proportional to the orbit's inclination and semi-major axis.^[35] The perturbing potential due to J_2 is approximated as

\Delta V \approx \frac{GM}{r} \left( \frac{R}{r} \right)^2 J_2 P_2 (\sin \theta),

where GM is Earth's gravitational parameter, [r](/page/R) is the geocentric distance, [R](/page/R) is the equatorial reference radius, \theta is the geocentric latitude, and P_2 (\sin \theta) = \frac{1}{2} (3 \sin^2 \theta - 1) is the associated Legendre function of degree 2.^[36]

Tesseral and Sectorial Harmonics

Tesseral harmonics refer to the components of the geopotential spherical harmonic expansion where the azimuthal order m satisfies $0 < m < n, with n being the degree; these terms introduce non-axisymmetric, longitude-dependent perturbations to the gravitational field. Unlike zonal harmonics, which vary only with latitude, tesseral harmonics produce east-west gravity variations that reflect the Earth's irregular mass distribution, such as deviations from rotational symmetry. A prominent example is the degree-2, order-2 coefficient C_{2,2}, which quantifies the equatorial triaxiality arising from the planet's oblate shape and internal density contrasts. The magnitude of C_{2,2} is approximately $2.44 \times 10^{-6}, an order of magnitude smaller than the zonal J_2 term but crucial for high-precision applications due to its influence on orbital dynamics.^[20]^[2] Sectorial harmonics, characterized by m = n, exhibit nodal lines along meridians and capture variations primarily with latitude, emphasizing circumferential mass asymmetries. These terms are particularly sensitive to long-wavelength geophysical features, such as the density contrasts between oceanic basins and continental crusts, which generate broad-scale gravity anomalies observable in global models. For instance, sectorial coefficients at low degrees contribute to the undulation of the geoid over major land-ocean boundaries, aiding interpretations of mantle convection and crustal thickness.^[2]^[37] The primary dynamic effects of tesseral and sectorial harmonics manifest as perturbations in satellite trajectories, notably along-track accelerations in low-Earth orbits, where these terms induce velocity changes on the order of micrometers per second over short arcs. In geostationary satellites, resonance with dominant tesseral harmonics—such as those commensurate with the Earth's rotation—leads to periodic longitude drifts and libration, necessitating east-west stationkeeping to maintain fixed positions relative to the surface. These resonances arise when the orbital period aligns with multiples of the sidereal day, amplifying the otherwise small perturbations into secular trends.^[38]^[39] Time-variable tesseral harmonics further complicate modeling, with the M2 tidal term serving as a key example; this principal lunar semidiurnal constituent generates periodic geopotential variations through ocean and solid Earth tides, represented as degree-2, order-2 tesseral coefficients with amplitudes around $10^{-10} that modulate over tidal cycles. Such effects are incorporated into modern geopotential models to account for sub-monthly gravity fluctuations impacting precise orbit propagation.^[20]

Computational Applications

Recursive Algorithms

Recursive algorithms are crucial for efficiently evaluating the spherical harmonic expansion of the geopotential, particularly for high-degree models where direct computation of each term would be prohibitive. These methods leverage recurrence relations to compute associated Legendre functions P_n^m(\sin \phi) and their derivatives iteratively, reducing the need for independent evaluations and enabling real-time applications in orbit determination and geophysical modeling. The core approach involves three-term recurrence relations that relate P_n^m to values at lower degrees n and orders m, allowing sequential computation from initial seed values.^[40] One widely used recursion for associated Legendre functions proceeds along increasing n for fixed m, using the relation:

(n - m + 1) P_{n+1}^m (x) = (2n + 1) x P_n^m (x) - (n + m) P_{n-1}^m (x),

with initial conditions P_0^0(x) = 1 and P_1^0(x) = x. These formulas enable the generation of the full grid of P_n^m up to degree N starting from zonal terms (m=0). However, forward recursion in m can introduce numerical instability due to alternating signs and growing magnitudes, particularly for high n. To mitigate this, sectorial-to-zonal recursion is employed, where functions are first computed along the sectorial diagonal (m = n) using:

P_n^n (x) = (-1)^n (2n - 1)!! (1 - x^2)^{n/2},

and then recursed downward in m for each fixed n via:

P_n^{m-1} (x) = \frac{(2n - 1) x P_n^m (x) - (n + m) P_n^{m+1} (x)}{n - m + 1}.

This backward recursion in m enhances stability by avoiding underflow in low-order terms and overflow in high-order ones, especially for arguments near the poles where x \approx \pm 1.^[22]^[41] For applications involving gridded geopotential fields, such as global gravity anomaly maps, collocation methods precompute the normalized associated Legendre functions and spherical harmonics at discrete latitude-longitude grid points. These values are stored in lookup tables, allowing rapid summation of the potential or its derivatives at any grid location without on-the-fly recursion. This approach is particularly effective for fixed-resolution outputs, where the precomputation overhead is amortized over multiple evaluations, though it requires careful management of storage for high-degree models.^[42] The overall computational complexity of evaluating the geopotential via recursive algorithms up to degree N is O(N^2), arising from the (N+1)^2 harmonic terms and the linear cost per term using recursion. Optimizations tailored for spacecraft applications, such as those addressing the Gibbs phenomenon—oscillatory artifacts near sharp gravitational features from truncated series—employ spectral filtering or localized basis adjustments to reduce ringing without increasing complexity significantly. For instance, Gibbs' algorithm resolves these artifacts by reconstructing spectrally accurate approximations from the first (N+1)^2 coefficients, improving accuracy in orbit propagation near equatorial discontinuities.^[43] Numerical errors in recursive computations primarily stem from rounding in floating-point arithmetic, exacerbated by the large dynamic range of associated Legendre functions (spanning orders of magnitude from $10^{-300} to $10^{300} for high n). These errors accumulate in forward recursions, leading to loss of precision in tesseral terms. Mitigation strategies include computing unnormalized intermediates throughout the recursion and applying normalization factors only at the final step, which preserves relative accuracy and avoids intermediate overflow or underflow. Extended-precision arithmetic or rescaling during recursion further bounds error growth to machine epsilon levels.^[44]^[45]

Orbit Propagation Techniques

In orbit propagation, the geopotential spherical harmonic model contributes to the equations of motion through the gravitational acceleration term, where the second derivative of the position vector \ddot{\mathbf{r}} is given by \ddot{\mathbf{r}} = -\nabla V + \mathbf{a}_{\text{other}}, with V representing the gravitational potential expanded in spherical harmonics and \mathbf{a}_{\text{other}} accounting for non-gravitational forces such as atmospheric drag or solar radiation pressure.^[46] The gradient \nabla V is computed from the potential V = \frac{\mu}{r} \sum_{n=2}^{N} \sum_{m=0}^{n} \left( \frac{R_e}{r} \right)^n \overline{P}_{nm}(\sin \phi) (\overline{C}_{nm} \cos m\lambda + \overline{S}_{nm} \sin m\lambda), where \mu is the Earth's gravitational constant, R_e is the reference radius, \overline{P}_{nm} are normalized associated Legendre functions, and \overline{C}_{nm}, \overline{S}_{nm} are the normalized harmonic coefficients.^[46] This formulation allows the geopotential perturbations to be integrated alongside the central two-body acceleration for accurate trajectory prediction over extended periods. Numerical integrators are essential for solving these perturbed differential equations, with explicit Runge-Kutta methods, such as the adaptive embedded Fehlberg schemes of orders 7 and 8, widely employed due to their balance of accuracy and computational efficiency in Cartesian coordinates.^[46] For instance, the state vector \mathbf{s} = [\mathbf{r}, \mathbf{v}] evolves via \dot{\mathbf{s}} = \mathbf{f}(t, \mathbf{s}), where step sizes are adjusted based on relative tolerances (e.g., $10^{-11}) to meet positional accuracies on the order of meters over multi-day propagations.^[46] Symplectic integrators, such as those based on Deprit's radial intermediary or Gauss-Legendre collocation, offer advantages for long-term stability in Hamiltonian systems like orbital mechanics, preserving energy and phase space structure better than non-symplectic Runge-Kutta methods, particularly when geopotential perturbations induce secular drifts. These methods are preferred in scenarios requiring simulations over years, reducing artificial dissipation in highly eccentric or resonant orbits. Perturbation approaches like Encke's method enhance efficiency by separating the Keplerian reference orbit from small deviations, formulating the equations for the deviation vector \mathbf{S} as \ddot{\mathbf{S}} + \frac{\mu \mathbf{S}}{r^3} = \mathbf{a}_p - \frac{\mu (\mathbf{r} + \mathbf{S})}{|\mathbf{r} + \mathbf{S}|^3} + 2 \boldsymbol{\omega} \times \dot{\mathbf{S}} + \boldsymbol{\omega} \times (\boldsymbol{\omega} \times \mathbf{S}) + \dot{\boldsymbol{\omega}} \times \mathbf{S}, where \mathbf{r} is the reference position, \boldsymbol{\omega} is the reference angular velocity, and \mathbf{a}_p includes geopotential accelerations from models like GEM-9.^[47] This technique is particularly effective for low-Earth orbits where perturbations are modest relative to the central force, allowing reinitialization of the reference orbit when deviations grow, thus maintaining numerical stability without the full integration cost of Cowell's method.^[47] High-degree spherical harmonic models pose computational challenges in orbit propagation, leading to the use of truncated expansions where the maximum degree N is selected based on altitude and required precision—typically lower (e.g., N \approx 20-50) for real-time onboard applications to minimize processing time, versus full-degree models (e.g., N > 200) for post-mission analysis. Truncation errors, assessed via omission and commission uncertainties, can be quantified using metrics like the Mahalanobis distance to ensure positional errors remain below mission thresholds, such as 1 km for lunar orbits at 500 km altitude. A representative case is the GPS constellation, where the zonal harmonic J_2 induces a nodal precession that is compensated by setting orbital periods approximately 4 seconds faster than half-sidereal to achieve daily ground track repeatability within ±2° longitude. Tesseral harmonics further perturb the semimajor axis through resonances, causing secular variations up to several meters and daily oscillations on the order of 5 seconds in repeat times, necessitating precise modeling for maintaining the 20,200 km altitude orbits' accuracy to within centimeters for navigation applications.

Specific Models

Historical Developments

The theoretical foundations of geopotential spherical harmonic models were laid in the mid-20th century by Weikko A. Heiskanen and Helmut Moritz, who advanced the mathematical framework for representing the Earth's gravity field using spherical harmonics in their seminal work on physical geodesy during the 1950s.^[48] Their contributions emphasized the expansion of the gravitational potential into harmonic coefficients derived from surface gravity measurements, providing a basis for subsequent global models.^[49] Early geopotential models relied primarily on terrestrial gravity data, with Richard H. Rapp developing a high-degree expansion to spherical harmonic degree 360 in 1986 using surface observations to capture finer gravitational structures.^[50] This model marked a significant advancement in resolution from prior low-degree efforts, though it depended heavily on a priori assumptions for higher-order coefficients due to incomplete global coverage.^[51] The advent of satellite observations in the 1970s ushered in the satellite era of geopotential modeling, beginning with the GEM10 model released in 1977 by NASA's Goddard Space Flight Center, complete to degree and order 30 and incorporating tracking data from satellites like Geos-3.^[52] Concurrently, the GRIM series emerged from a German-French collaboration starting in the mid-1970s, with GRIM1 (1976) and GRIM2 utilizing optical and laser ranging data from geodetic satellites to refine low-degree harmonics.^[53] These models improved accuracy over surface-only approaches by integrating orbital perturbations, yet remained limited to moderate degrees due to data sparsity.^[54] By the 1990s, pre-GRACE models achieved greater integration of data sources, exemplified by EGM96 (1996), a joint NASA and National Imagery and Mapping Agency effort complete to degree 360 that combined satellite tracking, altimetry, and over 2 million terrestrial gravity anomalies for enhanced global consistency.^[55] This model represented a milestone in blending disparate observations to reduce long-wavelength errors, though its determination of coefficients still referenced techniques from earlier harmonic analyses.^[3] Historical models faced inherent limitations, including low resolution for short-wavelength gravity features because satellite altitudes attenuated high-degree harmonics, and reliance on smoothing assumptions or incomplete datasets that introduced uncertainties in regional anomalies.^[54] These constraints underscored the need for denser, space-based gravity observations to surpass the capabilities of pre-2000 developments.^[52]

Modern and High-Resolution Models

Modern high-resolution geopotential models have advanced significantly since the early 2000s, leveraging satellite missions such as GRACE, GOCE, and their successors to achieve spatial resolutions approaching 10 km. The Earth Gravitational Model 2008 (EGM2008), developed by the National Geospatial-Intelligence Agency (NGA), represents a landmark in this progression, expanding to spherical harmonic degree and order 2190 through the integration of GRACE satellite gravimetry, terrestrial gravity anomalies, and satellite altimetry data. This model provides global gravity field information equivalent to a 5 arcminute grid, enabling applications in precise geodesy and oceanography with uncertainties below 20 cm for geoid heights in many regions.^[3] Building on EGM2008, the Earth Gravitational Model 2020 (EGM2020) is under development by the NGA, planned for release around 2025, and will incorporate additional datasets from GOCE gradiometry and extended GRACE observations, maintaining the high degree 2190 while aiming for improved accuracy in medium- to high-frequency components through refined least-squares collocation techniques. As of November 2025, EGM2020's development emphasizes ellipsoidal harmonics for better representation of Earth's non-spherical shape, with evaluations ongoing to confirm enhancements in global consistency over its predecessor.^[56]^[21] These static models serve as baselines but are complemented by dynamic updates to capture temporal variations. Time-variable geopotential models have become essential for monitoring Earth's mass redistribution, particularly through the GRACE Follow-On (GRACE-FO) mission, operational since 2018 and providing monthly spherical harmonic solutions up to degree 60 as of 2025. GRACE-FO data reveal seasonal fluctuations in continental water storage and ice sheet mass loss, such as annual changes exceeding 100 Gt in Greenland ice melt, by isolating non-tidal signals from inter-satellite ranging and accelerometer measurements. These models, processed by centers like the Center for Space Research (CSR), achieve geoid height sensitivities of 1-2 cm for basin-scale hydrology.^[57] The Gravity Observation Combination (GOCO) series exemplifies combined satellite-only approaches, with GOCO06s deriving a static field up to degree 300 from over 15 years of GRACE, GOCE, and supplementary GPS data, while incorporating temporal constraints for improved signal-to-noise ratios in the 100-300 km wavelength band. Subsequent releases like GOCO05c extend to degree 720 by blending altimetric and terrestrial inputs, offering a versatile framework for regional refinements.^[58] Looking ahead, proposed next-generation gravity missions, including concepts reaching Technology Readiness Level 6 by 2025, aim to combine multi-satellite arrays with laser interferometry for resolutions below 5 km and continuous temporal coverage. Despite these improvements, uncertainties persist: commission errors for low-degree coefficients have decreased to approximately $10^{-10} in normalized units, enhancing long-wavelength accuracy, yet polar regions suffer from incomplete orbital coverage, leading to omission errors up to 50% higher than equatorial zones in high-resolution models.^[59]

References

[1]
None
### Summary of Spherical Harmonic Representation of the Gravity Field Potential
[2]
[PDF] Chapter 2 - The Earth's Gravitational field
system spherical harmonics are the general solutions of Laplace's equation. aSurface spherical harmonics are at the surface of a sphere what a Fourier series is ...
[3]
The development and evaluation of the Earth Gravitational Model ...
Apr 19, 2012 · EGM2008 is a spherical harmonic model of the Earth's gravitational potential, developed by a least squares combination of the ITG-GRACE03S gravitational model.<|control11|><|separator|>
[4]
[PDF] Definition of Functionals of the Geopotential and Their Calculation ...
To compare or combine gravity anomalies or gravity disturbances derived from terrestrial measurements with those derived from gravity field models in spherical ...
[5]
Global spherical harmonic analysis by least-squares and numerical ...
It is also shown that the two-step formulation of global spherical harmonic computation was applied already by Neumann (1838) and Gauss (1839). Computational ...
[6]
Earth's gravity field | McGraw Hill's AccessScience
For the potential outside the Earth, an expansion into an infinite series (of spherical harmonics, discussed later in this article) is useful; the principal ...
[7]
[PDF] The Spherical Harmonics
ℓ (θ, φ) = ℓ(ℓ + 1)Y m ℓ (θ, φ) . That is, the spherical harmonics are eigenfunctions of the differential operator L2, with corresponding eigenvalues ℓ(ℓ + 1), ...
[8]
[PDF] Spherical Harmonics
Spherical harmonics are the Fourier series for the sphere. These functions can are used to build solutions to Laplace's equation and other differential ...
[9]
[PDF] Boundary Value Problems in Electrostatics II
Dec 23, 2000 · In particular, the first topic is the separation of variable method in spherical polar coordinates. 1 Laplace Equation in Spherical Coordinates.
[10]
[PDF] GRACE Gravity Model Suite GGM05
The coefficients for GGM05 are normalized according to the so-called “fully-normalized” convention, where the squared norm of a spherical harmonic over a unit ...
[11]
[PDF] egm96-geoid-use-global-geopotential-model.pdf
Spherical harmonic models of the Earth's external gravitational potential ("geopotential models") contain no explicit information about which level surface, out ...Missing: explanation | Show results with:explanation
[12]
[PDF] GRACE 327-734 (CSR-GR-03-01) Gravity Recovery and Climate ...
Dec 1, 2003 · Though the spherical harmonic expansion of the geopotential requires an infinite series of harmonics, practicality dictates that the ...
[13]
[PDF] Physical Geodesy - Universität Stuttgart
Jun 15, 2006 · i) write down Laplace's equation in spherical coordinates, ... This is equation is known as the fundamental equation of physical geodesy.
[14]
[PDF] Spherical Harmonics and Related Topics - Colorado State University
We look for solutions of Laplace's differential equation, which is. ∇2S = 0 ... The ∇2 operator can be expanded in spherical coordinates as: ∇2S = 1.
[15]
[PDF] ∫ d2 ˆs ∇1Ym
Solution of Laplace's Equation in Spherical Harmonics. This is the application of spherical harmonics that we need for potential theory. We consider the ...
[16]
[PDF] ICGEM – 15 years of successful collection and distribution of global ...
The images show one specific surface spherical harmonic of degree l and order m such as (a) tesseral (l = 9, m = 4), (b) sectorial (l = 9, m = 9), and (c) zonal ...
[17]
SHTools: Tools for Working with Spherical Harmonics - AGU Journals
May 12, 2018 · SHTools is a time and user-tested open-source archive of both Fortran 95 and Python routines for performing spherical harmonic analyses.Missing: sectorial fully
[18]
[PDF] 6 Geopotential (01 February 2018) - IERS Conventions Centre
Feb 1, 2018 · The IERS recommends the EGM2008 model as the conventional geopotential model, complete to degree and order 2159, with scaling parameters GM⊕ ...
[19]
World Geodetic System 1984 (WGS 84) - NGA - Office of Geomatics
The EGM2008 is provided as a set of normalized, geopotential coefficients complete to degree and order 2159, and contains additional spherical harmonic ...
[20]
[PDF] On the Numerical Problems of Spherical Harmonics - geodaesie.info
Normalized spherical harmon ics involve normalization factors for all degrees and orders. In order to preserve Eq. (1) valid, geopotential coefficients have to ...Missing: fully | Show results with:fully<|control11|><|separator|>
[21]
High‐Resolution GRACE Monthly Spherical Harmonic Solutions
Dec 2, 2020 · This study proposes a new regularization method. Transforming spatial constraints from GRACE-based filtered mass changes into the spectral domain and imposing ...Missing: formula | Show results with:formula
[22]
Gravity Field Model Determination Based on GOCE Satellite Point ...
Spherical harmonics triangle of the estimated geopotential coefficients when compared with EIGEN-6C4. To assess the performance of the recovered gravity ...
[23]
A High-Resolution Earth's Gravity Field Model SGG-UGM-2 from ...
This paper focuses on estimating a new high-resolution Earth's gravity field model named SGG-UGM-2 from satellite gravimetry, satellite altimetry, ...<|control11|><|separator|>
[24]
[PDF] Implementation of Parallel Least-Squares Alorithms for Gravity Field ...
The DMA method can solve for all geopotential coefficients complete to spherical harmonic degree 120 in roughly 30 minutes using 24 CPUs. The serial ...
[25]
Treatment of temporal aliasing effects in the context of next ...
Jul 27, 2017 · Temporal aliasing effects have a large impact on the gravity field accuracy of current gravimetry missions and are also expected to dominate ...
[26]
(PDF) Impact of short period, non-tidal, temporal mass variability on ...
Aug 6, 2025 · De-aliasing done with approximate models gave the greatest reduction in aliasing error for the mid-degrees and higher. For the atmosphere, the ...
[27]
[PDF] Description of the multi-approach gravity field models from Swarm ...
Jun 22, 2020 · Using GRACE data for comparison, we demonstrate that the geophysical signal in the Swarm GFMs is largely restricted to spherical harmonic ...
[28]
Combined monthly GRACE-FO gravity fields for a ... - Oxford Academic
The GRACE-FO combination is used to derive global grids of groundwater storage anomalies. To meet the user requirements and achieve optimal signal-to-noise ...
[29]
Bridging the gap between GRACE and GRACE Follow-On by ...
Sep 13, 2024 · We present an intermediate technique bridging the gap between the two missions allowing (1) for a continued and uninterrupted time series of mass observations.2.1 Hl-Sst Gravity Field... · 2.1. 1 Kinematic Orbit... · 3.1. 2 Time Series For...<|control11|><|separator|>
[30]
[PDF] Evaluation of the EGM2008 Geopotential Model for Egypt
Table 1 shows that EGM2008 model fits best to the point free-air data in view of the smallest standard deviation of the differences between the point free-air ...Missing: convention | Show results with:convention
[31]
Spherical Harmonic coefficients of Degree 2 - GRACE Tellus - NASA
The spherical harmonic of degree 2 and order 0 - C(2,0) - is due to the flattening of the Earth. Its technical name is 'Earth's dynamic oblateness'.
[32]
J2 Perturbation - a.i. solutions
The coefficients of each term in this series is described as Jk, of which J2, J3, and J4 are called "zonal coefficients." However, J2 is over 1000 times larger ...
[33]
Zonal Harmonic Gravity Model - Simulink - MathWorks
The Zonal Harmonic Gravity Model block calculates the zonal harmonic representation of planetary gravity at a specific location based on planetary ...Description · Ports · Parameters
[34]
[PDF] Survey and comparative analysis of current geophysical models
If the spherical harmonic is a function of two'variables, such as latitude and longitude (rather than only one, the lati tude) , it involves the so-called ...
[35]
Variation of spherical harmonic power as a function of harmonic ...
Oct 27, 2006 · For the continents, induced magnetization is probably the cause, while for the ocean basins, remanent magnetization is likely to be the cause, ...
[36]
Perturbations of satellite orbits by tesseral harmonics in the earth's ...
The effects of tesseral harmonics, which express the variation of the Earth's gravity with longitude, are investigated. Long-period variations due to these ...
[37]
East–west stationkeeping of satellite orbits with resonant tesseral ...
This class of orbit is typically in resonance with multiple Earth tesseral harmonics. Depending on the selected resonance (i.e. orbits with a 24 h period, 12 h ...
[38]
Recursion formulas of Legendre functions for use with nonsingular ...
The use of the classical and associated Legendre functions provides only one reasonably efficient recursion formula for the computation of the Legendre ...Missing: Holiner | Show results with:Holiner
[39]
[PDF] Normalization and Implementation of Three Gravitational ...
As part of Legendre's approach to the problem, he introduced a family of functions that eventually became known as Associated Legendre Functions (ALFs).
[40]
[PDF] Numerical Methods for Harmonic Analysis on the Sphere - DTIC
ficient algorithm for numerical analysis of data regularly sampled on the sphere. The algorithms presented In section 1 for implementing numerical quad-.
[41]
[PDF] EFFICIENT ON-ORBIT SINGULARITY-FREE GEOPOTENTIAL ...
The computation of a spherical harmonic geopotential can be very computationally intensive to run, even when using a truncated model. This can be due to the ...
[42]
Recurrence relations for fully normalized associated legendre ...
The definite integral of the associated Legendre functions by recursion is elaborated in detail in Fukushima (2012). An alternative recursion formula is derived ...
[43]
[PDF] DRAFT GSL Technical Report #1 - GNU.org
2.2 Spherical harmonic normalization. The spherical harmonic normalized ALFs Y m l (x) are defined with. Am l = r. 2l + 1. 4π. (9) so that. Y m l (x) = s. 2l + ...
[44]
https://www.researchgate.net/publication/292464284_Recurrence_relations_for_fully_normalized_associated_legendre_functions_and_their_derivatives_and_integrals
[45]
[PDF] thesis - DTIC
The Encke Method assumes that the spacecraft is traveling in a conical path called osculating orbit. This osculating orbit is determined by the position and ...
[46]
[PDF] Methods of Physical Geodesy - The Ohio State University
This report was prepared by Dr. W. A. Heiskanen, Director, and. Dr. H. Moritz, Research Associate, of the Institute of Geodesy,.
[47]
Full text of "Heiskanen Moritz 1967 Physical Geodesy"
HEISKANEN HELMUT MORITZ Contents 1 Fundamentals of Potential Theory 4/1-1. ... The geopotential number C is measured in geopotential units (g.p.u.), where ...
[48]
[PDF] The Representation of the Earth's Gravitational Potential in a ... - DTIC
In 1981 Rapp carried out a second expansion to degree 180 using a more complete set of a priori potential coefficients than used in the 1978 solution. Other ...
[49]
The development and analysis of geopotential coefficient models to ...
Dec 10, 1990 · History of Geophysics · Marine Geology ... The development and analysis of geopotential coefficient models to spherical harmonic degree 360.
[50]
Gravity model improvement using Geos 3 (GEM 9 and 10) - Lerch
Jul 10, 1979 · GEM 10 is a combination solution containing a global set of surface gravity anomalies along with the data in GEM 9. Radial errors of Geos 3 for ...
[51]
The GRIM 2 earth gravity field model - NASA ADS
A new geopotential model, GRIM 2, was computed, using most of the optical and laser observations made on major geodetic satellites as well as all ...
[52]
[PDF] THE EVOLUTION OF EARTH GRAVITATIONAL MODELS USED IN ...
This article summarizes the history of the tracking and instrumentation systems used, discusses the limitations and constraints of these systems, and reviews ...Missing: MacCullagh 1830
[53]
[PDF] The Development of the Joint NASA GSFC and the National Imagery ...
... Development of the Joint NASA GSFC and the. National Imagery and Mapping Agency (NIMA). Geopotential Model EGM96. F.G. Lemoine, Laboratory for Terrestrial ...
[54]
https://www.jhuapl.edu/Content/techdigest/pdf/V15-N04/15-04-Vetter.pdf
[55]
GRACE-FO Level-2 Monthly Geopotential Spherical Harmonics ...
FOR EXPERT USE ONLY. This dataset contains estimates of the total month-by-month geopotential of the Earth, derived from the Gravity Recovery and Climate ...
[56]
GOCO06s – a satellite-only global gravity field model
Jan 27, 2021 · GOCO06s is a satellite-only global gravity field model based on 15 years of satellite data, including temporal variations.
[57]
Short‐Period Mass Variations and the Next Generation Gravity Mission
Jan 21, 2025 · Those efforts are expected to demonstrate a Technology Readiness Level 6 (TRL 6) in 2025 (Dávila Álvarez et al., 2022). The Gravitational ...