The launch of artificial satellites (as early as in 1957), specifically the launch of the first laser tracked satellite, BeaconB, in 1964, has provided data sets which have allowed researchers to probe the long to medium components of the gravitational field of the Earth. In particular, observational data recorded at satellite laser ranging tracking stations have since been used to develop models that quantify the global longwavelength and mediumwavelength gravity field of the Earth. Currently, literature reviewing gravity field models with geophysical applications is scarce and not up to date. The most recent review paper was published more than a decade ago. In the interim, there has been an unprecedented increase in gravity field modelling, which can be attributed to the deployment of new and dedicated satellite missions. As a result, a number of existing geopotential models have been improved and new models have been developed. Each of these models differs in accuracy and spatialtemporal scale. This review extends the earlier review of gravity field models, by incorporating uptodate research efforts in geopotential modelling with geophysical applications in oceanography, hydrology, geodesy and solid Earth science.
The concept of the Earth’s gravity field is often described through gravity or gravitational potential. Alternatively, the definition of gravity can be viewed in terms of the cause and effect of
gravity. To this end, gravity could be described as an attraction that causes acceleration amongst objects on or near the Earth’s surface. The gravitational potential of the Earth is the quantity
of energy that is associated with the position of a unit mass in the gravitational field of the Earth.^{1} Earth itself is a complex dynamic system driven by many geophysical processes. These include the coupled atmosphere–ocean system, varying mass distribution of ice and the isostatic correction from
the glacial loading of the last Ice Age and mobile tectonic plates.^{2} In addition, internal mass distribution is often controlled by thermal convection of the core mantle.^{2} Some of
the geophysical processes taking place within Earth’s system act to redistribute the Earth’s mass thereby changing the motion of the solid Earth relative to the centre, as well as causing
spatial and timedependent variations of the gravitational field of the Earth. Observations of this variability of the Earth’s gravity field using artificial satellites via the Satellite Laser
Ranging (SLR) technique and other geodetic techniques, Global Navigation Satellite Systems (GNSS), Very Long Baseline Interferometry and Doppler Orbitography and Radiolocation Integrated by Satellite
(DORIS) can be used to study a wide variety of geophysical processes that involve changes in mass.^{3} The Earth’s gravity field plays an important role in understanding dynamic processes taking place within the Earth system. These processes include interactions between the cryosphere, hydrosphere,
atmosphere and ocean at spatial scales ranging from a few metres to continental and global scales. Temporal scales of these dynamic processes range from an hour to geological time.^{4} The gravity
field of the Earth can also be used to determine global ocean circulation which relates to global climate change. Gravitational field changes may be used in detecting mass shifts in the Earth’s interior,
which might be associated with movements on the Earth’s surface.^{4} Characteristics of the gravitational field are often defined from artificial satellite tracking data. In particular, artificial
satellites are used to detect longwavelength components of the gravitational potential. A gravitational potential is often expressed as a series expansion of spherical harmonics known as a global geopotential
model (GGM). In this paper, a general review of GGMs, released from 1990 to date, is presented.
The concept of Satellite Laser Ranging


Gravity field models are often derived by use of data collected from the SLR observations technique. This is a technique that measures the twoway travel time of a short laser pulse which is reflected by an
orbiting satellite. SLR as an observational method is possible through satellites equipped with retroreflectors made from glass prisms. An example of retroreflectors is illustrated in Figure 1. In a typical
SLR system, a transmitting telescope emits short laser pulses with energy between 10 mJ and 100 mJ at a repetition frequency ranging between 5 Hz and 20 Hz. Some modern systems have lower
power levels and higher firing rates (of up to 2 kHz). Laser pulses which illuminate any of the retroreflectors are reflected back to the ground station where they are collected via the receiving telescope
and detected by a photomultiplier or a solid state photo diode. The measurement of laser ranges from laser tracking stations to a retroreflector on an orbiting satellite, for example, the timeofflight (TOF)
is often measured by either a time interval counter or an epoch timer.

FIGURE 1: The satellite LAGEOS 1 contains an array of retroreflective mirrors
covering its surface.


A basic equation representing the approximate TOF is given by: where c is the speed of light, d is the roundtrip distance from the SLR station to the target satellite retroreflector and t is the TOF. In order to obtain the best possible range precision
from the ground station to the satellite, additional parameters and numerous corrections corresponding to internal delays in the transmission and detection systems need to be considered. Taking into account such
corrections, the basic range equation given by [Eqn 1] can be expanded as in [Eqn 2] as reported in Seeber^{5}: In [Eqn 2], ∆t is the measured TOF and is mostly affected by uncertainties in the signal identification. The preferred resolution for the measured TOF is often a few picoseconds. In addition, the
measured TOF needs to be tied to universal time (because of the satellite’s motion relative to the Earth). The ∆d_{0} term corresponds to the eccentric correction on the ground,
which is the intersection of the vertical axis and horizontal axis and is used as a reference point in the laser system. Similarly, ∆d_{S} corresponds to the eccentric correction at the
satellite and gives a geometrical relationship between the centre of the corner cube and the centre of mass of the satellite; the accuracy of this parameter is very difficult to obtain on satellites with irregular
shapes (e.g. satellites equipped with solar panels and antennas). The ∆d_{b }term in [Eqn 2] corresponds to the signal delay in the ground system – the geometric reference point 0 to the
electrical 0 point and is often not exactly at the same point; this correctional parameter is often determined through calibration with older systems that were calibrated with respect to a defined terrestrial target.
Furthermore, ∆d_{r} is the refraction correction as a result of atmospheric conditions which affect the propagation velocity of laser pulses. Laser pulses experience a delay in the lower part of
the atmosphere, which makes measurement of these parameters along the total path difficult. Therefore atmospheric models are used that incorporate variables such as SLR site pressure and temperature and are supported
by measured data at the laser site.^{5} Lastly, η are random systematic and observation errors related to unmodelled residual effects. The first SLR experiment campaign began in the 1960s with the development of Rubybased SLR stations tracking satellites such as the Beacon ExplorerB. Since then numerous satellite missions have been launched for
different applications, such as geodetic, Earth sensing and radio navigation, and a global network of SLR stations has been established, replacing the old BakerNunn optical camera tracking network.^{6} A
historical overview of such missions is summarised in Table 1. The current global network of SLR stations involved in satellite tracking consists of over 40 stations and their global distribution is depicted in Figure 2.
Most of the stations are located in the Northern Hemisphere leaving the Southern Hemisphere with weak coverage. In Africa, there are two stations: Helwan in Egypt and MOBLAS6 (see Figure 3) located at Hartebeesthoek
Radio Astronomy Observatory (HartRAO) in South Africa. The space geodetic fundamental station HartRAO is involved with the International Laser Ranging Service activities as well as the other services of the International
Association of Geodesy. This SLR tracking station is relatively isolated in Africa and more active than Helwan; hence HartRAO plays a very important role as far as data coverage is concerned.
TABLE 1: Timeline of artificial satellites that were tracked by global Satellite
Laser Ranging stations.


FIGURE 2: The global distribution of the International Laser Ranging Service (ILRS) tracking network.



FIGURE 3: The Satellite Laser Ranging station MOBLAS6 is located at the
Hartebeesthoek Radio Astronomy Observatory in South Africa.


Forces acting on an orbiting satellite


Precise satellite tracking measurements provide orbit solutions which can be utilised to derive geopotential models. For instance, the longwavelength gravity field models can be derived through SLR range
measurements from highaltitude satellites such as LAser GEOdynamics Satellite (LAGEOS), Stella and Starlette. On the other hand, the mediumwavelength gravity field models are often computed by use of SLR
tracking data from low Earthorbiting satellites such as the Challenging Minisatellite Payload (CHAMP), the Gravity Recovery and Climate Experiment (GRACE) and the Gravity Field and Steadystate Ocean Circulation
Explorer (GOCE). Orbits of such satellites are altered by various gravitational, nongravitational and other unmodelled forces. The motion of a satellite in an inertial reference frame perturbed by gravitational,
nongravitational and unmodelled forces is often expressed by [Eqn 3], as reported in Seeber^{5}: where is the position vector of the centre of mass of the satellite, ā_{g} is the sum of the gravitational forces acting on the
satellite, ā_{ng}_{ }is the sum of the nongravitational forces acting on the surface of the satellite
and ā_{emp} represents the unmodelled forces which act on the satellite because of either a functionally incorrect or an incomplete
description of the various forces acting on the satellite.^{5} The gravitational forces (ā_{g}) acting on an orbiting satellite
are composed of a series of perturbations, expressed as: where_{ }Ρ̅ _{geo} is the geopotential force as a result of the gravitational attraction of the Earth; Ρ̅ _{set }and Ρ̅
_{ot} are perturbations as a result of solid Earth tides and ocean tides, respectively; Ρ̅ _{rd}_{ }is the perturbation caused by the rotational deformationof
the Earth; Ρ̅ _{ smp} is the perturbation caused by the Sun, Moon and other planets; and Ρ̅_{rel} is the perturbation as a result of general
relativity.^{5} The nongravitational forces acting on an orbiting satellite are given by: where Ρ̅_{drag } is the atmospheric drag acting on a satellite; Ρ̅_{solar} is the perturbation as a result of solar radiation pressure;
Ρ̅_{earth} is the perturbation caused by Earth radiation pressure; and Ρ̅_{thermal} is the perturbation as a result of thermal radiation
imbalances resulting from nonuniform temperature distribution on different satellite surfaces. The forces described in [Eqn 4] and [Eqn 5], together with unmodelled forces, are solved for during geopotential
modelling and therefore are central, in particular, to the derivation of precise geopotential models.
Precise satellite orbit determination


Precise satellite orbit determination (POD) is one of the most essential applications of geopotential modelling. The POD process involves the estimation of position and velocity of an orbiting satellite at a
specific time epoch.^{7} POD is used for geolocation of the satellite sensors and to measure the gravity field and its variations in time. There are currently three ways in which a satellite orbit can
be calculated: dynamic, kinematic and reduced dynamic.
Dynamic orbit determination
Dynamic orbit determination^{7} requires a set of tracking observations and mathematical models acting on an orbiting satellite. Here the force and satellite models are used to compute
a model of satellite acceleration over a given time. A nominal trajectory is generated analytically or numerically by integrating the acceleration model. The orbit solution is compared with the
one predicted by the observations. Selected parameters of the force models acting on the satellite may be adjusted along with an initial satellite position and velocity in order to minimise the
difference between the actual observations and the predicted ranges (this difference is called OC residuals) in a leastsquares sense.^{7} The accuracy of the dynamic orbit determination
approach is highly dependent on the satellite force models.
Kinematic orbit determination
Kinematic orbit determination is a purely geometric technique that depends only on GNSS (e.g. GPS) measurements and cannot be used by SLR.^{8} It does not take into account the dynamic properties
(e.g. gravity field or air drag) of an orbiting satellite. Here the errors emanating from the satellite force models do not affect the accuracy of the kinematic orbit determination, but its accuracy does
depend on the availability and accuracy of GNSS data.
Reduceddynamic orbit determination
In the dynamic and kinematic methods, the accuracy of the solution may be reduced as a result of modelling errors and GNSS measurement noise, respectively. The reduceddynamic technique proposed by
Yunck et al.^{9} may be defined as a method that exhibits a combination of dynamic and kinematic components and that minimises the errors caused by each method. In the reduceddynamic orbit
determination approach, the kinematic components of the dynamic force models are introduced in the form of a process noise model containing two parameters – the correlation time constant T (which
defines the correlation in the dynamic model error over one update interval) and the dynamic model steadystate variance V. Accuracy in this method depends on proper adjustment of the two parameters.
Global geopotential models


The perturbing potential function for the solidbody mass distribution of the Earth is often expressed in terms of a spherical harmonic expansion, obtained when solving a Laplace equation in spherical
coordinates described in Tapley et al.^{10} by: Here (r,φ,λ) represent the magnitude of the radius vector, the latitude and the longitude, respectively; {C̅ _{nm}, S_{nm}} are fully
normalised spherical harmonic coefficients of degree n and order m; and Ρ̅_{nm }is a fully normalised associated Legendre function. The adopted gravity mass
constant GM is set to 398600.4415 km^{3}/s^{2 }in most recent geopotential models.^{11} A typical geopotential model is often described by {C̅ _{nm},
S_{nm}} spherical harmonic coefficients. The values of {C̅ _{nm}, S_{nm}} coefficients decrease as the degree increases.
For satellitebased global gravity field models the accuracy of the lower degree coefficients
is typically higher than the higher degree coefficients. A number of spherical harmonic models have been developed over the years by different research groups, for example, Ohio State University (OSU), GeoForschungs Zentrum (GFZ) Potsdam, Goddard Earth Models (GEM),
Joint Gravity Models (JGM), Texas Earth Gravity models (TEG) and European Improved Gravity Model of the Earth by New Techniques (EIGEN). The existing GGMs representing the Earth’s gravitational field can be
classified into three groups: satelliteonly, combined and tailored gravity field models. The satelliteonly GGMs are derived from the analysis of the orbits of artificial Earth satellites. Numerous factors have
been attributed to the inherent inaccuracies of the satelliteonly models. These factors include weakening of the gravitational field with altitude, precession of the Earthbased range measurements to the satellites,
the lack of continuous tracking data from the existing stations and difficulties in modelling nongravitational and third body perturbations.^{12} The satelliteonly models are often combined with terrestrial
gravity data and marine gravity anomalies to yield highdegree combined GGMs. The combined GGMs are subjected to the same deficiencies as in satelliteonly GGMs and also to other errors emanating from terrestrial
gravity anomalies.^{13} In tailored GGMs, the spherical harmonic coefficients of the satelliteonly or the combined models are often adjusted and extended to higher degrees by using previously used or
unused gravity data.^{14} Currently a number of GGMs have been derived and made available freely to the scientific community.^{15} A review of gravity field models derived between 1970 and 1997 can be found in Rapp^{16}.
This paper reviews developments undergone in the gravity field modelling for the last two decades (i.e. 1990–2010). Characteristics of these models are summarised in Table 2. The first considered model is a
combined gravity field model, GRIM4C1, reported by Schwintzer et al^{17}. This model was computed as a joint collaboration between the German Geodetic Research Institute (DGFI) and Groupe de Recherches de
Geodesie Spatiale (GRGS). The GRIM4C1 model was derived up to degree and order 50 in terms of spherical harmonics. It incorporated GRIM4S1 satellite solution, mean gravity anomalies and Seasat altimeter derived mean
geoid undulations. The OSU91A geopotential model was reported by Rapp et al^{18}. This model was an upgraded version of OSU89a and OSU89b. It was computed complete to degree and order 360 in terms of spherical
harmonics in a blended form. In the computation of the OSU91A, coefficients to degree 50 were based on a combined solution from the GEMT2 model, surface gravity data and GEOSAT altimeter data. The remaining coefficients
(51–360) were derived from a combined solution computed from terrestrial data, altimeterderived anomalies and topographic or isostatic anomalies.
TABLE 2: Summary of some of the global geopotential models released between
1990 and 2008. Geophysical applications of these models include gravity field,
satellite orbit determination, station coordinates, reduction of altimeter data,
Earth’s rotation and computation of geoid undulations.

The Joint Gravity Model 3 (JGM3) released in 1994 was reported by Tapley et al.^{19} This model was developed by NASA Goddard Space Flight Center (GSFC) and the University of Texas at Austin as part of the
Topex Poseidon (T/P) project. This combined model was derived by adding the geopotential coefficients from the prelaunch model, JGM1 and their associated error covariance with GPS, SLR, DORIS tracking of T/P, laser
ranging tracking of LAGEOS 2 and Stella and DORIS tracking of SPOT 2. The model was derived complete to degree and order 70. The GRIM4S4 and GRIM4C4 reported by Schwintzer et al.^{20} were developed jointly
by GFZ Potsdam and GRGS Toulouse/Grasse for requirements of geodetic and altimeter satellite missions. The GRIM4S4 model was derived solely from satellite tracking data complete to degree and order 70. On the other
hand, the GRIM4C4 model was derived based on a least squares adjustment involving a combined solution from the GRIM4S4 model and surface gravity data from gravimetric and altimeter measurements. This model was computed
complete to degree and order 72, corresponding to a spatial resolution of 555 km at the surface of the Earth.^{20} The GRIM4S4 and GRIM4C4 models were thought to be efficient for satellite orbit computations
especially with orbit altitudes exceeding about 800 km.^{20} The GFZ96 geopotential model which was an upgrade of the GFZ93 and GFZ95 models was reported to provide high resolution in the history of
GFZderived models.^{21} This combined model was computed from the then improved terrestrial data derived from a 3year ERS1 mean sea surface and PMG055 solution. The solution was also combined with
altimeterderived gravity anomalies and normal equations and potential coefficients of the GRIM4S4 model as the a priori model. The GFZ96 model was derived to degree and order 359. Lemoine et al.^{22} described the combined spherical harmonic model, EGM96, which is complete to degree and order 360 and corresponds to a global resolution of about 55 km. The EGM96 model was developed
based on a joint collaboration between NASA–GSFC, the National Imagery and Mapping Agency (NIMA) and the Ohio State University. EGM96 is a blend model in which three computational procedures were used. The spherical
harmonic coefficients from 2–70 were derived based on a least squares adjustment involving satellite tracking data, terrestrial data and altimeter data of the ocean surface from the T/P, ERS1 and GEOSAT missions
and fillin gravity anomalies in areas lacking data.^{12} From degree 71–359, the coefficients were computed from a combined solution based on normal equations derived from the satellite tracking data which
were used as a priori values. The remaining coefficients at degree 360 were taken from a quadrature combined solution derived from the a priori satellite model and ERS1/GEOSAT altimeterderived anomalies. The EGM96
geopotential model was believed to provide a more accurate reference surface for the topography, as well as an improved orbit determination for low orbiting satellites.^{22} The GRIM5C1 gravity field model
reported by Gruber et al.^{23} was derived in a German–French joint collaboration between GFZ Potsdam and GRGS Toulouse. The model was computed up to degree and order 120. It incorporated terrestrial
and airborne mean gravity anomalies, altimetric gravity anomalies from NIMA and mean gravity anomalies derived from the GRIM5S1 model. Most of the geopotential models released from 2000 onwards were derived mostly from CHAMP, GRACE and GOCE missions plus other satellites, terrestrial and altimeter data. Geopotential models generated from the
inclusion of the three satellite missions data are believed to be more accurate than prior models (e.g. they allow, with an unprecedented accuracy and resolution, the recovery of the mean sea surface topography
from the difference between an altimetrybased mean sea surface height model and the gravity model’s derived geoid).^{24} The first CHAMP geopotential model, EIGEN1, reported by Reigber et al.^{25},
was derived in a German–French collaboration complete to degree and order 119. This model was derived by use of GPS tracking and 3 months of onboard accelerometer data from CHAMP. The EIGEN1 geopotential model was
reported to resolve the geoid and gravity with an accuracy of about 20 cm and 1 mGal, respectively, at a halfwavelength resolution of 550 km.^{25} The EIGEN2 model reported by Reigber
et al.^{26} was also derived in collaboration between Germany and France. This satelliteonly model was derived complete to degree and order 140. The model incorporated gravity orbit perturbations exploiting
GPS CHAMP satellitetosatellite tracking and 6 months of onboard accelerometer data. The accuracy in terms of geoid and gravity for the EIGEN2 model was reported to be about 10 cm and 0.5 mGal, respectively. Like the CHAMP mission, the GRACE mission data set has enabled a homogeneous determination of the geopotential gravity field modelling. The first model that was derived from the GRACE data was a ‘satelliteonly
model’ called the GGM01S reported by Tapley et al.^{10} This model, derived to complete degree and order 120, incorporated GRACE tracking data spanning from April to November 2002 (summing to a total of
111 selected days) and used least squares adjustment. The authors reported error estimates to an accuracy of about 2 cm over the land and ocean regions. An improvement on the geopotential model GGM01, called GGM02,
was released in 2005. GGM02 exists both in the GRACEbased satelliteonly GGM02S and the combined model GGM02C.^{27} The combined geopotential model incorporated the GRACEonly model GGM02S with EGM96 plus
14 months of GRACE data spanning from April 2002 to December 2003. The GGM02C model was computed to maximum degree and order 200 in terms of spherical harmonics. Improvements by a factor of two were reported with
error estimates of less than 1 cm geoid height to spherical harmonic at degree 70. The satelliteonly model EIGENGL04S1, described by Foerste et al.^{28}, has a maximum degree and order of 150. It incorporated GRACEonly (EIGEN–GRACE04S) and GRACE–LAGEOS (EIGEN–GL04S)
solutions. EIGEN–GL04S1 was later combined with surface gravity data from altimetry over the oceans and gravimetry over the continents to derive a highresolution gravity model, EIGEN–GL04C, released in
2006.^{28} This combined gravity field model is an outcome of the joint gravity field processing between GRGS Toulouse and GFZ Potsdam. The satellite part of EIGEN–GL04C is based on GRACE and LAGEOS data
and the maximum degree and order of this model is 360 in terms of spherical harmonics. EIGEN5C^{29} was also a joint collaboration between GFZ Potsdam and GRGS Toulouse. EIGEN5C is an upgrade of EIGEN–GL04C and has a maximum degree and order of 360. The model is again a combination of
GRACE and LAGEOS tracking data combined with gravimetry and altimeter surface data. The combination of the satellite and surface data has been done by combining normal equations obtained from observation equations
for the spherical harmonic coefficients. The National GeospatialIntelligence Agency released the first ever global model capable of resolving the Earth’s gravity field beyond a spherical harmonic degree of 2000;
this model is called EGM2008.^{30} A description of this model can be found in Pavlis et al.^{30} The EGM2008 gravity field model has a maximum degree and order of 2159. It incorporates improved gravity
anomaly data, altimetryderived gravity anomalies and GRACEbased satellite solutions. It allows proper computation of quasigeoid heights, gravity anomalies and vertical deflections and has a spatial resolution of ~5
arc minutes or ~9 km in the latitudinal direction.^{30}
Improvements in gravity field models


Improvement in gravity field modelling in terms of accuracy and spatial resolution is necessary in order to understand the physics of the interior of the Earth; the dynamics of the ocean and the interaction
of continents; to study the sea levels of ice and oceans; and to better determine satellite orbits and height systems in science and engineering.^{31} Such improvements are expected owing to the availability
of qualitative data, especially from the low Earthorbiting satellites. Satellite missions such as CHAMP, GRACE and GOCE (launched in 2000, 2002 and 2009, respectively), are believed to have improved the spatial
resolution sensitivity and accuracy of the newly developed GGMs.^{32} The CHAMP, GRACE and GOCE missions were designed to resolve the longwavelength part of the gravity field and hence provide unprecedented
accuracy.^{32} In contrast to the sporadic tracking by the SLR network of the ILRS, the three satellite missions carry GPS receivers on board that allow continuous orbit tracking. Furthermore, these satellites
are equipped with accelerometers which provide direct measurements of the nonconservative forces (e.g. air drag). In the case of GOCE, six accelerometers are installed in a gradiometer arrangement which additionally
allows for direct measurement of the Earth’s gravity gradients, which gives an improvement in the mediumwavelength part of the gravity.Improvement in gravity field modelling has already been noticed in several models as the resolution in such models is increased to reach higher degree and order in terms of spherical harmonics. Such improvements can
be measured by studying characteristics of the GGMs based on several factors. For example, the behaviour of GGMs can be analysed by performing orbit adjustment tests on artificial satellites and GPS or levelling tests,
and by comparing the spectral behaviour of the models or ocean geoid.^{33} Whilst old geopotential models derived up to degree and order 70 can resolve spatial features (geoid computation) at a half wavelength
of about 290 km, models (particularly the most recent) computed up to degree and order 360 can resolve spatial features down to 55 km.^{34} Early evaluations of gravity field models by Zhang and
Featherstone^{35} reported that the OSU91A geopotential model provided better fits to the gravity field over the Australian region than did prior models. In contributions by Pearse and KearsIey^{36}
and Kirby et al.^{37}, the OSU91A model has been ousted by the EGM96 model, which reportedly gives better solutions to the computation of geoid heights. Evaluations of GGMs released between 1996 and 2002 by Amos and Featherstone^{12}, based on comparisons of gravity anomalies, freeair gravity anomalies, geoid heights and GPS or levelling tests, found that
EIGEN1S was the best satelliteonly GGM when applied in the Australian and New Zealand region, whilst the best combined GGM over the same region was reported to be PGM2000A. The quality of the GGM01 model was
assessed by Ellmann^{38} by comparing it with the combined gravity field model EGM96. It was reported that the GGM01 model gives better solutions of gravity anomaly and geoidal heights over Fennoscandia
(e.g. Finland, Germany, Norway and Sweden) and the Baltic Sea region. In a comparison study of 10 geopotential models (EGM96, GGM02C, GGM03S, ITGGRACE03, JEM01RL03B, EIGENGL04C, EIGEN5C/5S and EGM2008) using
geoid heights and GPS or levelling data points, the EGM2008 model was reported to provide the best solution compared to the other models at degree 360.^{33} A much improved solution was also reported for
EGM2008 when its coefficients were increased to degree 2190. A similar study evaluating the GGMs EGM96, EIGEN5C and EGM2008, based on the comparison of geoid heights to the GPS or levelling over Afyonkarahisar
in Western Turkey, has also confirmed the improvements in the EGM2008 model in the computation of geoid heights.^{39} Botai and Combrinck^{40} investigated the general improvement of 13 gravity field
models released between 1990 and 2008 based on LAGEOS data and they found that gravity field modelling had improved during this period by a factor of two (Figure 4). In this review, the GRACEonly gravity field model,
AIUBGRACE01S, has been shown to provide lower root mean square orbital errors compared to the other tested models.

FIGURE 4: Time series of the mean root mean square (RMS) values for 13 evaluated global geopotential models based on the LAGEOS satellite.


Geophysical applications of gravity field modelling


According to Newton’s law, changes in the gravity field are evidence of changes in mass and density distribution in the Earth system. Any movement of masses in, on or above the Earth will therefore
introduce variations in the gravity field of the Earth.^{2,41} Temporal variations in Earth’s gravity field are often in the order of 10^{6 }N/kg (variation from the mean) and occur on
a scale ranging from hours to thousands of years.^{42} Such temporal variations are caused by several phenomena that redistribute mass, which include tides caused by the Sun and Moon, and postglacial
rebound caused by isostatic correction. Surface mass changes in the atmosphere, oceans, hydrosphere and cryosphere are dominated by seasonal and interannual variations, whilst processes such as isostatic
glacial recovery and sealevel change give rise to longterm secular or quasisecular signatures.^{43} Several studies have investigated the longterm and the seasonal variations of the Earth’s gravity field using data collected from different satellite missions.^{44,45,46} In particular, the
lower order harmonic component of the gravity field with n = 2 and m = 0 (hereafter J_{2}), which characterises the gravitational oblateness of the Earth has attracted a lot of interest
from the scientific community. Early studies of J_{2}, for example, that by Yoder et al.^{47}, showed a secular increase in J_{2} that was consistent with a steady migration of mass from
low latitudes towards high latitudes, resulting in a linearly decreasing trend. Such a trend was thought to be related to postglacial rebound, the Earth’s ongoing response to the removal of the ice loads at
the end of the last Ice Age. However, longterm studies by Cox and Chao^{48} have discovered that J_{2} started to increase from about 1997. The detected increasing trend reported in Cox and Chao^{48} is believed to have turned again, with J_{2} once more decreasing from late 1997. Several mechanisms have been suggested to be the cause
in the sudden changes of the J_{2} coefficient. These mechanisms include processes involved in the surge in subpolar glacial melting and mass shifts in the Southern, Pacific and Indian Oceans.^{1}
In addition to the increasing trend of the J_{2} coefficient, Nerem et al.^{49} reported that the J_{2} coefficient might be exhibiting seasonal variability as a result of a combination of
atmospheric pressure variations and variations in the distribution of water in the oceans and on land. Furthermore, Dickey et al.^{2} detected interannual variability in the J_{2} which they
attributed to climatically driven oscillations in the ocean, storage of water, snow and ice on land, and also partly to the consequence of the effects of anelasticity on the 18.6year solid Earth tide, as suggested
by Benjamin et al.^{50} Earth orientation changes, often represented by polar motion, Xequatorial and Yequatorial polar components in a geographical reference frame and variations in the length of day (LOD), are often explained by studying
variations of atmospheric and/or oceanic angular momentum. Such variations are caused by the exchange of angular momentum between the solid Earth and its geophysical fluid envelope. Eubanks^{51} found that
variations in the Earth’s rotation rate corresponded to changes in LOD and amounted to a few parts in 10^{8}. Studies by Ponsar et al.^{52} suggested that variations in LOD are caused by the
interaction between the Earth’s core and mantle. Similar studies by Gross et al.^{53} related the LOD variations with tidal variations, exhibiting periods between 12 h and 18.6 years. Such variations
were believed to be as a result of the deformation of solid Earth and changes in the strength and direction of the winds. Geopotential models can be used to derive the geoid, the equipotential surface of the Earth’s gravity field that corresponds closely with mean sea level in the open oceans, ignoring oceanographic effects as well
as the geoidal height which is the separation between the geoid and the ellipsoid.^{54} Determination of the geoid has been one of the main research areas in Geodesy for decades, because of its various
applications, which include vertical data for orthometric heights, understanding of ocean circulation patterns and dynamics, refinement of satellite orbits and the modelling of geodynamical phenomena. To this end,
geoid heights at any point on the Earth’s surface can be determined with an accuracy ranging from 30 cm to a few metres.^{55} A number of researchers have addressed the precise determination of
geoid height on local and regional scales for oceanographic and geophysical applications.^{56,57} At a local scale, the geoid can be determined by a combination of GPSderived heights and levelled heights,
through gravimetric and geometric approaches. For instance, the quasigeoid for southern Africa based on SLRderived geopotential gravity models has been reported by Merry^{58}. In Merry’s^{58}
study, gravity data for South Africa was combined with different geopotential models (based on the removerestore procedure) to derive a quasigeoid model, UCT2006. In addition, quasigeoids produced in South Africa
by use of different GGMs were compared with GPS or levelling data points to assess the suitability and reliability of the considered models. Merry^{58} concluded that the UCT2006 model gives a better solution
(with a root mean square fit of 15 cm) compared to the EGM96, EIGENCG03C and GGM02C models. A slight improvement of 4 cm was also reported when the UCT2006 geoid model was allowed to tilt in two directions
(north–south and east–west). Global gravity change has also attracted particular attention in the scientific community as it is often related to global sealevel changes.^{59,60} The sources of global sealevel rise often involve the
redistribution of mass from the continents to the ocean. The use of gravity field measurements allows the discrimination of several sources through the continual monitoring of geoid changes on both global and regional
scales as well as on basin scales. Gravity field solutions can be used to numerically estimate components such as thermal expansion (eustatic) and freshwater influx influencing global sea level. Measurements of temporal
gravity variations can also be used to determine water storage change in the hydrological system. Since the launch of the GRACE mission in 2002, numerous articles assessing the potential of GRACE recovering hydrological
signals have been published. For example, Andersen and Hinderer^{61} have investigated the potential of inferring interannual gravity field changes caused by continental water storage change, as determined from
GRACE observations between 2002 and 2003. Contributions from continental water storage change were compared to the output from global hydrological models. Andersen et al.^{62} and Neumeyer et al.^{63}
correlated large scale hydrological events with the estimated change in the gravity field for certain areas of the world to an accuracy of 0.4 µGa corresponding to 9 mm of water. On a regional scale,
Winsemius et al.^{64} compared hydrological model outputs for the Zambezi river basin with estimates derived from GRACE. Monthly storage depths produced by the hydrological model displayed larger amplitudes
and were partly out of phase compared to the estimates based on GRACE data.
The continuous design and deployment of satellite missions dedicated to gravity field measurements and the availability of highprecision data have led to the availability of gravity information with
unprecedented spatialtemporal resolution and accuracy. In particular, the advent of satellite data has made it possible to determine the gravity field of the Earth via modelling. To this end, these data
sets are the basis of robust gravity field modelling with more than 100 gravity models being released to the scientific community since the early 1960s. Research dedicated to assessing the accuracy of the
existing gravity field models has been reported in the literature. Different gravity field models could be characterised by various degrees of spatialtemporal resolution. Such improvements are as a result
of the availability of quantitative and qualitative SLR and terrestrial gravity data. With the development of gravity field models, dedicated review papers that report on the chronology of gravity field
modelling for geophysical applications have been lacking. Review papers known thus far are more than a decade old and therefore cannot provide a complete account of uptodate gravity field modelling activities.
This review has explored the various gravity field modelling efforts with specific geophysical applications. It is concluded that gravity field modelling algorithms have improved over time partly due to the
availability of specialised gravity mission satellite data and partly because of the advancement of technology and Earth and orbit system modelling techniques or approaches.
The authors wish to thank the anonymous reviewers for their constructive and detailed comments that helped to improve the quality of the manuscript.
Competing interests
The authors declare that they have no financial or personal relationship(s) which may have inappropriately influenced them in writing this article.
Authors’ contributions
C.M.B. conceptualised the review, prepared the plotting scripts and wrote the manuscript. L.C. made conceptual contributions and edited the manuscript.
1. Blakely RJ. Potential theory in gravity and magnetic applications. Cambridge, NY: Cambridge University Press; 1995.
http://dx.doi.org/10.1017/CBO97805115498162. Dickey JO, Marcus SL, De Viron O, Fukumori I. Recent Earth oblateness variations: Unraveling climate and postglacial rebound effects. Science. 2002;298:1975–1977.
http://dx.doi.org/10.1126/science.1077777,
PMid:12471253
3. Dickey JO, Bentley CR, Bilham R, et al. Satellite gravity and the geosphere: Contributions to the study of the solid Earth and its fluid envelope. Washington DC: National Academic Press; 1997. 4. Forsberg R, Sideris MG, Shum CK. The gravity field and GGOS. J Geodynamics. 2005;40:387–393.
http://dx.doi.org/10.1016/j.jog.2005.06.014 5. Seeber G. Satellite geodesy: Foundations, methods, and applications. 2nd ed. Berlin: Walter de Gruyter; 2003.
http://dx.doi.org/10.1515/9783110200089 6. Combrinck L. Satellite Laser Ranging. In: Xu G, editor. Sciences of geodesy I – Advances and future directions. Berlin: Springer Verlag, 2010; p. 301–336. 7. Yunck TP. Orbit determination. In: Parkinson BW, Spilker JJ, editors. Global positioning system: Theory and applications, Vol 2. Progress in Astronautics and Aeronautics. Reston, VA: American
Institute of Aeronautics and Astronautics, 1997; p. 567–585. 8. Colombo OL, Luthcke SB. Kinematic point positioning of a LEO, with simultaneous reduceddynamic orbit estimation. Paper presented at: ION GNSS 2004. Proceedings of the 17th International Technical
Meeting of the Satellite Division of the Institute of Navigation; 2004 Sept 21–24; Long Beach, CA, USA. 9. Yunck TP, Bertiger WI, Wu SC, et al. First assessment of GPSbased reduced dynamic orbit determination on TOPEX/Poseidon. Geophys Res Lett. 1993;21:2179–2182. 10. Tapley BD, Bettadpur S, Watkins M, Reigber C. The gravity recovery and climate experiment: Mission overview and early results. Geophys Res Lett. 2004;31:L09607.
http://dx.doi.org/10.1029/2004GL019920 11. Ellmann A, Jurgenson H. Evaluation of a GRACEbased combined geopotential model over the Baltic countries. Geod Cartograf. 2008;34(2):35–44. 12. Amos MJ, Featherstone WE. Comparisons of global geopotential models with terrestrial gravity field data over New Zealand and Australia. Geomatics Res Australas. 2003;78:67–84. 13. Heck B. An evaluation of some systematic error sources affecting terrestrial gravity anomalies. Bull Geodesique. 1990;64:88–108.
http://dx.doi.org/10.1007/BF02530617 14. Kearsley AHW, Forsberg R. Tailored geopotential models – Applications and shortcomings. Manuscripta geodaetica. 1990;15:151–158. 15. International Centre for Global Earth Models [homepage on the Internet]. c2011 [cited 2011 Aug 04]. Available from: http://icgem.gfzpotsdam.de/ICGEM/ICGEM.html 16. Rapp RH. Past and future developments in geopotential modelling. In: Forsberg R, Feissl M, Dietrich R, editors. Geodesy on the move. Berlin: Springer, 1997; p. 58–78. 17. Schwintzer P, Reigber C, Massmann FH, et al. A new Earth gravity field model in support of ERS1 and SPOT2, Final report to the German Space Agency (DARA) and the French Space Agency (CNES). Munich: German
Geodetic Research Institute; Toulouse: Groupe de Recherches de Geodesie Spatiale; 1991. 18. Rapp RH, Wang YM, Pavlis NK. The Ohio State 1991 geopotential and sea surface topography harmonic coefficient models. Report 410. Columbus, OH: Department of Geodetic Science and Surveying, Ohio State
University; 1991. 19. Tapley B, Watkins M, Ries J, et al. The joint gravity model 3. J Geophys Res. 1996;101(B12):28029–28049.
http://dx.doi.org/10.1029/96JB01645 20. Schwintzer P, Reigber C, Bode A, et al. Longwavelength global gravity field models: GRIM4S4, GRIM4C4. J Geodesy. 1997;71:189–208.
http://dx.doi.org/10.1007/s001900050087 21. Gruber T, Anzenhofer M, Rentsch M. Improvements in high resolution gravity field modeling at GFZ. In: Segawa J, Fujimoto H, Okubo S, editors. Gravity, geoid and marine geodesy. Berlin:
Springer, 1997; p. 445–452. 22. Lemoine FG, Kenyon SC, Factor JK, et al. The development of the joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96. NASA Technical Paper NASA/TP1998206861.
Greenbelt, MD: Goddard Space Flight Center; 1998. 23. Gruber T, Bode A, Reigber C, et al. GRIM5C1: Combination solution of the global gravity field to degree and order 120. Geophys Res Lett. 2000;27(24):4005–4008.
http://dx.doi.org/10.1029/2000GL011589 24. Dobslaw H, Schwintzer P, Barthelmes F, et al. Geostrophic ocean surface velocities from TOPEX altimetry, and CHAMP and GRACE satellite gravity models. Scientific technical report 04/07. Potsdam:
GeoForschungs Zentrum; 2004. 25. Reigber C, Balmino G, Schwintzer P. A high quality global gravity field model from CHAMP GPS tracking data and accelerometry (EIGEN1S). Geophys Res Lett. 2002;29(14):1–4.
http://dx.doi.org/10.1029/2002GL015064 26. Reigber C, Schwintzer P, Neumayer KH, et al. The CHAMPonly Earth gravity field model EIGEN2. Adv Space Res. 2003;31(8):1883–1888.
http://dx.doi.org/10.1016/S02731177(03)001625 27. Tapley B, Ries J, Bettadpur S, et al. GGM02 – An improved Earth gravity field model from GRACE. J Geodesy. 2005;79:467–478.
http://dx.doi.org/10.1007/s001900050480z 28. Foerste C, Schmidt R, Stubenvoll R, et al. The GFZ/GRGS satellite and combined gravity field models EIGENGL04S1 and EIGENGL04C. J Geodesy. 2006;82(6):331–346. 29. Foerste C, Flechtner F, Schmidt R. A new global combined highresolution GRACEbased gravity field model of the GFZGRGS cooperation. Geophys Res Abstr. 2008;10:EGU2008A03426. 30.Pavlis NK, Holmes SA, Kenyon SC, Factor JK. An Earth gravitational model to degree 2160: EGM2008. Vienna: EGU General Assembly; 2008. 31. Rummel R, Balmino G, Johnhannessen J, Visser P, Woodworth P. Dedicated gravity field missions – Principles and aims. J Geodyn. 2002;33:3–20.
http://dx.doi.org/10.1016/S02643707(01)000503 32. Featherstone WE. Improvement to longwavelength Australian gravity anomalies expected from the CHAMP, GRACE and GOCE dedicated satellite gravimetry missions. Explor Geophys (Collingwood, Aust). 2003;34:69–76.
http://dx.doi.org/10.1071/EG03069 33. Foerste C, Stubenvoll R, Koenig R, et al. Evaluation of EGM2008 by comparison with other recent global gravity field models. Newton’s Bull. 2009;4:26–37. 34. Moore P, Zhang Q, Alothman A. Recent results on the modeling of the spatial and temporal structure of the Earth’s gravity field. Phil Trans R Soc. 2006;A364:1009–1026.
http://dx.doi.org/10.1098/rsta.2006.1751 35. Zhang KF, Featherstone WE. The statistical fit of high degree geopotential models to the gravity field of Australia. Geomatics Res Australas. 1995;63:1–18. 36. Pearse MB, Kearsley AHW. Analysis of EGM models in New Zealand. Int Geoid Serv Bull. 1996;6:203–213. 37. Kirby JF, Featherstone WE, Kearsley AHW. Tests of the DMA/GSFC geopotential models over Australia. Int Geoid Serv Bull. 1998;7:213. 38. Ellmann A. The geoid for the Baltic countries determined by the least squares modification of Stokes’ formula. Geodesy Report No. 1061. Stockholm: Royal Institute of Technology,
Department of Infrastructure; 2004. 39. Yilmaz I, Yilmaz M, Gullu M, Turgut B. Evaluation of recent global geopotential models based on GPS/leveling data over Afyonkarahisar (Turkey). Sci Res Essays. 2010;5(5):484–493. 40. Botai MC, Combrinck L. Investigating the accuracy of gravity field models using Satellite Laser Ranging data. S Afr J Geol. 2011;114(3–4):539–544. 41. Cox CM, Au A, Boy JP, Chao BF. Timevariable gravity: Using SatelliteLaserRanging as a tool for observing longterm changes in the Earth system. In: Noomen R, Klosko S, Noll C, Pearlman M, editors.
Proceedings of the 13th International Workshop on Laser Ranging. Greenbelt, MD: NASA Goddard Space Flight Center; 2003. 42. Tapley BD, Bettadpur S, Ries JC, Thompson PF, Watkins MM. GRACE measurements of mass variability in the Earth system. Science. 2004;305:503–505. 43. Wahr J, Molenaar M, Bryan F. Time variability of the Earth’s gravity field: Hydrological and oceanic effects and their possible detection using GRACE. J Geophys Res. 1998;103:30205–30230.
http://dx.doi.org/10.1029/98JB02844 44. Klosko S, Chao B. Secular variations of the zonal gravity field, global sea level and polar motion as geophysical constraints. Phys Chem Earth. 1998;23(9–10):1091–1102.
http://dx.doi.org/10.1016/S00791946(98)001499 45. Chen JL, Wilson CR. Low degree gravitational changes from Earth rotation and geophysical models. Geophys Res Lett. 2003;30(24):2257.
http://dx.doi.org/10.1029/2003GL018688 46. Chen JL, Wilson CR, Tapley BD. Interannual variability of lowdegree gravitational change, 1980–2002. J Geodesy. 2005;78:535–543.
http://dx.doi.org/10.1007/s001900040417y 47. Yoder CF, Williams JG, Dickey JO, Schutz BE, Eanes RJ, Tapley BD. Secular variations of Earth’s gravitational harmonic J_{2} coefficient from Lageos and the nontidal
acceleration of Earth rotation. Nature. 1983;303:757–762. http://dx.doi.org/10.1038/303757a0 48. Cox CM, Chao BF. Detection of a largescale mass redistribution in the terrestrial system since 1998. Science. 2002;297:831–832.
http://dx.doi.org/10.1126/science.1072188,
PMid:12161652
49. Nerem RS, Eanes RJ, Thompson PF, Chen JL. Observations of annual variations of the Earth’s gravitational field using satellite laser ranging and geophysical models Geophys Res Lett. 2000;27:1783–1786.
http://dx.doi.org/10.1029/1999GL008440 50. Benjamin D, Wahr J, Ray RD, Egbert GD, Desai SD. Constraints on mantle anelasticity from geodetic observations, and implications for the anomaly. Geophys J Int. 2006;165:3–16.
http://dx.doi.org/10.1111/j.1365246X.2006.02915.x 51. Eubanks TM. Variations in the orientation of the Earth, in contributions of space geodesy to geodynamics: Earth dynamics. Geodyn Series volume 24. Washington DC: American Geophysical Union, 1993; p. 1–54. 52. Ponsar S, Dehant V, Holme R., Jault D, Pais A, Van Hoolst T. The core and fluctuations in the Earth’s rotation, in Earth’s core: Dynamics, structure, rotation. Washington DC: Geodyn,
American Geophysical Union, 2003; p. 251–261. 53. Gross RS, Fukumori I, Menemenlis D. Atmospheric and oceanic excitation of the Earth’s wobbles during 1980–2000. J Geophys Res. 2003;109:B01406.
http://dx.doi.org/10.1029/2003JB002432 54. Eckman M. What is the geoid? In: Vermeer M, editor. Coordinate systems, GPS, and the geoid. Report 95:5 of the Finnish Geodetic Institute. Masala: Finnish Geodetic Institute, 1998; p.
49–51. 55. Rapp RH. Global models for the 1 cm geoid – Present status and near term prospects. In: Sanso F, Rummel R, editors. Geodetic boundary value problems in view of the one centimeter
geoid. Berlin: Springer, 1997; p. 273–311. 56. Kiamehr R, Sjoberg LE. Comparison of the qualities of recent global and local gravimetric geoid models in Iran. Stud Geophys Geod. 2005;49:289–304.
http://dx.doi.org/10.1007/s1120000500117 57. Chandler G, Merry CL. The South African geoid 2010: SAGEOID10. Position IT. 2010(June):29–33. 58. Merry C. Evaluation of global geopotential models in determining the quasigeoid for South Africa. Surv Rev. 2007;39(305):180–192. 59. Cazenave A, Nerem RS. Presentday sea level change: Observations and causes. Rev Geophys. 2004;42:RG3001.
http://dx.doi.org/10.1029/2003RG000139
60. Woodworth PL. Some important issues to do with longterm sea level change. Phil Trans R Soc A. 2006;364:787–803.
http://dx.doi.org/10.1098/rsta.2006.1737,
PMid:16537140
61. Andersen OB, Hinderer J. Global interannual gravity changes from GRACE: Early results. Geophys Res Lett. 2005;32:L01402.
http://dx.doi.org/10.1029/2004GL020948 62. Andersen OB, Seneviratne SI, Hinderer J, Viterbo P. GRACEderived terrestrial water storage depletion associated with the 2003 European heat wave. Geophys Res Lett. 2005;32:L18405.
http://dx.doi.org/10.29/2005GL023574 63. Neumeyer J, Barthelmes F, Dierks O, et al. Combination of temporal gravity variations resulting from superconducting gravimeter (SG) recordings, GRACE satellite observations and global hydrology
models. J Geodesy. 2006;79(10–11):573–585. http://dx.doi.org/10.1007/s0019000500148 64. Winsemius HC, Savenije HHG, Van den Hurk NC, Zapreeva EA, Klees R. Assessment of gravity recovery and climate experiment (GRACE) temporal signature over the upper Zambezi. Water Resour Res. 2006;42:2–8. 65. Torrence M. Lageos1, 2 [homepage on the Internet]. No date [cited 2012 Feb 08]. Available from:
http://ilrs.gsfc.nasa.gov/satellite_missions/list_of_satellites/lag1_general.html 66. International Laser Ranging Service [homepage on the Internet]. c2009 [updated 2009 Feb 26, cited 2012 Feb 08]. Available from:
http://ilrs.gsfc.nasa.gov/images/ilrsmap.jpg 67. Combrinck L, Haupt W. MOBLAS6 at HartRAO [homepage on the Internet]. No date [cited 2012 Feb 08]. Available from:
http://www.hartrao.ac.za/geodesy/NASA_Network_station_report.Mob6.htm
