Global gravity field modeling based on GOCE and complementary gravity data

International Journal of Applied Earth Observation and Geoinformation 35 (2015) 120–127

Contents lists available at ScienceDirect

International Journal of Applied Earth Observation and Geoinformation journal homepage: www.elsevier.com/locate/jag

Global gravity field modeling based on GOCE and complementary gravity data Thomas Fecher ∗ , Roland Pail, Thomas Gruber Institute of Astronomical and Physical Geodesy, Technische Universität München, Arcisstraße 21, 80333 München, Germany

a r t i c l e

i n f o

Article history: Received 3 June 2013 Accepted 30 October 2013 Available online 2 December 2013 Keywords: GOCE Combination Gravity field model High performance computing

a b s t r a c t A combined high-resolution global gravity field model up to degree/order (d/o) 720, including error estimates in terms of a full variance–covariance matrix, is determined from GOCE (Gravity field and steady-state Ocean Circulation Explorer) and complementary gravity field data. GOCE observations, highly accurate in the low to medium wavelength part (∼d/o 40–220), are supplemented by GRACE (Gravity Recovery and Climate Experiment) with high accuracy in the low wavelength part (∼d/o 2–150), and altimetric and terrestrial gravity field observations to enhance the spectral resolution of the combined gravity field model. The theory of combining different data sets by least-squares techniques, applying optimum weighting strategies, is illustrated. Full normal equation systems are used to enable stochastic modeling of all individual observations. High performance computing techniques are applied in order to handle normal equations of enormous size (about 2 TB). The quality of the resulting gravity field solution is analyzed by comparisons with independent gravity field models and GPS/leveling observations, and also in the frame of the computation of a mean dynamic topography. The validation shows that the new combined model TUM2013C achieves the quality level of established high-resolution models. Compared to EGM2008, the improvements due to the inclusion of GOCE are clearly visible. © 2013 Elsevier B.V. All rights reserved.

1. Introduction In the framework of ESA’s dedicated gravity field satellite mission GOCE (Drinkwater et al., 2003) various kinds of gravity products are provided to the scientific community. They are generated by the GOCE High-level Processing Facility (HPF), a collaboration of 10 European scientific institutions under contract of ESA. Examples for GOCE data products are gravity gradients in different reference frames (Bouman and Fuchs, 2012) or GOCE based satellite gravity field models (Pail et al., 2011). At present the fourth generation of GOCE gravity field models have been released based on 32 months of GOCE data. The combination of GOCE with data from the US–German gravity field mission GRACE (Tapley et al., 2004), potentially complemented by CHAMP (CHAllanging Mini-satellite Payload; Reigber et al., 2002)) and satellite laser ranging (SLR) data, yields the actual state-of-the-art for satellite-only models. The GOCE models are basically enriched by the excellent performance of the GRACE information in the long wavelength part, which is shown in Fig. 1. The GOCO03s satellite-only gravity field model (compare Pail et al., 2010), a gravity model containing ITG-GRACE2010S

∗ Corresponding author. Tel.: +49 89 289 231 85. E-mail address: [email protected] (T. Fecher). 0303-2434/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jag.2013.10.005

(Mayer-Gürr et al., 2010) normal equations and 18 months of GOCE gradiometer data (as well as CHAMP and SLR information), combines the advantages of GRACE and GOCE. This can be seen from their formal errors, which are displayed in terms of degree error medians (Siemes, 2008): sn = medianm {C¯ nm ; S¯ nm }

(1)

Nevertheless, as the gravity signal is due to Newton’s law attenuated with increasing satellite height, satellite-only gravity field models are restricted in their spectral (and as a result spatial) domain. This is indicated by the blue curve in Fig. 1, where the intersection with the gray signal curve shows a signal-to-noise ratio of 1. To achieve higher spectral resolutions, which is important for various scientific but also engineering and navigation applications, terrestrial and altimetric gravity data have to be included. As Pavlis et al. (2012) point out, a global high resolution gravity model can amongst others serve as a reference and a base for the modeling of regional geoids or can provide geoid heights for the usage in the framework of orthometric height determination. In the pre-satellite gravity field mission era several gravity field models containing altimetric and terrestrial data such as the OSU or the GRIM series have been released. The most prominent model at this time was probably the EGM96 (Lemoine et al., 1998) with a spectral resolution up to degree/order (d/o) 360. Nowadays the most

T. Fecher et al. / International Journal of Applied Earth Observation and Geoinformation 35 (2015) 120–127

121

fully-normalized associated Legendre functions P¯ nm , the spherical harmonic degree n and order m, and the spherical harmonic coefficients (C¯ nm , S¯ nm ), which parameterize the gravity field model ( implies that the coefficients of the normal field are subtracted). Other relevant anomalous gravity quantities, such as the second derivatives of the anomalous potential Tij , corresponding to the gravity gradients measured by GOCE after subtracting the normal field contribution, or gravity anomalies g are derivatives or linear combinations of the derivatives of T. As an example, the radial component Tzz yields

 n+3

GM  R ∂2 T = r ∂2 r R3 ∞

Tzz (, , r) =

(n + 1)(n + 2)

n=0

×

n 

P¯ nm (cos )[C¯ nm cos(m) + S¯ nm sin(m)],

m=0

Fig. 1. Degree (error) medians of satellite-only models based on GRACE, GOCE gradiometry, and the combined model GOCO03s. (For interpretation of the references to color in the citation of this figure, the reader is referred to the web version of this article.)

(3) and the gravity anomaly in spherical approximation is

2. Theoretical background and strategy The Earth’s gravity field is generally described in terms of a spherical harmonic series through the anomalous potential n  n+1 

GM  R R r ∞

T (, , r) =

n=0

(2)

with the spherical coordinate triplet (, , r), the Earth mean radius R, the product of gravitational constant and Earth mass GM, the

(n − 1)

n=2

×

n 

P¯ nm (cos )[C¯ nm cos(m) + S¯ nm sin(m)].

m=0

(4) Since in practice the summation in Eqs. (2)–(4) cannot be infinite, the spherical harmonic series has to be truncated at a specific degree nmax . This degree nmax defines the spatial resolution of the gravity field model by the relation max =

180◦ . nmax

(5)

max describes the spatial dimension of structures on the Earth surface, so that they can be resolved by the corresponding gravity field model. Vice versa, the spatial resolution of the observations defines the maximal degree nmax which can be determined from an observed data set, representing the well-known Nyquist theorem on the sphere (e.g. Colombo, 1981). coefficients, The optimal estimation of the spherical harmonic   also referred to as unknown parameters x = C¯ nm , S¯ nm , from the observations l is the main goal of gravity field modeling. As method for the determination of the unknowns least squares techniques based on a Gauß–Markov model are applied, where the functionals of the anomalous potential serve as the observations l. In the functional model l + v = Aˆx

(6)

v is the vector of the residuals and A the so-called design matrix. It contains the partial derivatives of the spherical harmonic series, and provides the linear (or linearized) connection of the unknowns and the observations. The accuracy information of the observations is summarized in the covariance matrix ll of the observations and is referred to as stochastic model. Correspondingly the weight matrix of the observations is P=

P¯ nm (cos )[C¯ nm cos(m)

m=0

+ S¯ nm sin(m)]

g(, , r) = −

 n+2

2 ∂T GM  R − T= 2 R r ∂r R ∞

prominent high-resolution gravity field model is the EGM2008 (Pavlis et al., 2012) with a resolution up to degree 2190 and order 2159. Beside EGM2008 GFZ’s and CNES’ EIGEN series is established in the field of global high resolution gravity field modeling (Förste et al., 2008), where the newest EIGEN models, such as the EIGEN6C2 (Förste et al., 2011), contain, compared to EGM2008, not only GRACE but also GOCE information. GOCE products contribute to several scientific disciplines. The GOCE geoid serves as important reference surface in oceanography, allowing the calculation of mean dynamic topography and ocean circulation (Knudsen et al., 2011) or the modeling of geostrophic currents (Bingham et al., 2011; Bingham, 2015). Also for the unification of global height systems GOCE plays an important role (Rummel, 2013; Woodworth et al., 2013; Gatti et al., 2013; Gerlach and Fecher, 2013). In the field of geophysics GOCE products can be used to get insight into geophysical and geological structures of the Earth’s interior (Braitenberg et al., 2011; Braitenberg, 2015) as for example the Moho determination (Reguzzoni and Sampietro, 2012; Sampietro et al., 2013; van der Meijde et al., 2013), and lithospheric modeling (Hosse et al., 2013; Bouman et al., 2015). This manuscript describes the combination of GOCE satellite data with complementary gravity field information applying optimal weighting, in order to derive a high-resolution gravity field model with a spectral resolution of d/o 720. Section 2 provides an overview of the basic theory, the functional model and the combination strategy, followed by Section 3 with a description of the complementary data sets as well as the stochastic modeling. Section 4 deals with the computational challenges of high-resolution gravity field modeling based on full normal equations, while Section 5 presents the results of a global combined gravity field model. Finally Section 6 shows validation results for the computed model by independent data sources.

1 02

−1

ll

(7)

with 02 denoting the a priori variance. The solution, which fulfills the condition vT Pv = min

(8)

122

T. Fecher et al. / International Journal of Applied Earth Observation and Geoinformation 35 (2015) 120–127

is regarded as optimal solution under the L2-norm in the least squares adjustment process. It is given by xˆ = (APA)−1 APl = N−1 n

(9)

where N is denoted as normal equation matrix, and n as right-hand side. The residuals v can be used to determine a posteriori the variance ˆ 02 =

vT Pv lT Pl − lT PAx = r r

(10)

with r denoting the redundancy (number of observations minus number of unknowns), which can be used to obtain the covariance matrix of the unknowns ˙xx = ˆ 02 N−1

xˆ =

1 ˆ 02 1

x = N−1 ni i

AT1 P1 A1 + · · · +

−1 

1 ˆ 02 k

ATk Pk Ak

1 ˆ 02 1

AT1 P1 l1 + · · · +

1 ˆ 02 k

2 −1 ll  = 02  P−1  = 0 i Ni .

,

= N−1 n

(12) where k indicates the number of different data sets respectively normal equations. ˆ 02 i are now denoted as variance components (Koch and Kusche, 2002) which are responsible for the correct relative weighting of the different data sets and are determined by an iterative procedure. Usually, it requires several iterations until

 xˆ =

1 ˆ 02



×

AT1 P1 A1 + · · · + 1

1 ˆ 02

ˆ 02

AT1 P1 l1 + · · · + 1

ˆ 02



=

vT P v

P 



1 ˆ 02

P x

.

(18)



r

T

=



r = u − trace

(ˆx − x ) P (ˆx − x ) r

1 02 

(19)

 P N−1

,

(20)

where u is the number of prior information elements. In our combination scheme, the normal equations of the satellite gravity field missions GRACE and GOCE are included as “prior information” via this formulation. Altimetric and terrestrial data is given as gravity anomalies, but not in form of discrete point values, as the Eqs. (2)–(4) imply, but in form of block-mean values. The modified Eq. (4) for block-means (Gruber, 2000) reads

 n+2

1 GM  R r S R2 ∞

m

−1

1

The variance component of the prior information can be estimated by

g(1 , 2 , 1 , 2 , r) = n 

(17)

The original normal equation matrix Ni serves as weight matrix for the prior information. Extending Eq. (12) by this additional component expressed by Eqs. (15) and (17) yields

with

 ATk Pk lk

(16)

and the stochastic model

(11)

containing the accuracy information of the estimated coefficients and their correlations among each other. The combination procedure corresponds to a summation of the full normal equations of the individual components, assuming uncorrelated observation groups, i.e., in our case terrestrial, altimetric and satellite data sets



with the prior information

(n − 1)·

n



2 2 2 + 1  2 + 1  PI nm (1 , 2 ) C¯ nm cos m sin m + S¯ nm sin m sin m m 2 2 m 2 2



(21)

with the final result is achieved. The covariance matrix of the unknowns is now the inverse of the sum of all correctly weighted normal equations. The variance components of the individual observations groups can be estimated analogously to (10) by ˆ 02

i

=

vTi Pi vi ri

T

=

T

xˆ Ni xˆ − 2ˆx ni + lTi Pi li , ri

(13)

with their redundancy factor



ri = ni − trace

1 02 i



−1

Ni N

,

(14)

where ni is the number of observations of the individual observation groups and Ni = ATi Pi Ai . Especially in the case of the satellite gravity field missions, frequently the individual normal equations have been processed by external institutions, so that the originally used observations li and their ni are unknown. In this case these normal equations can be used directly, and the formulation has to be modified according to Koch and Kusche (2002), by interpreting these already assembled normal equations as prior information. The functional model referred only to the prior information reads x + v = I xˆ

(15)

S = (cos 1 − cos 2 ) = (2 − 1 )(cos 1 − cos 2 )

(22)

and PI nm denoting the integral of the fully-normalized associated Legendre functions within the boundaries  1 and  2 . Full normal equations become quite large with increasing degree. Their (square) dimension d is related to the maximum degree by d = (nmax + 1)2 . The corresponding normal equations of a satellite-only gravity field model such as the GOCO03s have already a dimension of 63,000 × 63,000 elements and require a working memory of 30 GB, for a model up to d/o 360 the storage of the normal equations requires already 125 GB. As the gravity data at our disposal have at least a spatial resolution of 15 , our model’s spectral resolution should be, according to Eq. (5), d/o 720, and the corresponding normal equations require a working memory of 2 TB. As in former times computational resources were limited, there exist several strategies to avoid full normal equations, by making use of the sparse matrix structure, such as block-diagonal techniques. They are presented for example in Colombo (1981), Sneeuw (1994), Gruber (2001), Pavlis et al. (2012), Reguzzoni and Sansò (2012). To attain a block-diagonal structure, besides an order-dependent ordering scheme of the normal equation elements three conditions have to be fulfilled: (i) the data set has to be consistent and complete, (ii) the weighting of the individual observations must be

T. Fecher et al. / International Journal of Applied Earth Observation and Geoinformation 35 (2015) 120–127

123

Fig. 3. Fill-in strategy for the leftover areas.

Fig. 2. Block-mean value gravity anomaly data sets. (For interpretation of the references to color in the citation of this figure, the reader is referred to the web version of this article.)

independent of the longitude, and (iii) the weighting of the individual observations must be equator-symmetric. However, since we want to enable an individual stochastic modeling of all components and even assign individual weights to every block-mean value of the terrestrial and satellite altimetry data, these conditions cannot be fulfilled, and thus full normal equations have to be used. 3. Data sets The satellite gravity components of TUM2013C (TUM = Technische Universität München; C = combined) include GRACE and GOCE-gradiometry normal equations, as they have also been included, among other components, in the GOCO03S satellite-only gravity field model. In the case of GRACE, these are the ITG-GRACE2010S (Mayer-Gürr et al., 2010) normal equations, which base on a data period from August 2002 to August 2009. The GOCE-gradiometry normal equations include GOCE gradients from November 2009 to April 2011. In addition to these satellite normal equations, a global 15 min block mean value gravity anomaly data set is generated through homogenization of different data sets (Fig. 2 and Table 1). Over the oceans the DTU10 gravity anomalies are used, and ArcGP gravity anomalies in the Northern polar areas. Moreover we gained access to several not freely available gravity anomaly data sets: these are EGG2008 (Europe), Ausgeoid09 (Australia), and block-mean value gravity anomalies from the United States and Canada. The spatial resolution of single data sets in some areas was higher than 15 , and therefore they were down-sampled to 15 . The spectral signal content above d/o 720 was reduced from all data sets by applying filtering in spectral domain. As Fig. 2 shows, some areas are remaining, where no terrestrial or altimetric data sets are at our disposal (in green: leftover). Since the normal equation system would become unstable, if we just neglect these areas, a fill-in strategy was applied, which does not make use of synthetic gravity anomalies from other global gravity field models (e.g. EGM2008). The fill-in strategy is illustrated in Fig. 3 and Table 1.

For filling the gaps, we decided to use NIMA96 30 block-mean value gravity anomalies (Kenyon and Pavlis, 1996), even if the quality of the data set is in general not as good as of the data sets mentioned before. NIMA96 was the terrestrial database used for the EGM96 model (Lemoine et al., 1998). The main problem of the NIMA data set is that due to its spatial resolution of 30 the corresponding spectral resolution is limited to d/o 360, compare Eq. (5). Therefore the gap between d/o 360 and d/o 720 has to be bridged. This can be done by topographic gravity anomalies, compare Pavlis et al. (2012). Topographic gravity anomalies can be synthesized from digital terrain models by the calculation of the gravitational attraction of rectangular prisms, which represent the topographic masses, at surface points (Hirt and Kuhn, 2012). Topographic anomalies are calculated from the global topographic/isostatic model RWI TOIS2012 (Grombein et al., 2013), which is parameterized in terms of spherical harmonic coefficients up to d/o 1800. For this work, we apply only the topographic part up to a maximum d/o of 720. To avoid an edged transition at d/o 360, we use topographic anomalies with signal content from d/o 120 up to d/o 720 and set up the normal equations correspondingly. Thus we achieve a weighted average of NIMA and RWI TOIS2012 in the fill-in regions. As Fig. 3 shows, there are few areas, where no NIMA data is available. Here the gaps are filled with 15 block-mean value gravity anomalies from a regularized GRACE/GOCE combination, whose GRACE and GOCE normal equations correspond to the normal equations used for TUM2013C. Combining all these gravity anomalies of the different regions, a global coverage of gravity anomalies is achieved. At this point it should be mentioned, that for all calculations the heights that were delivered together with the gravity anomaly values have been used, and correspondingly the normal equations refer to topography (variable r in Eq. (21)). In the estimation process all data sets are kept separately as observation groups, and normal equations are assembled individually. Finally, they are optimally combined according to the description in Section 2. The accuracy of the data sets is the basis for the stochastic modeling. Fig. 4 shows the errors assigned to the block-mean values, forming the variance–covariance matrix ll of the terrestrial gravity anomaly grid, which is assumed to be diagonal, because no information on correlations among the observations is available.

Table 1 Gravity anomaly data sets. Data set ArcGP Ausgeoid09 – DTU10 EGG2008 NIMA96 RWI TOIS2012 –

Region Northern polar land areas Australia Canada Ocean areas Europe Leftover Leftover USA

Data provider ArcGP team Featherstone and coworkers Veronneau and Huang Anderson and Knudsen Denker Lemoine and coworkers Grombein and Seitz Kenyon

Original resolution 

5 1 2 1 15 30 15 15

Min. degree 0 0 0 0 0 0 120 0

Max. degree 720 720 720 720 720 360 720 720

124

T. Fecher et al. / International Journal of Applied Earth Observation and Geoinformation 35 (2015) 120–127

Fig. 4. Errors assigned to the block-mean values from Fig. 2 (mGal). (For interpretation of the references to color in the citation of this figure, the reader is referred to the web version of this article.)

In case of the DTU10 gravity anomalies the metric DTU10 error information (sea surface height errors), resulting from altimetric applications, is linearly transformed to gravity anomaly errors. For the US and European data sets error information is available. For all other data sets, where no error information is available, empirical error information can be obtained by comparison with EGM2008. As the US authorities have the widest and probably best terrestrial data material, which is partially restricted and not freely accessible, but has entered in wide parts the EGM2008, it can be used to estimate the quality of terrestrial data sets. For topographic gravity anomalies empirical values were used as stochastic model. These initial error assumptions are modified by the variance components during the iterative estimation process, yielding the final stochastic model. Fig. 4 shows the final errors for the terrestrial data grid after variance component estimation, which included also the satellite components GRACE and GOCE. The final weights for the fill-in values of NIMA and the topographic anomalies, which exceed the range of the color bar of Fig. 4, are 6 < =  NIMA < = 22 mGal and  TOPO = 20 mGal, respectively. Since the normal equations are numerically stable, no regularization has to be applied.

libraries. Additionally BLACS, a linear algebra orientated message passing interface, is applied. In the program for accumulating the normal equations, each of the up to 3500 total cores calculates independently a sub-matrix of the normal equation system, which are then stored on the disk. The calculation of the sub-matrices is performed by dyadic products of the rows of the design matrix by the BLAS library, which is integrated in MKL. The runtime of this program is varying depending on the number of observations. The program which combines the individual normal equation systems uses MPI-IO2 routines for an effective reading and writing. Two 2 TB normal equations are combined by up to 3500 cores in less than 1 min. The final step to achieve a gravity field solution, i.e. solution and inversion of the normal equations, is performed by a separate program using SCALAPACK routines which are also integrated in MKL, and has to be carried out by all cores together, because extensive communication among the nodes is required. Applying 2000 cores a 2 TB normal equation system can be solved and inverted in less than 3 h. 5. Combined global gravity field model up to d/o 720 The data sets described in Section 3 can be processed based on the theory formulated in Section 2, resulting, after several iteration steps at LRZ’s SuperMUC, in a combined high resolution global gravity field model up to d/o 720 based on full normal equations. Fig. 5 shows the results in terms of degree variances 1  2 2 [C¯ nm + S¯ nm ], 2n + 1 n

n2 =

(23)

m=0

4. Computational challenge

which are used here in contrary to the degree medians in Eq. (1), as none of the displayed solutions which are based on global data sets is affected by the polar gap (the blue curve of GOCE gradiometry in Fig. 1 is affected by the polar gap and is therefore difficult to compare with the other solutions in term of degree variances). Differences of several gravity models to the EGM2008 signal are displayed. The green curve, the difference of the new solution TUM2013C to GOCO03S indicates that the combined solution agrees with GRACE and GOCE in the low degrees, but the influence of the satellite information is decreasing with increasing degree. In contrast, the blue and the orange curve, the differences of the new solution TUM2013C and GOCO03S to EGM2008 highlight the increasing influence of the terrestrial data with higher degrees. This

Gravity field modeling using full normal equation systems is a computational challenge. The size of the full normal equations of a gravity field model resolved up to d/o 720 is about 2 TB. Besides hard disk storage and some other issues, the most challenging task is to keep the normal equation system in the working memory. This can only be achieved by using massive parallel techniques and supercomputers. Here gravity field modeling crosses the field of high performance computing (compare e.g., Wittwer, 2006; Fecher et al., 2011; Brockmann and Schuh, 2011). The computations for this work were performed on SuperMUC. SuperMUC, operated by the Leibniz Supercomputing Center in Garching near Munich and fully operational since summer 2012, is equipped with 155,000 cores, and achieves a peak performance of 3 Petaflop/s. A general parallel file system (GPFS) allows the hard disk storage of several full normal equation systems and, as the name implies, a parallel and hence fast access to the normal equations. Our calculations are performed on Sandy Bridge-EP Intel Xeon processors with 2 GB memory per core. Our applications run with IBM’s MPI implementation, and Intel’s Math Kernel Library (MKL) provides all required mathematical

Fig. 5. Difference of TUM2013C to other models in terms of degree variances.

T. Fecher et al. / International Journal of Applied Earth Observation and Geoinformation 35 (2015) 120–127

125

Table 2 Orbit validation: RMS of 12 h GOCE orbit fits (m) based on different gravity field models up to d/o 180.

2013-01-01 2013-01-02 2013-01-03 2013-01-04 2013-01-05 2013-01-06

TUM2013C

ITG-GRACE 2010s

EGM2008

EIGEN-6C2

0.12 0.07 0.05 0.06 0.10 0.08

0.13 0.07 0.06 0.06 0.11 0.08

0.14 0.09 0.07 0.06 0.12 0.10

0.12 0.07 0.06 0.06 0.10 0.08

no high accurate gravity data is available, benefit most from the satellite information. 6. Model validation and quality analysis

Fig. 6. Difference of TUM2013C and EGM2008 in terms of geoid heights (m).

result is in agreement with the theoretical considerations in Section 2. The systematic difference of the new combination model TUM2013C to EGM2008 in terms of geoid heights (Fig. 6) shows the typical structures in South America, Africa and Asia, which are caused by the improvements from GOCE data (e.g., Pail et al., 2011). These improvements underline the importance of the GOCE mission especially for several continental regions, where no highquality terrestrial gravity measurements are available. In areas where the terrestrial data sets are quite accurate, as in North America, Europe or Australia, nearly no differences are visible. This proves that GOCE is consistent with the high-accuracy terrestrial measurements in these regions. Besides the improvements from GOCE some bubble structures can be seen in South America, Asia and Africa, resulting from the limited performance of the fill-in data sets. Based on the full variance–covariance matrix of the combined gravity field solution, which has a dimension of 519,841 × 519,841 elements, covariance propagation to geoid heights can be performed. If this is done separately for the satellite-only contribution and the combined model, the ratio of both results corresponds to redundancy factors, which express the relative contribution of the satellite information to the combined solution. This ratio is shown in Fig. 7 at d/o 190. Again it is evident, that those regions, where

In previous sections we have shown, how GOCE information can be supplemented by additional gravity field information to determine high resolution combination gravity field models. Here we assess the quality of our combination model TUM2013C. In this quality analysis we compare the new model to the established high resolution gravity field models EGM2008 and EIGEN-6C2. The first step of the quality analysis is an orbit validation. We consider the RMS of 12 h GOCE orbit fits up to d/o 180. This means that an orbit is fitted to kinematic positions. The results are summarized in Table 2. In addition to the high resolution models the RMS of the GRACE-only solution ITG-GRACE2010s is shown, as GRACE performs best in the very low harmonic degrees, which have very strong influence on satellite orbits. The RMS of our solution is smaller than the RMS of ITGGRACE2010s and EGM2008, demonstrating the positive influence of the GOCE information. Compared to EIGEN-6C2, the result is nearly identical, because both models include GOCE information. Another strategy of external validation is the comparison with independent “direct geoid observations” derived from long-term GPS observations (delivering ellipsoidal heights h) and spirit leveling (providing orthometric heights H), because their difference is the geoid: N = H − h. Following the methodology described in Gruber et al. (2011), the models were truncated at a certain maximum degree and order n, and filled up beyond with the combined gravity field model EGM2008 (Pavlis et al., 2012) in order to reduce omission errors. It should be emphasized, that this approach gives EGM2008 for these comparisons an advantage, because it is the only model which is consistently used up to its maximum degree. Table 3 summarizes the RMS of geoid differences at d/o 720 for different regions where GPS/leveling data sets are available. The results show that TUM2013C achieves the same quality level as the reference models. In some areas, where terrestrial input data sets show insufficiencies and systematic errors, there is still some room for improvement. Errors can be seen for example in the Australian and European data sets with a magnitude of up to 15 respectively 25 mGal in regions with rough topography. They have been identified to result from deficiencies of the processing of the input block mean values. These data sets are actually under investigation by the producers. GPS-leveling comparisons at d/o 180 (not Table 3 RMS of geoid height differences (mm) between gravity field models and GPS/leveling observations in selected regions (number of stations in brackets).

Fig. 7. Relative contribution of satellite information to the combined solution at d/o 190 in %.

Australia (197) Germany (675) Canada (430) Europe (1233) USA (5168) UK (177)

TUM2013C

EGM2008

EIGEN-6C2

0.243 0.040 0.102 0.213 0.265 0.064

0.240 0.034 0.100 0.210 0.259 0.060

0.237 0.045 0.101 0.216 0.259 0.068

126

T. Fecher et al. / International Journal of Applied Earth Observation and Geoinformation 35 (2015) 120–127

Fig. 8. Mean dynamic topography (m) derived from sea surface height GFSC00.2 and different gravity field models up to d/o 720.

shown here) demonstrate moreover, that the good performance of GOCE is not negatively affected by the combination with terrestrial data. This is an indication for a good stochastic modeling and relative weighting between satellite components and the terrestrial grid. Finally, mean dynamic topographies (MDTs) were calculated based on all high-resolution gravity field models. The mean dynamic topography is the difference between the mean sea surface (here: GSFC00.2 (Wang, 2001)) and the geoid height. Fig. 8 shows the results for an evaluation of the MDT up to d/o 720 (see Gruber, 2009). It demonstrates that the quality of TUM2013C is comparable to EGM2008, while EIGEN-6C2 shows bumpy features which are suspected to be numerical artifacts resulting from a non-optimum combination. 7. Conclusions and outlook In this study GOCE was combined with complementary gravity field observations using variance-component estimation and least squares techniques based on full normal equations, resulting in the combined gravity field model TUM2013C. It is resolved up to d/o 720 and delivered together with the full parameter variance–covariance matrix. We put emphasis on a realistic stochastic modeling of all components. The use of full normal equation systems enables an individual weighting of all observations and components. LRZ’s supercomputing system SuperMUC was used to handle these very large normal equation systems. As data sources only GRACE, GOCE and available terrestrial and

altimetric data sets have been used, but no external gravity field models or gravity anomalies synthesized from them. Therefore, with the resulting gravity model TUM2013C for the very first time a high-resolution combined model, which is totally independent of EGM2008, has been produced. Validation shows that our model achieves the quality level of the established high resolution gravity field models. Compared to EGM2008, the improvements due to the inclusion of GOCE data are clearly visible, and TUM2013C outperforms EIGEN-6C2 in the MDT validation. An inclusion of additional GOCE data sets is planned for the future. The next generation of GOCE gravity field models will persuade by amazing improvements, due to the reprocessing of the GOCE gravity gradients (Stummer et al., 2012), and the lowering of the satellite’s altitude (Floberghagen, 2013). Including these data sets will also push combination models a next step forward. Further investigations on the preprocessing of the available terrestrial gravity anomaly data sets and the elimination of systematic inconsistencies, which could be identified during this work, shall also enable further improvements. Acknowledgements We gratefully acknowledge computational resources and support granted by Leibniz Supercomputing Center under project pr32qu. Steve Kenyon, Heiner Denker, Christian Hirt, Ole Anderson, Marc Verenneau and Jianling Huang, Thomas Grombein and Kurt Seitz gave us access to their data sets, what we cordially appreciate. This work would not be possible without the amazing work of the

T. Fecher et al. / International Journal of Applied Earth Observation and Geoinformation 35 (2015) 120–127

GOCO project team (colleagues from TU Graz, Uni Bonn, Uni Bern and ÖAW). We thank Markus Heinze for calculating the orbit validation results. We kindly acknowledges the support given by the Institute for Advanced Study (IAS)/TUM. References Bingham, R.J., Knudsen, P., Andersen, O., Pail, R., 2011. An initial estimate of the North Atlantic steady-state geostrophic circulation from GOCE. Geophysical Research Letters 38 (1), L01606. Bingham, R., 2015. A comparison of GOCE and drifter-based estimates of the North Atlantic steady-state surface circulation. Int. J. Appl. Earth Obs. Geoinform. 35, 140–150. Bouman, J., Fuchs, M., 2012. GOCE gravity gradients versus global gravity field models. Geophysical Journal International 189 (2), 846-850, http://dx.doi.org/10.1111/j.1365-246X.2012.05428.x. Bouman, J., Ebbing, J., Meekes, S., Fattah, R., Fuchs, M., Gradmann, S., Haagmans, R., Lieb, V., Schmidt, M., Dettmering, D., Bosch, W., 2015. Goce gravity gradient data for lithospheric modelling. Int. J. Appl. Earth Obs. Geoinform. 35, 16–30. Braitenberg, C., 2015. Goce observations for mineral exploration in Africa and across continents. Int. J. Appl. Earth Obs. Geoinform. 35, 88–95. Brockmann, J.M., Schuh, W.-D., 2011. Use of Massive Parallel Computing Libraries in the Context of Global Gravity Field Determination from Satellite Data. In: Ouwehand, L. (Ed.), Proceedings of the 4th international GOCE User Workshop. ESA Publication SP-696, ESA/ESTEC. , ISBN 978-92-9092-260-5, ISSN 1609-042X. Colombo, O.L., 1981. Numerical Methods for Harmonic Analysis on the Sphere, OSU Report 310. Dept. of geodetic Science, The Ohio State University, Columbus, OH. Drinkwater, M.R., Floberghagen, R., Haagmans, R., Muzi, D., Popescu, A., 2003. GOCE: ESA’s first Earth Explorer Core mission. In: Beutler, G., Drinkwater, M.R., Rummel, R., von Steiger, R. (Eds.), Earth Gravity Field from Space – From Sensors to Earth Sciences. Space Sciences Series of ISSI, vol. 17. Kluwer Academic Publishers, Dordrecht, The Netherlands, pp. 419–432, ISBN: 1-4020-1408-2. Fecher, T., Gruber, T., Pail, R., 2011. The gravity field – an important parameter for Earth observation. inSiDE 9 (2), 26–31, Gauss Centre for Supercomputing (HLRS, LRZ, JSC). Floberghagen, R., 2013. International Journal of Applied Earth Observation and Geoinformation (in this issue). Förste, C., Bruinsma, S., Shako, R., Marty, J.-C., Flechtner, F., Abrikosov, O., Dahle, C., Lemoine, J.-M., Neumayer, K.H., Biancale, R., Barthelmes, F., König, R., Balmino, G., 2011. EIGEN-6 – a new combined global gravity field model including GOCE data from the collaboration of GFZ-Potsdam and GRGS-Toulouse. Geophysical Research Abstracts 13, EGU2011-3242-2, EGU General Assembly. Förste, C., Flechtner, F., Schmidt, R., Stubenvoll, R., Rothacher, M., Kusche, J., Neumayer, K.H., Biancale, R., Lemoine, J.-M., Barthelmes, F., Bruinsma, S., König, R., Meyer, U., 2008. EIGEN-GL05C – a new global combined high-resolution GRACEbased gravity field model of the GFZ-GRGS cooperation. Geophysical Research Abstracts 10, EGU2008-A-03426, SRef-ID: 1607-7962/gra/EGU2008-A-03426. Gatti, A., Reguzzoni, M., Venuti, G., 2013. The height datum problem and the role of satellite gravity models. Journal of Geodesy 87 (1), 15–22, http://dx.doi.org/10.1007/s00190-012-0574-3. Gerlach, C., Fecher, T., 2013. Approximations of the GOCE error variance–covariance matrix for least-squares estimation of height datum offsets. Journal of Geodetic Science 2012, 247–256, http://dx.doi.org/10.2478/v10156-011-0049-0, Nr. 2, Heft 4, ISSN 2081-9943. Grombein, T., Seitz, K., Heck, B., 2013. Topographic-isostatic reduction of GOCE gravity gradients. In: Proceeding of the XXV General Assembly of the International Union of Geodesy and Geophysics, Melbourne, Australia, International Association of Geodesy Symposia. Gruber, T., 2000. Hochauflösende Schwerefeldbestimmung aus Kombination von terrestrischen Messungen und Satellitendaten über Kugelfunktionen. GeoForschungsZentrum Potsdam, Dissertation, Scientific Technical Report STR 00/16, ISSN 1610-0956. Gruber, T., 2001. High resolution gravity field modeling with full variance–covariance matrices. Journal of Geodesy 75 (9–10), 505–514, http://dx.doi.org/10.1007/s001900100202, Springer, ISSN 0949-7714. Gruber, T., 2009. Evaluation of the EGM2008 gravity field by means of GPSlevelling and sea surface topography solutions. External quality evaluation reports of EGM08. Newton’s Bulletin (4), 3–17, Bureau Gravimétrique International (BGI)/International Geoid Service (IGeS), ISSN 1810-8555. Gruber, T., Visser, P.N.A.M., Ackermann, C., Hosse, M., 2011. Validation of GOCE gravity field models by means of orbit residuals and geoid comparisons. Journal of Geodesy 85 (11), 845–860, http://dx.doi.org/10.1007/s00190-011-0486-7, Springer, ISSN 0949-7714. Hirt, C., Kuhn, M., 2012. Evaluation of high-degree series expansions of the topographic potential to higher-order powers. Journal of Geophysical Research – Solid Earth 117, B12407, http://dx.doi.org/10.1029/2 2012JB009492. Hosse, M., Pail, R., Horwath, M., Köther, N., Gutknecht, B.D., 2013. Combined regional gravity model of the Andean convergent subduction zone and its application to

127

lithospheric modelling in active plate margins. Surveys of Geophysics (submitted for publication). Kenyon, S.C., Pavlis, N.K., 1996. The development of a global surface gravity data base to be used in the joint DMA/GSFC geopotential model. In: Proceedings of IAG Symposium No. 116, Global Gravity Field and its Temporal Variations, Springer Verlag. Knudsen, P., Bingham, R., Andersen, O., Rio, M.H., 2011. A global mean dynamic topography and ocean circulation estimation using a preliminary GOCE gravity model. Journal of Geodesy 85 (11), 861–879. Koch, K.R., Kusche, J., 2002. Regularization of geopotential determination from satellite data by variance components. Journal of Geodesy 76, 259–268. Lemoine, F.G., Kenyon, S.C., Factor, J.K., Trimmer, R.G., Pavlis, N.K., Chinn, D.S., Cox, C.M., Klosko, S.M., Luthcke, S.B., Torrence, M.H., Wang, Y.M., Williamson, R.G., Pavlis, E.C., Rapp, R.H., Olson, T.R., 1998. The Development of the Joint NASA GSFC and the National Imagery and Mapping Agency (NIMA) Geopotential Model EGM96. NASA Goddard Space Flight Center, Greenbelt, MD, 20771 USA. Mayer-Gürr, T., Eicker, A., Kurtenbach, E., Ilk, K.-H., 2010. ITG-GRACE: global static and temporal gravity field models from GRACE data. In: Flechtner, F., Gruber, T., Güntner, A., Mandea, M., Rothacher, M., Schöne, T., Wickert, J. (Eds.), System Earth via Geodetic-Geophysical Space Techniques. , pp. 159–168, http://dx.doi.org/10.1007/978-3-642-10228-8 13. Pail, R., Bruinsma, S., Migliaccio, F., Förste, C., Goiginger, H., Schuh, W.-D., Höck, E., Reguzzoni, M., Brockmann, J.M., Abrikosov, O., Veicherts, M., Fecher, T., Mayrhofer, R., Krasbutter, I., Sansò, F., Tscherning, C.C., 2011. First GOCE gravity field models derived by three different approaches. Journal of Geodesy 85 (11), 819–843, http://dx.doi.org/10.1007/s00190-011-0467-x, Springer, ISSN 09497714. Pail, R., Goiginger, H., Schuh, W.-D., Höck, E., Brockmann, J.M., Fecher, T., Gruber, T., Mayer-Gürr, T., Kusche, J., Jäggi, A., Rieser, D., 2010. Combined satellite gravity field model GOCO01S derived from GOCE and GRACE. Geophysical Research Letters 37, EID L20314, http://dx.doi.org/10.1029/2010GL044906, American Geophysical Union, ISSN 0094-8276. Pavlis, N.K., Holmes, S.A., Kenyon, S., Factor, J.K., 2012. The development and evaluation of the earth gravitational model 2008 (EGM2008). Journal of Geophysical Research, http://dx.doi.org/10.1029/2011JB008916. Reguzzoni, M., Sansò, F., 2012. On the combination of high-resolution and satellite-only global gravity models. Journal of Geodesy 86 (6), 393–408, http://dx.doi.org/10.1007/s00190-011-0526-3. Reguzzoni, M., Sampietro, D., 2012. Moho estimation using GOCE data: a numerical simulation. In: Kenyon, S., Pacino, M.C., Marti, U. (Eds.), Proceedings of the 2009 IAG Symposium. Buenos Aires, Argentina, 31 August–4 September 2009, IAG Symposia 136, Springer Berlin Heidelberg, http://dx.doi.org/10.1007/978-3-642-20338-1 25. Reigber, C., Balmino, G., Schwintzer, P., Biancale, R., Bode, A., Lemoine, J.M., Koenig, R., Loyer, S., Neumayer, H., Marty, J.C., Barthelmes, F., Perossanz, F., 2002. A high quality global gravity field model from CHAMP GPS tracking data and accelerometry (EIGEN-1S). Geophysical Research Letters 29, 14, http://dx.doi.org/10.1029/2002GL015064. Rummel, R., 2013. Height unification using GOCE. Journal of Geodetic Science 2012 (2), 355–362, http://dx.doi.org/10.2478/v10156-011-0047-2, Heft 4, Versita, ISSN 2081-9943. Sampietro, D., Reguzzoni, M., Braitenberg, C., 2013. The GOCE estimated Moho beneath the Tibetan Plateau and Himalaya. In: Rizos, C., Willis, P. (Eds.), Proceedings of the IAG General Assembly. Melbourne, Australia, June 28–July 2 2011, IAG Symposia 139, Springer Berlin Heidelberg. Siemes, C., 2008. Digital Filtering Algorithms for Decorrelation Within Large Least Squares Problems. Universität Bonn, Dissertation, urn:nbn:de:hbz:5N-13749. Sneeuw, N., 1994. Global spherical harmonic analysis by least-squares and numerical quadrature methods in historical perspective. Geophysical Journal International 118 (3), 707–716, http://dx.doi.org/10.1111/j.1365-246X.1994.tb03995.x, ISSN 0956-540X. Stummer, C., Siemes, C., Pail, R., Frommknecht, B., Floberghagen, R., 2012. Upgrade of the GOCE Level 1b gradiometer processor. Advances in Space Research 49 (4), 739–752, http://dx.doi.org/10.1016/j.asr.2011.11.027. Tapley, B.D., Bettadpur, S., Watkins, M., Reigber, C., 2004. The gravity recovery and climate experiment: mission overview and early results. Geophysical Research Letters 31 (9), L09607, http://dx.doi.org/10.1029/2004GL019920, American Geophysical Union. van der Meijde, M., Julia, J., Assumpcao, M., 2013. Gravity derived Moho for South America. Tectonophysics: The International Journal of Integrated Solid Earth Sciences. Wang, Y.M., 2001. GSFC00 mean sea surface, gravity anomaly, and vertical gravity gradient from satellite altimetry. Journal of Geophysical Research 106, 31167–31174. Wittwer, T., 2006. An Introduction to Parallel Programming. Textbooks from VSSD., ISBN 978-90-71301-78-0. Woodworth, P.L., Hughes, C.W., Bingham, R.J., Gruber, Th., 2013. Towards worldwide height system unification using ocean information. Journal of Geodetic Science 2012 (2), 302–318, http://dx.doi.org/10.2478/v10156-012-0004-8, Heft 4, Versita, ISSN 2081-9943.