Forces acting on an orbiting satellite Precise satellite tracking measurements provide orbit solutions which can be utilised to derive geopotential models. For instance, the long-wavelength gravity field models can be derived through SLR range measurements from high-altitude satellites such as LAser GEOdynamics Satellite (LAGEOS), Stella and Starlette. On the other hand, the medium-wavelength gravity field models are often computed by use of SLR tracking data from low Earth-orbiting satellites such as the Challenging Minisatellite Payload (CHAMP), the Gravity Recovery and Climate Experiment (GRACE) and the Gravity Field and Steady-state Ocean Circulation Explorer (GOCE). Orbits of such satellites are altered by various gravitational, non-gravitational and other unmodelled forces. The motion of a satellite in an inertial reference frame perturbed by gravitational, non-gravitational and unmodelled forces is often expressed by [Eqn 3], as reported in Seeber 5: http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/662/3390 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 non-gravitational 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: http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/662/3391 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 non-gravitational forces acting on an orbiting satellite are given http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/662/3392 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 non-uniform 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. 7POD 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 7requires 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 O-C residuals) in a least-squares sense. 7The 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. 8It 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. Reduced-dynamic 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 reduced-dynamic technique proposed by Yunck et al. 9may 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 reduced-dynamic 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 steady-state variance V. Accuracy in this method depends on proper adjustment of the two parameters. Global geopotential models The perturbing potential function for the solid-body 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. 10by: http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/662/3393 Here (r,φ,λ) represent the magnitude of the radius vector, the latitude and the longitude, respectively; {C̅ nm, Snm} 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 km3/s2 in most recent geopotential models.11 A typical geopotential model is often described by {C̅ nm, Snm} spherical harmonic coefficients. The values of {C̅ nm, Snm} coefficients decrease as the degree increases. For satellite-based 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: satellite-only, combined and tailored gravity field models. The satellite-only GGMs are derived from the analysis of the orbits of artificial Earth satellites. Numerous factors have been attributed to the inherent inaccuracies of the satellite-only models. These factors include weakening of the gravitational field with altitude, precession of the Earth-based range measurements to the satellites, the lack of continuous tracking data from the existing stations and difficulties in modelling non-gravitational and third body perturbations. 12The satellite-only models are often combined with terrestrial gravity data and marine gravity anomalies to yield high-degree combined GGMs. The combined GGMs are subjected to the same deficiencies as in satellite-only GGMs and also to other errors emanating from terrestrial gravity anomalies. 13In tailored GGMs, the spherical harmonic coefficients of the satellite-only or the combined models are often adjusted and extended to higher degrees by using previously used or unused gravity data. 14Currently a number of GGMs have been derived and made available freely to the scientific community. 15A 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 al17. 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 GEM-T2 model, surface gravity data and GEOSAT altimeter data. The remaining coefficients (51–360) were derived from a combined solution computed from terrestrial data, altimeter-derived anomalies and topographic or isostatic anomalies. 2Summary 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, Earths rotation and computation of geoid undulations. http://sajs.co.za/index.php/SAJS/article/downloadSuppFile/662/3399

The Joint Gravity Model 3 (JGM3) released in 1994 was reported by Tapley et al. 19This 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. 20were 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. 20The GRIM4S4 and GRIM4C4 models were thought to be efficient for satellite orbit computations especially with orbit altitudes exceeding about 800 km. 20The GFZ96 geopotential model which was an upgrade of the GFZ93 and GFZ95 models was reported to provide high resolution in the history of GFZ-derived models. 21This combined model was computed from the then improved terrestrial data derived from a 3-year ERS-1 mean sea surface and PMG055 solution. The solution was also combined with altimeter-derived 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. 22described 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, ERS-1 and GEOSAT missions and fill-in gravity anomalies in areas lacking data. 12From 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 ERS-1/GEOSAT altimeter-derived 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. 22The GRIM5C1 gravity field model reported by Gruber et al. 23was 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 altimetry-based mean sea surface height model and the gravity model’s derived geoid). 24The first CHAMP geopotential model, EIGEN-1, 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 on-board accelerometer data from CHAMP. The EIGEN-1 geopotential model was reported to resolve the geoid and gravity with an accuracy of about 20 cm and 1 mGal, respectively, at a half-wavelength resolution of 550 km. 25The EIGEN-2 model reported by Reigber et al. 26was also derived in collaboration between Germany and France. This satellite-only model was derived complete to degree and order 140. The model incorporated gravity orbit perturbations exploiting GPS CHAMP satellite-to-satellite tracking and 6 months of on-board accelerometer data. The accuracy in terms of geoid and gravity for the EIGEN-2 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 ‘satellite-only model’ called the GGM01S reported by Tapley et al. 10This 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 GRACE-based satellite-only GGM02S and the combined model GGM02C. 27The combined geopotential model incorporated the GRACE-only 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 satellite-only model EIGEN-GL04S1, described by Foerste et al. 28, has a maximum degree and order of 150. It incorporated GRACE-only (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 high-resolution gravity model, EIGEN–GL04C, released in 2006. 28This 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. EIGEN-5C 29was also a joint collaboration between GFZ Potsdam and GRGS Toulouse. EIGEN-5C 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 Geospatial-Intelligence 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. 30The EGM2008 gravity field model has a maximum degree and order of 2159. It incorporates improved gravity anomaly data, altimetry-derived gravity anomalies and GRACE-based satellite solutions. It allows proper computation of quasi-geoid heights, gravity anomalies and vertical deflections and has a spatial resolution of ~5 arc minutes or ~9 km in the latitudinal direction. 30Improvements 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. 31Such improvements are expected owing to the availability of qualitative data, especially from the low Earth-orbiting 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. 32The CHAMP, GRACE and GOCE missions were designed to resolve the long-wavelength part of the gravity field and hence provide unprecedented accuracy. 32In 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 non-conservative 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 medium-wavelength 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. 33Whilst 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. 34Early evaluations of gravity field models by Zhang and Featherstone 35reported 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 36and 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, free-air gravity anomalies, geoid heights and GPS or levelling tests, found that EIGEN-1S was the best satellite-only 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 38by 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, ITG-GRACE03, JEM01-RL03B, EIGEN-GL04C, EIGEN-5C/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. 33A much improved solution was also reported for EGM2008 when its coefficients were increased to degree 2190. A similar study evaluating the GGMs EGM96, EIGEN-5C 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. 39Botai and Combrinck 40investigated 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 GRACE-only gravity field model, AIUB-GRACE01S, has been shown to provide lower root mean square orbital errors compared to the other tested models. 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 sea-level change give rise to long-term secular or quasi-secular signatures. 43Several studies have investigated the long-term 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 J2), which characterises the gravitational oblateness of the Earth has attracted a lot of interest from the scientific community. Early studies of J2, for example, that by Yoder et al.47, showed a secular increase in J2 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, long-term studies by Cox and Chao 48have discovered that J2 started to increase from about 1997. The detected increasing trend reported in Cox and Chao 48is believed to have turned again, with J2 once more decreasing from late 1997. Several mechanisms have been suggested to be the cause in the sudden changes of the J2 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 J2 coefficient, Nerem et al. 49reported that the J2 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 J2 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.6-year solid Earth tide, as suggested by Benjamin et al. 50 Earth orientation changes, often represented by polar motion, X-equatorial and Y-equatorial 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 51found 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. 52suggested that variations in LOD are caused by the interaction between the Earth’s core and mantle. Similar studies by Gross et al. 53related 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. 54Determination 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. 55A 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 GPS-derived heights and levelled heights, through gravimetric and geometric approaches. For instance, the quasi-geoid for southern Africa based on SLR-derived geopotential gravity models has been reported by Merry 58. In Merry’s 58study, gravity data for South Africa was combined with different geopotential models (based on the remove-restore procedure) to derive a quasi-geoid model, UCT2006. In addition, quasi-geoids 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 58concluded that the UCT2006 model gives a better solution (with a root mean square fit of 15 cm) compared to the EGM96, EIGEN-CG03C 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 sea-level changes. 59,60 The sources of global sea-level 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 61have 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. 62and Neumeyer et al. 63correlated 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. 64compared 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. Summary The continuous design and deployment of satellite missions dedicated to gravity field measurements and the availability of high-precision data have led to the availability of gravity information with unprecedented spatial-temporal 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 spatial-temporal 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 up-to-date 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. Acknowledgements 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.http://dx.doi.org/10.1017/CBO97805115498162.http://dx.doi.org/10.1126/science.1077777124712533.DickeyJOBentleyCRBilhamRet al.Satellite gravity and the geosphere: Contributions to the study of the solid Earth and its fluid envelope4.http://dx.doi.org/10.1016/j.jog.2005.06.0145.http://dx.doi.org/10.1515/97831102000896.CombrinckL Satellite Laser Ranging7.YunckTP Orbit determinationParkinson BW, Spilker JJ Global positioning system: Theory and applications8.ColomboOLLuthckeSBKinematic point positioning of a LEO, with simultaneous reduced-dynamic orbit estimation9.YunckTPBertigerWIWuSCet alFirst assessment of GPS-based reduced dynamic orbit determination on TOPEX/Poseidon Geophys Res Lett.1993212179218210. http://dx.doi.org/10.1029/2004GL01992011.EllmannAJurgensonHEvaluation of a GRACE-based combined geopotential model over the Baltic countriesGeod Cartograf2008342354412.AmosMJFeatherstoneWEComparisons of global geopotential models with terrestrial gravity field data over New Zealand and AustraliaGeomatics Res Australas200378678413.http://dx.doi.org/10.1007/BF0253061714.KearsleyAHWForsbergRTailored geopotential models – Applications and shortcomingsManuscripta geodaetica19901515115815. International Centre for Global Earth Models [homepage on the Internet]. c2011 [cited 2011 Aug 04]. Available from: http://icgem.gfz-potsdam.de/ICGEM/ICGEM.html 16. RappRHPast and future developments in geopotential modellingForsberg R, Feissl M, Dietrich RGeodesy on the move. Berlin: Springer1997587817. SchwintzerPReigberCMassmannF-Het alA new Earth gravity field model in support of ERS-1 and SPOT-2, 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 Spatiale199118. RappRHWangYMPavlisNKThe Ohio State 1991 geopotential and sea surface topography harmonic coefficient modelsReport 410. Columbus, OH: Department of Geodetic Science and Surveying, Ohio State University199119.http://dx.doi.org/10.1029/96JB0164520.http://dx.doi.org/10.1007/s00190005008721.GruberTAnzenhoferMRentschMImprovements in high resolution gravity field modeling at GFZIn: Segawa J, Fujimoto H, Okubo S, editors. Gravity, geoid and marine geodesy. Berlin: Springer199744545222.LemoineFGKenyonSCFactorJKet alThe development of the joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) geopotential model EGM96NASA Technical Paper NASA/TP1998206861. Greenbelt, MD: Goddard Space Flight Center199823.http://dx.doi.org/10.1029/2000GL01158924. DobslawHSchwintzerPBarthelmesFet alGeostrophic ocean surface velocities from TOPEX altimetry, and CHAMP and GRACE satellite gravity modelsScientific technical report 04/07. Potsdam: GeoForschungs Zentrum200425. http://dx.doi.org/10.1029/2002GL01506426.http://dx.doi.org/10.1016/S0273--1177(03)00162-527. http://dx.doi.org/10.1007/s00190-005-0480-z28.FoersteCSchmidtRStubenvollRet alThe GFZ/GRGS satellite and combined gravity field models EIGEN-GL04S1 and EIGEN-GL04CJ Geodesy200682633134629.FoersteCFlechtnerFSchmidtRA new global combined high-resolution GRACE-based gravity field model of the GFZ-GRGS co-operationGeophys Res Abstr200810:EGU2008A-0342630.PavlisNKHolmesSAKenyonSCFactorJKAn Earth gravitational model to degree 2160: EGM2008Vienna: EGU General Assembly200831.http://dx.doi.org/10.1016/S0264-3707(01)00050-332. http://dx.doi.org/10.1071/EG0306933. FoersteCStubenvollRKoenigRet alEvaluation of EGM2008 by comparison with other recent global gravity field modelsNewton’s Bull20094263734. http://dx.doi.org/10.1098/rsta.2006.175135. ZhangKFFeatherstoneWEThe statistical fit of high degree geopotential models to the gravity field of AustraliaGeomatics Res Australas19956311836.PearseMBKearsleyAHWAnalysis of EGM models in New ZealandInt Geoid Serv Bull1996620321337.KirbyJFFeatherstoneWEKearsleyAHWTests of the DMA/GSFC geopotential models over AustraliaInt Geoid Serv Bull1998721338.EllmannAThe geoid for the Baltic countries determined by the least squares modification of Stokes’ formulaGeodesy Report No. 1061. Stockholm: Royal Institute of Technology, Department of Infrastructure200439. YilmazIYilmazMGulluMTurgutBEvaluation of recent global geopotential models based on GPS/leveling data over Afyonkarahisar (Turkey)Sci Res Essays20105548449340. BotaiMCCombrinckLInvestigating the accuracy of gravity field models using Satellite Laser Ranging dataS Afr J Geol. 20111143–453954441.CoxCMAuABoyJ-PChaoBFTime-variable gravity: Using Satellite-Laser-Ranging as a tool for observing long-term changes in the Earth systemIn: Noomen R, Klosko S, Noll C, Pearlman M, editors. Proceedings of the 13th International Workshop on Laser Ranging. Greenbelt, MD: NASA Goddard Space Flight Center200342.TapleyBDBettadpurSRiesJCThompsonPFWatkinsMMGRACE measurements of mass variability in the Earth systemScience200430550350543.http://dx.doi.org/10.1029/98JB0284444. http://dx.doi.org/10.1016/S0079-1946(98)00149-945.http://dx.doi.org/10.1029/2003GL01868846.http://dx.doi.org/10.1007/s00190-004-0417-y47. http://dx.doi.org/10.1038/303757a048.ttp://dx.doi.org/10.1126/science.10721881216165249.http://dx.doi.org/10.1029/1999GL00844050.http://dx.doi.org/10.1111/j.1365-246X.2006.02915.x51.EubanksTMVariations in the orientation of the Earth, in contributions of space geodesy to geodynamics: Earth dynamicsGeodyn Series volume 24. Washington DC: American Geophysical Union199315452.PonsarSDehantVHolmeRJaultDPaisAVan HoolstTThe core and fluctuations in the Earth’s rotation, in Earth’s core: Dynamics, structure, rotationWashington DC: Geodyn, American Geophysical Union200325126153. http://dx.doi.org/10.1029/2003JB00243254.EckmanMWhat 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–5155.RappRHGlobal models for the 1 cm geoid – Present status and near term prospectsSanso F, Rummel RGeodetic boundary value problems in view of the one centimeter geoid. Berlin: Springer199727331156.http://dx.doi.org/10.1007/s11200-005-0011-757.ChandlerGMerryCLThe South African geoid 2010: SAGEOID10Position IT. 2010(June):29–3358.MerryCEvaluation of global geopotential models in determining the quasigeoid for South Africa Surv Rev.20073930518019259. http://dx.doi.org/10.1029/2003RG00013960.http://dx.doi.org/10.1098/rsta.2006.17371653714061.http://dx.doi.org/10.1029/2004GL02094862.http://dx.doi.org/10.29/2005GL02357463.http://dx.doi.org/10.1007/s00190-005-0014-864.WinsemiusHCSavenijeHHGVan den HurkNCZapreevaEAKleesRAssessment of gravity recovery and climate experiment (GRACE) temporal signature over the upper ZambeziWater Resour Res2006422865. TorrenceMLageos-1, -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.CombrinckLHauptWMOBLAS-6 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