Improved daily GRACE gravity field solutions using a Kalman smoother

Journal of Geodynamics 59–60 (2012) 39–48

Contents lists available at SciVerse ScienceDirect

Journal of Geodynamics journal homepage: http://www.elsevier.com/locate/jog

Improved daily GRACE gravity field solutions using a Kalman smoother E. Kurtenbach a,∗ , A. Eicker a , T. Mayer-Gürr a , M. Holschneider b , M. Hayn b , M. Fuhrmann b , J. Kusche a a b

Institute of Geodesy and Geoinformation, Bonn University, Bonn, Germany Institute of Mathematics, Potsdam University, Bonn, Germany

a r t i c l e

i n f o

Article history: Received 30 September 2010 Received in revised form 6 February 2012 Accepted 12 February 2012 Available online 22 February 2012 Keywords: GRACE Daily gravity field Kalman smoother ITG-Grace2010

a b s t r a c t Different GRACE data analysis centers provide temporal variations of the Earth’s gravity field as monthly, 10-daily or weekly solutions. These temporal mean fields cannot model the variations occurring during the respective time span. The aim of our approach is to extract as much temporal information as possible out of the given GRACE data. Therefore the temporal resolution shall be increased with the goal to derive daily snapshots. Yet, such an increase in temporal resolution is accompanied by a loss of redundancy and therefore in a reduced accuracy if the daily solutions are calculated individually. The approach presented here therefore introduces spatial and temporal correlations of the expected gravity field signal derived from geophysical models in addition to the daily observations, thus effectively constraining the spatial and temporal evolution of the GRACE solution. The GRACE data processing is then performed within the framework of a Kalman filter and smoother estimation procedure. The approach is at first investigated in a closed-loop simulation scenario and then applied to the original GRACE observations (level-1B data) to calculate daily solutions as part of the gravity field model ITG-Grace2010. Finally, the daily models are compared to vertical GPS station displacements and ocean bottom pressure observations. From these comparisons it can be concluded that particular in higher latitudes the daily solutions contain high-frequent temporal gravity field information and represent an improvement to existing geophysical models. © 2012 Elsevier Ltd. All rights reserved.

1. Introduction Since March 2002 the twin satellite mission GRACE (Tapley et al., 2004) has measured the Earth’s gravity field and its temporal variations with unprecedented accuracy. Various products and solutions have been presented in recent years. The standard procedure is the representation of the gravity field variations in terms of mean fields calculated over a certain time span resulting in monthly (Watkins and Yuan, 2007; Bettadpur, 2007; Flechtner et al., 2010), 10-day (Bruinsma et al., 2010) or weekly (Flechtner et al., 2010) solutions. But there are various mass variation phenomena occurring on shorter time scales, for example the atmospheric variations or the barotropic motion of the ocean. In order to recover these fast gravity field variations as detailed as possible, it is reasonable to increase the temporal resolution with the goal of calculating daily GRACE solutions. This increase in temporal resolution, however, results in less observations per time span and therefore a reduced redundancy in the parameter estimation process. This leads to a decreasing accuracy of the estimated parameters with decreasing time span when the solutions are calculated individually. It can be

∗ Corresponding author. E-mail address: [email protected] (E. Kurtenbach). 0264-3707/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.jog.2012.02.006

assumed, however, that the gravity field does not change arbitrarily from one epoch to the next but that the information about the spatial and temporal correlation patterns can be approximated from geophysical models. Utilizing this knowledge, the temporal resolution can be enhanced without losing spatial information within the framework of a Kalman smoother estimation procedure. The Kalman smoother includes the stochastic prior information derived from geophysical models and the daily GRACE observations in a joint estimation process and delivers an updated state of the gravity field for each day. The stochastic information is introduced in terms of the process model which formulates a prediction of the current state resulting from the state of the previous time step. The process model is constructed from spatial and temporal covariance matrices derived from the output of the geophysical models. In Kurtenbach et al. (2009) a first attempt was carried out to derive daily gravity field solutions from GRACE level-1B data. In this publication the assumption was made that the gravity field does not change at all from one day to the next. The uncertainty of this assumption was modeled in terms of a covariance matrix derived only from the hydrological model WGHM (Döll et al., 2003; Hunger and Döll, 2008). In this first approach a stationary, homogeneous, and isotropic stochastic process on the sphere was assumed. In the investigations presented here, an improved version will be introduced which takes into account the full correlation structure

40

E. Kurtenbach et al. / Journal of Geodynamics 59–60 (2012) 39–48

between two subsequent days, represented by the spatio-temporal covariance matrices. Besides the hydrological model, atmospheric and oceanic variations will be taken into account as well to derive a more realistic model of the process dynamic. In Section 2 the Kalman smoother approach will be explained. In Section 3 a simulation scenario will be set up that uses geophysical models (different from those used to formulate the process model) to simulate daily GRACE observations for one year. This closed-loop scenario allows a variety of investigations to evaluate the potential of the method and to understand it in more detail. Especially the contribution of the GRACE observations to individual daily state updates will be investigated. In Section 4 the Kalman smoother approach will be applied to GRACE level-1B data resulting in daily GRACE solutions for the time span 2002–08 to 2009–08. These daily models are validated using individual data sets of ocean bottom pressure observations and vertical GPS station displacements. 2. The Kalman smoother approach The temporal variations of the gravity field of the Earth can be approximated by a finite spherical harmonics expansion with a time dependant set of potential coefficients: n  n+1 

GM  R R r N

V (, ϕ, r, t) =

n=0

anm (t) Ynm (, ϕ),

(1)

m=−n

with V(, ϕ, r, t) being the gravitational potential at a given location (, ϕ, r) and point in time t, GM representing the Earth’s gravitational constant, and R being the reference radius of the Earth. The Ynm are the spherical harmonic functions and anm (t) stand for the potential coefficients depending on the time t. In the following, these coefficients for each day will be estimated (Section 2.3) from a combination of daily GRACE observations (Section 2.1) and prior information about the temporal and spatial correlation patterns in terms of a process model which is derived from geophysical models (Section 2.2). 2.1. Observation model: GRACE level-1B data processing The standard processing strategy used at the University of Bonn is based on the functional model described in, for example, MayerGürr (2006) and Mayer-Gürr et al. (2007). It is based on a Fredholm integral equation of the first kind, which represents the solution of a boundary value problem to Newton’s equation of motion formulated for short arcs of the satellite orbit. The observation equations for the GRACE observations of one specific point in time t (i.e. a particular day) can be formulated as follows: lt + vt = At xt

with



C lt , lt



= Rt ,

(2)

with the GRACE observations lt , the residuals vt , the design matrix At , the covariance matrix of the observations Rt , and the unknown state vector xt = (anm )Tt , considered constant for one epoch (e.g. one day). From these assumptions it follows that the least-squares estimation procedure is solution of the normal equations: AT R−1 At xt = ATt R−1 t lt .

t  t 

 

Nt

nt

(3)

This reads explicitly xˆ t = N−1 t nt .

(4)

However, the stand-alone daily solution of this ill-posed problem is very inaccurate due to a low number of observations and the dominant observation noise. Fig. 1 shows the gravity field solution for

Fig. 1. One daily stand-alone GRACE solution (top) and its standard deviations in the spatial domain (bottom).

one arbitrary day using only one day of observations up to spherical harmonic degree N = 40 (top part). It becomes obvious that the solution bears little resemblance with the true gravity field features and that it is dominated by noise, which is confirmed by the very large standard deviation of the estimate shown in the bottom plot. 2.2. The process model To obtain a reasonable gravity field model from one day of GRACE data, additional information is required. Up to now, no temporal correlations have been introduced into the processing strategy. But it can be safely predicted that the estimate of the gravity field will not change arbitrarily from one day to the next because temporal gravity field variations are driven by geophysical processes which are not random. Therefore one can assume the simplified case that the time evolution of the gravity field can be written as the following first order Markov process: xt = Bxt−1 + w.

(5)

This represents a prediction of the gravity field coefficients from day t − 1 to the current day t. The prediction is characterized by the matrix of the process dynamic B. The normally distributed noise vector w∼N(0, Q) with the covariance matrix Q represents the accuracy of the prediction. In Kurtenbach et al. (2009) the simple assumption, that there are no changes between two states (B = I) with a covariance matrix Q derived from the hydrological model WGHM (Döll et al., 2003; Hunger and Döll, 2008), was introduced. Hereby a homogeneous and isotropic stochastic process on the sphere was assumed. In this section an improved version will be presented that takes into account the full correlation structure between two subsequent

E. Kurtenbach et al. / Journal of Geodynamics 59–60 (2012) 39–48

41

days, represented asymptotically in the stationary regime by the covariance matrix:



C

xt xt−1





 T 

=

 

,

(6)

with the auto-covariance matrix:  := C {xt , xt } ,

(7)

and the cross-covariance matrix:





 := C xt , xt−1 .

(8)

The matrix B is supposed to model the complete physics of the temporal behavior of the Earth’s gravity field. Such a geophysical model would be very complex and is therefore assumed to be unknown in the following. Instead, an estimate for B under consideration of Eqs. (7) and (8) is given by Moritz (1980) with ˆ =  −1 , B

(9)

which minimizes the prediction error w. Its covariance matrix then reads: −1

Q =  −  

T

 .

(10)

In practice, the matrices (7) and (8) of the Earth’s gravity field are unknown and therefore an empirical version,  and  , of these matrices has to be derived, which can be calculated using the output of geophysical models. These model outputs describe the mass distribution for each epoch t and are generally provided as point values on a global grid. In a first step these values are reduced by their mean, trend, and annual signal and then converted to a spherical harmonic expansion as used in Eq. (1). This results in a time series of T daily sets of spherical harmonic coefficients denoted by mt = (anm )Tt . These state vectors can then be applied to approximate the matrix  by the empirical auto-covariance matrix: 1 mi mTi , T T

≈=

(11)

i=1

and to derive the empirical cross-covariance matrix between two subsequent days: 1  mi mTi−1 , T −1 T

 ≈  =

(12)

i=2

as approximation of  . This way one can estimate the a-priori (i.e. before any satellite observations are used) dynamics of the Earth’s gravity field by xt =  

−1

xt−1 + w

(13)

with the prediction error: −1

w∼N(0,  −  

T  ).

Fig. 2. GRACE Kalman filter analysis approach: the system model (upper part) contains the process model and observation model which are merged in the two-step Kalman filter procedure (lower part), shown here in terms of normal equations as given in Eq. 3. +

−

updated state xˆ t−1 is predicted to the state xˆ t of the current epoch t: −

+

xˆ t = Bˆxt−1 + w,

(15)

using the pre-derived process model B. In many applications deterministic, physical knowledge about the dynamics of the process is introduced into the process model. In contrast to this the present investigation uses stochastic information only, which tends to be a much weaker type of information. After predicting the state, its covariance matrix will be error propagated to the next epoch: + T P− t = BPt−1 B + Q.

(16)

In the second filter step, new observations of epoch t are included and the prediction is improved by these observations. In terms of the representation of the observations introduced in Eq. (3), the update of the state reads: +

−

−

ˆ t ), xˆ t = xˆ t + P+ t (nt − Nt x (14)

In the next section the use of this prior model to assimilate additional observations using the Kalman filter technique is presented. 2.3. Merging observations and process model by Kalman filter

with the update of the covariance matrix of the state: − −1 P+ + Nt )−1 . t = ((Pt )

(18)

To obtain more insight into what is going on during the update procedure, Eq. (17) can be re-arranged to +

A common tool to estimate the state xt of a process given by (5) under the consideration of observations at each epoch t is the Kalman filter (see i.e. Gelb, 1974; Simon, 2006). In the following, the standard equations will be rewritten to fit the representation of the observations as introduced in Eq. (3). The standard Kalman filter approach is divided into two steps: prediction and update. In the prediction step the previously

(17)

−

− −1 ˆ t + P+ ˆt xˆ t = P+ t Nt x t (Pt ) x

(19)

with the GRACE-only solution xˆ t from Eq. (4). Here one can see that the Kalman filter weights the two types of pseudo-observations xˆ t − + and xˆ t in a least-squares sense to derive the updated state xˆ t . The weight matrix for the GRACE-only solution is then given by −1 WtGRACE := ((P− + Nt )−1 Nt . t )

(20)

42

E. Kurtenbach et al. / Journal of Geodynamics 59–60 (2012) 39–48

Fig. 3. Differences in terms of equivalent water heights between simulated refer− + ence signal xAOH and prediction xˆ t (top, RMS: 2.1 cm) and update xˆ t (bottom, RMS: 1.3 cm) for one day.

Fig. 4. Standard deviations in terms of equivalent water heights [cm] of prediction − + xˆ t (top) and update xˆ t (bottom) of the states used in Fig. 3.

3. Simulation study The main diagonal of this weight matrix will be used later on to approximate the contribution of one day of GRACE observations to the updated state vector. A scheme of the approach is shown in Fig. 2. The system model is composed of the process model (Section 2.2), derived from stochastic prior information from geophysical model output, and the observation model described in Section 2.1. In the lower part of Fig. 2 the two step Kalman filter approach (prediction and update) is illustrated.

2.4. The Kalman smoother In the Kalman filter approach as described above, the estimated state (gravity field) at the current time t only depends on the earlier states, which results in the method being suitable for real-time applications. But as today, GRACE processing is not required to provide (near-)real-time results, we suggest to include all available information into the determination process. Therefore it is reasonable to implement a Kalman smoother instead of a Kalman filter procedure. This can be achieved by simple forward and backward filtering, and subsequent calculation of the weighted minimum-variance mean of the two solutions. It was done for the real GRACE data analysis of Section 4 in terms of the computationally more efficient RTS smoother (named after the authors of Rauch et al., 1965), for more details see Simon (2006). In the simulation study, only the (forward) Kalman filter was used to investigate the performance of the process model.

In order to evaluate the potential of the Kalman filter method and to understand the approach in more detail, a simulation scenario was set up first to test the procedure in a closed-loop environment. 3.1. Setup Mass variations for the year 2006 were simulated from geophysical models. The Global Land Data Assimilation System (GLDAS, Rodell et al., 2004) was used to obtain the hydrological variations (total water storage), the MOG2D ocean model (Carrère and Lyard, 2003) and the NCEP atmospheric model (Kalnay et al., 1996) were evaluated for deriving the oceanic and atmospheric pressure

+

Fig. 5. Contribution of GRACE observations to estimated state xˆ t in the spectral domain (in percent).

E. Kurtenbach et al. / Journal of Geodynamics 59–60 (2012) 39–48

43

measurements, and star camera observations. Each of the data sets was generated with a sampling rate of 5 s and corrupted with white noise (GPS  = 2 cm and K-Band range-rates  = 0.2 ␮m/s). As described in Section 2.2, implementing the GRACE Kalman filter approach requires knowledge of the spatial and temporal correlation patterns of the expected gravity field signal for building the process model. This knowledge is represented by the empirical covariance matrices  and  , which are derived from the output of geophysical models. Because in reality the models provide only an approximation of the true mass variations, the correlation patterns are intentionally derived from different models than those that were used for generating the simulated reference signal. In this way, too strong correlations between the geophysical a-priori information (in terms of the covariance matrices) and the input signal will be avoided. The models used to build the process model are the WaterGAP Hydrology Model (WGHM, Döll et al., 2003; Hunger and Döll, 2008), the ocean circulation model OMCT (Dobslaw and Thomas, 2007), and ECMWF atmospheric model data (Research Department ECMWF, 2008). Using this input information, a series of daily spherical harmonic coefficients for degrees n = 2 . . . 40 was derived applying the Kalman filter process as described in Section 2.3. Some representative results of the simulation study will be described below.

Fig. 6. Contribution of GRACE observations from Fig. 5 propagated to the spatial domain: strong correspondence to GRACE ground tracks (in red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

variations. From this signal, a long-term mean, an annual and semi-annual signal and a linear trend were reduced. This results in a reference signal, in the following referred to as xAOH , which contains only the short-term variations with daily resolution. The static gravity field model ITG-Grace03s was added to the temporal variations, after converting them from mass to potential, and daily GRACE observations were simulated for one year. The simulated data set comprises GRACE orbits, K-band ranges, accelerometer

3.2. Simulation results The simulation study provides as the state vector a set of spherical harmonic coefficients for each day, with the result being a

−

+

Fig. 7. Left: Time series of simulated reference signal xAOH (black), predicted state xˆ t (blue) and updated state xˆ t (red) for two different points. Right: zoom-in for three months. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

44

E. Kurtenbach et al. / Journal of Geodynamics 59–60 (2012) 39–48

−

+

Fig. 8. Temporal characteristics of predicted signal xˆ t (left) and updated signal xˆ t (right): (a) correlation coefficient, (b) error RMS (in cm water heights), both with respect to simulated reference signal xAOH , and (c) amount of reference signal that can be explained by predicted and updated signal (in percent), respectively.

two-dimensional map of the gravity field for each point in time. As these three dimensions cannot easily be illustrated in 2D figures, the result will be shown first for only one time epoch. In a second step the complete time series will be illustrated for one grid point. Finally, the time series will be investigated for the entire map. 3.2.1. Investigating one arbitrary day The Kalman filter output for one arbitrary day of the simulation scenario is analyzed in Fig. 3 by calculating the difference of the solution compared to the reference input signal for each point in space. The left graph shows the differences for just the predicted signal − xˆ t (i.e. without adding the GRACE observations of the specific day), whereas in the right part of the figure the differences are plotted + for the update xˆ t (i.e. including observations of the specific day). Adding the GRACE observations of just one single day significantly reduces the differences, the RMS decreases from 2.1 cm to 1.3 cm. Fig. 4 shows the standard deviations (formal errors) of the states − + xˆ t (top) and xˆ t (bottom), calculated from the covariance matrices given in Eqs. (16) and (18), respectively, both propagated to the spatial domain. Here it can be observed that the standard deviations of the update are significantly reduced compared to those of the prediction. Such an improvement caused by introducing the GRACE observations is remarkable keeping in mind that the GRACE

observations by themselves have the very large standard deviations illustrated in the lower part of Fig. 1. It is consistent with what can be seen in Fig. 3. To further investigate the contribution of GRACE, the main diagonal of Eq. 20 is used to approximate the contribution of one day of GRACE observations on the updated state vector. A similar interpretation of Eq. 20 based on statistical redundancies in the estimation procedure can be found in Bouman (2000). The contribution of one exemplary day to the estimated state + vector xˆ t of the respective day is illustrated in Fig. 5 for each of the spherical harmonic coefficients. It becomes clear that due to the low data coverage provided by only one day of satellite data the contribution is larger in the lower spherical harmonic coefficients (upper middle part of the triangle). This implies the obvious conclusion that longer wavelengths can more easily be recovered from few observations. It can also be observed that above an order of approximately m = 15 there is no more contribution from GRACE left. This can be understood from the orbit configuration with 15 revolutions of the satellite per day resulting in 30 crossings of the equator. Fig. 6 shows the contributions of the GRACE observations of the same day propagated from the frequency domain to the space domain. Here it can be noticed that the contribution is not equally distributed around the Earth. In regions towards the poles the

E. Kurtenbach et al. / Journal of Geodynamics 59–60 (2012) 39–48

45

Fig. 9. Left: OBP recorder observations in the northern Atlantic ocean (black), daily GRACE solutions (red) and AOD1B model (green). Right: zoom-in on a smaller time span. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

observation density is much higher and – due to the polar orbits – the contribution is much stronger compared to the regions around the equator. Furthermore a strong correlation between the contributions and the satellite ground tracks (plotted as red lines) can be observed. GRACE contributes only in those areas were it has collected information on the specific day.

3.2.2. Time series for individual grid points From the daily solutions, time series of mass variations (in terms of equivalent water heights) can be generated for each point in space. Fig. 7 illustrates two such time series, one for a point at higher latitude, where the GRACE observations are dense and one for a point at the equator, where the observation density is lower. In each of the plots the black line represents the reference signal and the red line the outcome of the Kalman filter procedure, i.e. the update for each day. The blue line shows the temporal evolution when for each − day only the prediction step xˆ t is plotted. It can be observed that for the higher latitudes (Fig. 7a) the red curve follows the input signal quite well and close as compared to the blue line. This illustrates the impact of the GRACE observations from only one day of data. The correlation coefficient improves from 0.73 between reference signal and prediction to 0.95 between reference signal and updated signal. Also the error RMS of 7.8 cm between reference signal and prediction improves to 4.3 cm for the updated signal. For the point at the equator (Fig. 7b), the signal itself is much smaller and thus more difficult to be measured by GRACE. In combination with the less dense satellite observations this leads to the red curve fitting

the black curve not as well as it is the case for a point at higher latitudes. Also the impact of the observations (difference between red and blue curve) is less prominent. The correlation coefficient and the error RMS improve for the updated signal from 0.67 and 4.9 cm to 0.75 and 4.3 cm, respectively.

3.2.3. Time series for the entire grid Such calculations can be carried out for each point in space, as shown in Fig. 8a where the correlation coefficients are plotted globally. The left part shows the correlation of the reference signal with the prediction, in the right part the daily observations are included in the estimate as well. Obviously the correlations for the update are significantly larger than for the prediction, which illustrates the influence of the daily observations. As the correlation coefficients are not sensitive to differences in amplitude between the two analyzed signals, the error RMS with respect to the reference signal was investigated as well. It provides a measure for the absolute differences between the reference signal and the Kalman filter results. These values are presented in Fig. 8b, where the error RMS is plotted for the prediction (left) and the update (right). Again, the update performs better with smaller RMS values. It can be observed, however, that the RMS depends on the magnitude of the signal, as it is generally larger in areas with large mass variations, such as in the northern latitudes (with large atmospheric variations) and the Antarctic Circumpolar Current. Therefore as a third measure Fig. 8c illustrates the signal reduction which explains the ratio of error RMS and signal RMS. It can be interpreted as information about

46

E. Kurtenbach et al. / Journal of Geodynamics 59–60 (2012) 39–48

Fig. 10. Comparison daily GRACE solutions with radial GPS station movements. Large figure: amount of the GPS signal that can be explained by GRACE. Top figures: Comparison of time series for individual stations for GPS (black), GRACE (red) and AOD1B (green), left: GLPS (Puerto Ayora), right: NRIL (Norilsk). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

how many percent of the reference signal can be explained by the predicted and updated signal, respectively.

information content of the daily GRACE models, they were compared to independent data sets as described below.

4. Real GRACE L1B data analysis

4.1. Comparison with ocean bottom pressure (OBP) data

The Kalman smoother approach as described in Section 2 was applied to calculate daily gravity field solutions from GRACE level1B data, covering the time span 2002–08 to 2009–08. For each day a set of spherical harmonic coefficients was derived for degrees n = 2 . . . 40 as differences to the static gravity field ITG-Grace2010s, using the standard background models described in Mayer-Gürr et al. (2010). For the generation of the process model the same geophysical models were used as described in the simulation scenario in Section 3, that is total water storage and pressure output from the models WGHM, OMCT, and ECMWF model. In order to guarantee that the GRACE solutions are not biased towards the model values themselves but that only the stochastic behavior is exploited, the model output from the years 1976 to 2000 (i.e. outside the GRACE time span) was used. The daily solutions described by the present paper are part of the GRACE gravity field model ITG-Grace2010 (Mayer-Gürr et al., 2010) and can be downloaded at http://www.igg.unibonn.de/apmg/index.php?id=itg-grace2010. To evaluate the

Ocean bottom pressure recorders represent an independent type of measurement for the integrated mass column of oceanic and atmospheric variations. Therefore, they can be used as a validation tool for the mass variations derived from GRACE. In Fig. 9 the comparison was performed for two OBP recorder representing the conditions in different latitudinal regions. One is located in the northern Atlantic (ABPR 1) and one given at a location in the equatorial Atlantic region (P128 2). The black lines show the temporal evolution of the OBP data as provided by Macrander et al. (2010). These in situ measurements are compared to the atmosphere and ocean de-aliasing product (AOD1B, green line) used in the GRACE L1B data analysis (Flechtner, 2007) and to the daily GRACE models obtained by the Kalman smoother approach (red line) described here. The AOD1B model represents the knowledge of the temporal high-frequent mass variations on a global scale without the use of GRACE. Again it can be seen that for the higher latitude (Fig. 9a) the daily GRACE solutions show better agreement to the independent OBP data set than the modeled atmospheric

E. Kurtenbach et al. / Journal of Geodynamics 59–60 (2012) 39–48

47

Fig. 11. Improvement of correlation coefficients of AOD time series against GPS time series after adding daily GRACE solutions.

and oceanic variations of AOD1B. In the equatorial region the correspondence between GRACE and OBP is not good due to the smaller signal and the sparse data coverage. In Kurtenbach (2011) the same investigation was carried out for a larger number of OBP stations world wide. While for the majority of the stations the correlation of GRACE with OBP improves compared to AOD, this is not the case for all of the stations. Here it also has to be kept in mind that OBP measurements represent a very local type of measurement containing very local gravity field information and therefore the significance for GRACE validation might be limited. A more detailed error assessment including an evaluation of the performance of the OBP recorders will have to be subject to further investigations.

4.2. Comparison with GPS station movements Mass variations at the Earth’s surface result in geometrical deformations of the Earth’s crust which can be measured by GPS receivers. Therefore, the global network of permanent GPS stations delivers a further set of independent observations which can be used as comparison for GRACE gravity field models. The vertical station displacements of the reprocessed time series of the International GPS Service (IGS), see Tesmer et al. (2011), were compared on a daily basis to the GRACE Kalman solutions transformed to vertical loading using the load Love numbers of Gegout (2005). The results are displayed in Fig. 10. The large figure shows the percentage of the signal of each of the stations that can be explained by the GRACE observations, i.e. it illustrates the reduction of the RMS values of the GPS observations when the GRACE solutions are subtracted. The smaller figures on top show the time series for two GPS stations. The in situ GPS observations are given by the black curve, the GRACE time series by the red curve and as a comparison the AOD1B model is plotted in the green line as well. The right figure represents an example for a station with very good agreement (NRIL, Norilsk). It is situated in northern Siberia, where the data coverage by GRACE is very good and the daily signal of the radial displacements is large due to strong atmospheric mass variations. In contrast to this, the correlation of the time series is very poor in the left figure (both for the Kalman solutions and for the AOD1B product) which shows a station near the equator in the Pacific ocean (GLPS, Puerto Ayora). Here the signal itself is very small and the data coverage by GRACE is less dense. To give a further idea of the global results, Fig. 11 shows the improvement of the correlation coefficient between GPS and GRACE compared to the respective correlation between GPS and AOD for the different IGS stations. There are only a few stations where the correlation decreases when adding the daily GRACE solutions, while most of the stations experience an improvement of the correlation.

5. Conclusions and outlook An improved approach to calculating daily GRACE gravity field models in the framework of a Kalman filter and smoother estimation procedure has been introduced. The daily solutions are not derived individually, but the full spatial-temporal correlation patterns obtained from the output of geophysical models is introduced as additional information. The approach was tested in a closed-loop simulation scenario and then applied to GRACE level-1b data analysis resulting in daily solutions as part of the new GRACE gravity field model ITG-Grace2010. Finally, the results were evaluated using independent data sets, i.e. ocean bottom pressure data and vertical GPS station displacements. The investigation presented here showed that the contribution of the GRACE data to the daily solutions is particularly large in areas where the observation density is high at each specific day. This is specifically the case along the satellite’s orbit and in the polar regions. Independent comparisons have shown that in particular in those regions the daily solutions provide additional temporal gravity field information and represent an improvement to existing models as shown by the comparison to the AOD1B atmosphere and ocean mass variations. Thus it can be concluded that Kalman filter and smoother procedure represents a very effective way to determine the temporal high-resolution gravity field variations. Further validation of the results with external data sets are in progress, one example comparing the daily GRACE solutions to altimetry measurements can be found in Bonin and Chambers (2011). Since these daily variations cannot be modeled by longer temporal means and are not ideally represented in the applied AOD1B de-aliasing product, they will cause aliasing errors in monthly and even weekly gravity field solutions. As mentioned above, the daily Kalman solutions capture more information about these mass variations. Therefore, they can be introduced as an improved dealiasing product; a strategy which was carried out very successfully in the calculation of the monthly and static ITG-Grace2010 solutions. As mentioned above, the contribution of the daily GRACE observations is not uniform around the Earth but large in those areas with a denser data coverage. This suggests to explore the advantages of the approach even more when the characteristics of the gravity field representation can be tailored to the given data set. This will be particularly straightforward using space-localizing basis functions instead of the globally defined spherical harmonics. A first implementation of the GRACE Kalman Filter approach using the localizing radial basis functions introduced by Eicker (2008) has successfully been carried out and will be further improved in the future.

48

E. Kurtenbach et al. / Journal of Geodynamics 59–60 (2012) 39–48

Acknowledgement We thank the Deutsche Forschungsgemeinschaft (DFG) for the funding of this work within the framework of the priority programme SPP1257 “Mass transport and mass distribution in the system Earth”. References Bettadpur, S., 2007. Utcsr level-2 processing standards document for level-2 product release 0004. CSR Publ. GR-03-03. Bonin, J., Chambers, D., 2011. Evaluation of high-frequency oceanographic signal in grace data: implications for de-aliasing. Geophys. Res. Lett., 38. Bouman, J., 2000. Quality Assessment of Satellite-based Global Gravity Field Models. Ph.D. thesis. Publications on Geodesy 48, Netherlands Geodetic Commission. Bruinsma, S., Lemoine, J.M., Biancale, R., Valés, N., 2010. Cnes/grgs 10-day gravity field models (release 2) and their evaluation. Adv. Space Res. 45, 587–601. Carrère, L., Lyard, F., 2003. Modeling the barotropic response of the global ocean to atmospheric wind and pressure forcing—comparisons with observations. Geophys. Res. Lett., 30. Dobslaw, H., Thomas, M., 2007. Simulation and observation of global ocean mass anomalies. J. Geophys. Res., 112. Döll, P., Kaspar, F., Kaspar, B., 2003. A global hydrological model for deriving water availability indicators: model tuning and validation. J. Hydrol. 270, 105–134. Eicker, A., 2008. Gravity Field Refinement by Radial Basis Functions from In-situ Satellite Data. Ph.D. thesis. University of Bonn. Flechtner, F., 2007. AOD1B Product Description Document for Product Releases 01 to 04. Technical Report. Geoforschungszentrum Potsdam. Flechtner, F., Dahle, C., Neumayer, K.H., Koenig, R., Foerste, C., 2010. The release 04 champ and grace Eigen gravity field models. In: Gruber, F., Guentner, T., Mandea, A., Rothacher, M., Wickert, M.J. (Eds.), Satellite Geodesy and Earth System Science Observation of the Earth from Space. Springer. Gegout, P., 2005. Load love numbers. http://gemini.gsfc.nasa.gov/aplo/Load Love2 CM.dat. Gelb, A., 1974. Applied Optimal Estimation. The MIT Press, Cambridge, Massachusetts and London, England. Hunger, M., Döll, P., 2008. Value of river discharge data for global-scale hydrological modeling. Hydrol. Earth Syst. Sci. 12, 841–861. http://www.hydrol-earth-systsci.net/12/841/2008/.

Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L., Iredell, M., Saha, S., White, G., Woollen, J., Zhu, Y., Chelliah, M., Ebisuzaki, W., Higgins, W., Janowiak, J., Mo, K.C., Ropelewski, C., Wang, J., Leetmaa, A., Reynolds, R., Jenne, R., Joseph, D., 1996. The ncep/ncar 40-year reanalysis project. Bull. Am. Meteorol. Soc. 77, 437–471. Kurtenbach, E., 2011. Entwicklung eines Kalman-Filters zur Bestimmung kurzzeitiger Variationen des Erdschwerefeldes aus Daten der Satellitenmission GRACE. Ph.D. thesis. University of Bonn. Kurtenbach, E., Mayer-Gürr, T., Eicker, A., 2009. Deriving daily snapshots of the earth’s gravity field from grace l1b data using kalman filtering. Geophys. Res. Lett., 36. Macrander, A., Böning, C., Boebel, O., Schröter, J., 2010. Validation of grace gravity fields by in-situ data of ocean bottom pressure. 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. Springer. Mayer-Gürr, T., 2006. Gravitationsfeldbestimmung aus der Analyse kurzer Bahnbögen am Beispiel der Satellitenmissionen CHAMP und GRACE. Ph.D. thesis. University of Bonn, Germany. Mayer-Gürr, T., Eicker, A., Ilk, K.H., 2007. Itg-grace02s: a grace gravity field derived from range measurements of short arcs. In: Proceedings of the 1st International Symposium of the International Gravity Field Service “Gravity Field of the Earth”, Istanbul, pp. 193–198. Mayer-Gürr, T., Kurtenbach, E., Eicker, A., 2010. Itg-grace2010 gravity field model. http://www.igg.uni-bonn.de/apmg/index.php?id=itg-grace2010. Moritz, H., 1980. Advanced Physical Geodesy. Wichmann, Karlsruhe. Rauch, H., Tung, F., Striebel, C., 1965. Maximum likelihood estimates of linear dynamic systems. AIAA J. 3, 1445–1450. Research Department ECMWF, 2008. Ifs documentation. Part iii. Dynamics and numerical procedures (cy33r1). http://www.ecmwf.int/research/ifsdocs. Rodell, M., Houser, P., Jambor, U., Gottschalck, J., Mitchell, K., Meng, C.J., Arsenault, K., Cosgrove, B., Radakovich, J., Bosilovich, M., Entin, J., Walker, J., Lohmann, D., Toll, D., 2004. The global land data assimilation system. Bull. Am. Meteorol. Soc. 85, 381–394. Simon, D., 2006. Optimal State Estimation: Kalman, H Infinity, and Nonlinear Approaches. Wiley-Interscience. Tapley, B., Bettadpur, S., Watkins, M., Reigber, C., 2004. The gravity recovery and climate experiment: mission overview and early results. Geophys. Res. Lett., 31. Tesmer, V., Steigenberger, P., van Dam, T., Mayer-Gürr, T., 2011. Vertical deformations from homogeneously processed grace and global GPS long-term series. J. Geodesy. Watkins, M., Yuan, D.N., 2007. Jpl level-2 processing standards document for level-2 product release 04. ftp://podaac.jpl.nasa.gov/pub/grace/doc/.