Simulation of the time-variable gravity ﬁeld by means of coupled geophysical models

. Time variable gravity ﬁelds, reﬂecting variations of mass distribution in the system Earth is one of the key parameters to understand the changing Earth. Mass variations are caused either by redistribution of mass in, on or above the Earth’s surface or by geophysical processes in the Earth’s interior. The ﬁrst set of observations of monthly variations of the Earth gravity ﬁeld was provided by the US / German GRACE satellite mission beginning in 2002. This mission is still providing valuable information to the science community. However, as GRACE has outlived its expected lifetime, the geoscience community is currently seeking suc-cessor missions in order to maintain the long time series of climate change that was begun by GRACE. Several studies on science requirements and technical feasibility have been conducted in the recent years. These studies required a realistic model of the time variable gravity ﬁeld in order to perform simulation studies on sensitivity of satellites and their instrumentation. This was the primary reason for the European Space Agency (ESA) to initiate a study on “Monitoring and Modelling individual Sources of Mass Distribution and Transport in the Earth System by Means of Satellites”. The goal of this interdisciplinary study was to create as realistic as possible simulated time variable gravity ﬁelds based on coupled geophysical models, which could be used in the simulation processes in a controlled environment. For this purpose global atmosphere, ocean, continental hydrology and ice models were used. The coupling was performed by using consistent forcing throughout the models and by including water ﬂow between the di ﬀ erent domains of the Earth system. In addition gravity ﬁeld changes due to solid Earth processes like continuous glacial isostatic adjustment (GIA) and a sudden earthquake with co-seismic and post-seismic signals were modelled. All individual model results were combined and converted to gravity ﬁeld spherical harmonic series, which is the quantity commonly used to describe the Earth’s global gravity ﬁeld. The result of this study is a twelve-year time-series of 6-hourly time variable gravity ﬁeld spherical harmonics up to degree and order 180 corresponding to a global spatial resolution of 1 degree in latitude and longitude. In this paper, we outline the input data sets and the process of combining these data sets into a coherent model of temporal gravity ﬁeld changes. The resulting time series was used in some follow-on studies and is available to anybody interested.

. Time variable gravity field sources in terms of spatial resolution and time variability including GRACE sensitivity and goals for future gravity field missions. the solid Earth and environmental surface processes within the Earth's outer fluid envelop (excluding tides). It shall be noted that "realistic" is understood in the sense that all spectra of variability shall be included in a comprehensive timevariable gravity field model in order to perform simulations approximating reality. A simulation time period of 12 yr spanning the years 1995 to 2006 was chosen. The resulting predictive model of the time variable gravity field is unique in the sense of signal sources incorporated, coupling of geophysical models and a decadal simulation period. No similar data sets have been published yet in the geoscience community. Here, we address the development of the realistic Earth model, which is described in detail below. The second step was to use the realistic Earth mass model to develop satellite orbits and configurations to monitor the changes such that we are able to separate the different components of the mass variability. The twin-satellite GRACE mission has shown that time variable gravity field is observable from space on a monthly base with spatial resolutions of a few hundred kilometres (Tapley et al., 2004). GRACE observes inter-satellite distance changes with a few µm accuracy, which are caused to a large extent by irregularities of the mass distribution in, on or above the Earth. Studies on future gravity field mission constellations target on higher temporal and spatial resolution by enhancing the observation accuracy and/or by implementing new satellite configurations. The subsequent research in which different satellite orbits and configurations were tested has already been described in Visser et al. (2010).
An overview of time variable gravity field sources in terms of spatial resolution and time variability is shown in Fig. 1. The horizontal axis defines the spatial resolution, while the vertical axis specifies the time variability from instantaneous to static. The bubbles define signals from various sources and their estimated space-time variability. The colours of the bubbles classify the sources for different domains. In addition the box identifies, what signals are sensitive to the GRACE mission (see below). The two dotted lines indicate somehow the signals to be observable by future gravity field missions, without specifying the requirements in detail.
Estimates of individual mass variations are provided via the different geophysical models for each component of the Earth mass system, namely the atmosphere, the continental hydrology, the continental ice sheets, the ocean and the solid Earth. As gravity is capable of observing only the integrated effect of the mass distribution at a specific time, simulating the time variable gravity field from the separate geophysical models needs to be made consistent and the models need to be coupled. Consistency means that geophysical model results need to be converted, such that they provide the same quantity for mass distribution at a specific time (e.g. pressure, equivalent water height or gravitational attraction), that the geospatial sampling is identical (e.g. geographical grid) and that they all refer to the same reference values (e.g. a mean field). The procedure applied to secure consistency is described in some detail in Sect. 7. More important, however, is the coupling of the geophysical models that together define the water cycle. For a realistic Earth model, we need to insure that mass flowing from one component to the other domain is not "lost" in the system or conversely that it is not counted twice. The coupling was performed by imposing that the forcing data be consistent over all models, i.e. all our surface mass models use forcing parameters from the same source (in our case the ECMWF atmospheric model) (see Sect. 2) and that water flowing from the continents (including continental ice) to the oceans is taken into consideration by the models. By using an identical forcing model, we guarantee that water exchange between the atmosphere and the Earth surface (continents and oceans) is taken into account intrinsically.
For coupling the geophysical models on the ground, specific adaptations were made, including the harmonization of land/ocean masks and the provision of outgoing water mass from one domain to another. Detailed descriptions of the geophysical models and the adaptations made for coupling them are provided in Sects. 3-5 for the continental hydrology, the continental ice and the ocean, respectively.
For the solid Earth processes, global isostatic adjustment was included as a continuous process and an earthquake as a sudden event. Their predicted impact on the gravity field was taken into account. Details on this aspect of the modelling will be provided in Sect. 6. Finally, Sect. 8 summarizes the work performed and provides some conclusions to be drawn from the results obtained.

Atmosphere
The atmosphere is one of the main contributors to mass variations in the Earth system. Atmospheric mass variations are present for periods ranging from a few hours to decadal. The primary manifestation of atmospheric mass variations, are the distribution and propagation of the pressure systems around the globe. However high frequency mass variability also occurs in the form of diurnal and the semi-diurnal atmospheric tides.
For modelling the atmospheric mass variations as well as for forcing the geophysical models, a global atmospheric model was applied. For this study we used the operational, the re-analysis (ERA-40) and the related short-term forecast models driven by the European Centre for Medium-Range Weather Forecasts (ECMWF). The latter model needs to be used for some parameters not available in the other two model versions. In order to identify the impact of using atmospheric data with different consistency over time, we used atmospheric parameters from the operational and the reanalysis model for different periods. Reanalysis data (ERA-40) were used for the years 1995 to 1999, while for the time span from 2000 until 2006 data from the ECMWF operational analysis were used. The list of atmospheric parameters and the geophysical models for which they are subsequently used is provided in Table 1. In general, we used the 6-hourly time series of the atmospheric parameters indicated. All atmospheric parameters were extracted from the ECMWF data base according to the requirements of the geophysical models. More details about the hydrology, ocean and continental ice models are provided in the subsequent Sects. 3 to 5.
For the atmospheric mass distribution per time step, the required parameters were extracted as 1 × 1 • equi-angular grid files in the GRIB format (standard format used within the atmospheric community). The conversion from the original representation as stored in the ECMWF archive into equiangular grids is performed directly during data base extraction process. No further interpolation or conversion is required. For this simulation, we used the simplified surface pressure approach disregarding the vertical structure of the atmosphere (cf. Flechtner, 2007). For the optimizing satellite orbits and configurations one would need in addition the geopotential surface and the multi-layer atmospheric temperature and specific humidity data. For the combination of the different mass fields as described in Sect. 7 we finally have available 6 hourly surface pressure grids in units of Pascal from the ERA-40 reanalysis model (for the years 1995 to 1999) and from the operational analysis (for years 2000 to 2006).

Continental hydrology
We used the global hydrological model PCR-GLOBWB (Van Beek and Bierkens, 2009;Bierkens and Van Beek, 2009) to estimate terrestrial water storage anomalies. PCR-GLOBWB calculates for each grid cell (0.5 × 0.5 • globally) and for each time step (daily) the water storage in two vertically stacked soil layers and an underlying groundwater layer, as well as the water exchange between the layers and between the top layer and the atmosphere (rainfall, evaporation and snow melt). The model also calculates canopy interception and snow storage. Sub-grid variability is taken into account by considering separately tall and short vegetation, open water, different soil types and the area fraction of saturated soil and the frequency distribution of groundwater depth based on the surface elevations of the 1 × 1 km Hydro1k data set. Fluxes between the lower soil reservoir and the groundwater reservoir are mostly downward, except for areas with shallow groundwater tables, where fluxes from the groundwater reservoir to the soil reservoirs are possible (i.e. capillary rise) during periods of low soil moisture content. The total specific runoff of a cell consists of saturation excess surface runoff, melt water that does not infiltrate, runoff from the second soil reservoir (interflow) and groundwater runoff (baseflow) from the lowest reservoir. To calculate river discharge, specific runoff is accumulated along the drainage network by means of kinematic wave routing including storage effects and evaporative losses from lakes, reservoirs and wetlands.

Adaptation of model made for this study
Forward modelling with PCR-GLOBWB returned daily fields of terrestrial freshwater storage with a resolution of 0.5 • for the period 1957-2006. The model was forced with ECMWF meteorological data for precipitation, air temperature and actual evaporation. For the period 1957-1999 inclusive, daily surface fields were obtained from ERA-40 reanalysis forecasts, from 2000 onwards from forecasts of the Operational Archive. To preserve consistency in the meteorological forcing, actual evaporation from the ECMWF forecasts was imposed directly on the model and partitioned over the different surface areas within each cell on the basis of cover fraction. Following the parameterization of the model, bare soil evaporation was drawn from the upper soil layer whereas transpiration by vegetation was broken down over both soil layers on the basis of root fractions. The imposed evaporative flux was not reduced, unless the absolute soil water storage was insufficient to sustain it.
Since glacier melt has only limited influence on discharge, glacier dynamics were insufficiently represented in PCR-GLOBWB to resolve variations in water storage over this type of terrain realistically. To amend this, four regions with significant terrestrial glacier mass were included separately in the model in addition to the land ice masses of Greenland and Antarctica, which were modelled separately by BGC (see below). The melt of the land ice mass of Greenland and Antarctica was supposed to drain directly to the ocean and the corresponding areas masked out. For the four glaciated regions (Alaska, Alps, Karakoram, and Himalayas), melt contributes to streamflow in PCR-GLOBWB. Glacier melt was assumed to follow a secular trend of 1 % per year from 1990 onwards starting with initial volumes of 7500, 1500, 4000 and 4000 km 3 for the respective areas (J. Bamber, personal communication, 2008). To obtain a finer representation of glaciers in these regions, the glacier volumes were downscaled proportional to the glacier coverage within each 0.5 • cell as represented by the FAO digital soil map of the world (FAO, 2011). At the beginning of 1990, the glacier volumes in these four regions were initialized at the estimated volumes, with variations arising due to precipitation and evaporation over time, and the secular melt rate applied. This melt rate was further broken down into monthly values on the basis of the temperature climatology of the ERA-40 dataset with the melt rate being proportional to the deviation above the mean annual temperature. Prior to this date, the melt rate was assumed to be zero for these four regions. An additional problem was encountered in those areas where the long-term temperature of the EMCWF data did not exceed 0 • C and that were not identified as glaciated areas a priori. The corresponding cells were identified entirely as glaciers and the mass balance of precipitation, evaporation and glacier melt evaluated with the latter amounting to 100 % of the accumulated volume per year.
Outside of the glaciated areas, stores in PCR-GLOBWB were initialized by a spin-up run in which the ERA-40 climatology was repeated until a dynamic steady-state in discharge was obtained. Anomalies include variations in the following stores of the terrestrial part of the hydrological cycle: open freshwater bodies (river stretches and active volumes of lakes and reservoirs), active groundwater storage, snow cover, interception storage, and soil moisture. Continental water stores, which were not considered, include ground ice, fossil groundwater bodies and inactive lake volumes and salt-water intrusions into estuaries or permeable coastal reservoirs or any anthropogenic abstractions or diversions of freshwater.

Description of resulting hydrology mass fields
Including the glaciated areas, anomalies in the hydrology mass fields represent the seasonal and inter-annual variations realistically over large areas. Small inconsistencies arise due to the fact that the spatial resolution of the ECMWF forcing is finer for the Operational Archive than for the ERA-40 reanalysis; at these transitions abrupt changes in the topographic effect on temperature leads to sudden snow melt and variations in the available storage. Overall, however, the quality of the hydrology mass fields is good over the Northern Hemisphere where the forecasts are well constrained by meteorological information, the amount of precipitation is large, and the resulting conversion into runoff is not unduly influenced by model errors. Elsewhere the quality can be affected by a poorer meteorological forcing and by biased and uncertain variables and unresolved processes (e.g. simple representation of snow and glacier melt). In addition, the ERA-40 reanalysis demonstrates variable skill in predicting precipitation for different zones and for different periods (Troccoli and Kållberg, 2004) and its water balance not closed. Particularly troubling are the long-term trends over the 14 % of the terrestrial surface that has no direct link to the oceans. Such endorheic basins are particularly sensitive to small variations in evaporative losses (and errors therein) as the accumulated runoff can only be lost to the atmosphere. Thus, gains or losses over these areas will directly affect the simulated storage (Fig. 2).
Negative fluxes over North America and Europe are relatively small when the corresponding area is taken into consideration. They are more substantial for the arid continents of Australia and Africa over the simulation period but in these water-limited environments this will not have large repercussions on the total storage. Where the trends are positive, for example for the Caspian Sea and Asia, the accumulations are significant as these areas correspond to 10 % of the total terrestrial land mass. For the Caspian Sea, the increases from the mid-1990s onwards can be linked to increased discharge from the Volga due to a stronger North Atlantic Oscillation (NAO) and should not be discredited as a mere model artefact.
For South America, however, the strong positive trend is erroneous. In this study, such spurious trends in water storage have been removed from the hydrology mass storage fields as a post-processing step. Yet, drawing from the experience of this study it would be advisable to force the model with a more consistent meteorological dataset and to impose reference potential evaporation. For each surface within PCR-GLOBWB, the specific potential evaporation can then be calculated using correction factors (or crop factors) (Allen et al., 1998) and scaled to the actual rate on the basis of water availability within the hydrological model. Tuning of the correction factors can then reduce spurious long-term trends. Although the a posteriori removal of trends largely resorts in the same effect, tuning the correction factors is slightly more transparent and may preserve regional and temporal detail that is currently lost.

Background
The purpose of the ice mass flux models was not to, necessarily, reproduce the observed mass trends but to simulate realistic seasonal and secular behaviour. The primary focus was on the contribution of the two largest ice masses: the Greenland and Antarctic ice sheets (GrIS, AIS) in terms of forward modelling. For four glacierized regions, the European Alps, Himalaya, Karakoram and Gulf of Alaska ( Fig. 4), we prescribed plausible secular mass loss trends that covered a range of signal magnitude from ∼15 to 75 Gt yr −1 . These fluxes were subsequently incorporated into the global hydrology model (see Sect. 3). In the rest of this Section we focus on how the ice sheet fluxes were prescribed.
Two processes that have different spatial and temporal behaviours control the ice sheet fluxes: surface mass balance (SMB) and ice dynamics (ID). The former is determined by the balance between solid precipitation and runoff and responds almost instantaneously to changes in external forcing. Inter-annual variability in SMB is large compared to, for example, secular trends in mass loss and it is important, therefore, that this is adequately simulated in addition to any underlying trends.
Ice dynamics is controlled by the thermal and basal regime of the ice mass and changes in boundary conditions such as at the grounding line of marine terminating glaciers. An important and interesting question is whether it will be possible to separate SMB from ID signals with sufficiently high  resolution (both in time and space) gravity data and our aim was, therefore, to ensure that both these signals were realistically simulated to explore this issue.

Short description of the model and adaptation to the study
For the ice sheets, a time series of surface mass balance fluxes were calculated for the period 1995-2005 using ECMWF reanalysis data (ERA-40 up to 2001) and the operational analysis beyond. The data were produced on a 5 km polar stereographic grid with 6-hourly time step for Green-land and a daily time step for Antarctica, as the diurnal cycle is only important for determining melt in Greenland. The SMB, which comprises accumulation-ablation, was calculated using a regional climate model for Antarctica (Van de Berg et al., 2006) and a simpler downscaling of the ECMWF data combined with an snow diagenesis model for Greenland (Bougamont et al., 2005). Superimposed on the SMB were monotonic, secular trends in ice dynamics (i.e. solid ice fluxes into the ocean), reflecting changes that are known to be taking place from a range of satellite observations (Rignot and Kanagaratnam, 2006;Rignot et al, 2008). The magnitude of these changes are representative of the real world signals but are also designed to capture a range of signal amplitudes and spatial scales. An example of the deviation of the year 2000 SMB from the eleven-year mean is shown in Fig. 3, which also indicates the pattern of ice dynamic losses imposed for both ice sheets. The areas shaded black indicate regions where a secular trend in ice dynamics was imposed based on the 50 m a −1 velocity contour for the outlet glacier. The spatial pattern and magnitude of the trends are comparable with recent observations of regional mass loss from, for example, SAR interferometry. In addition to the spatially distributed mass changes over the ice sheets we also calculated the flux of ice at the margins entering the ocean grid cells. This was done with the aid of the OMCT ocean/land mask. For Antarctica the fluxes were calculated annually as there is no seasonal or daily variation so that the daily value was 1/365 of the annual value. For Greenland the margin fluxes were calculated on a daily time step due the seasonality of surface runoff (cf. Fig. 5). Both Greenland and Antarctica were excluded from the hydrology model as explained in Sect. 3.

Description of resulting Greenland and Antarctic ice mass fields
The daily mass trends integrated over the whole of the Greenland Ice Sheet are shown in Fig. 5 starting on 1 January 1995.
Positive values indicate snowfall and an increase in SMB that is evident in winter and spring. During summer, surface ablation dominates and results in a negative SMB. The mean is less than zero even when there is no surface ablation due to continuous solid ice discharge, which was assumed to have no significant seasonality. The cumulative mass trends are shown in Fig. 6 and indicate the amplitude of the seasonal cycle in SMB which is around 300 Gt and is in good agreement with results from   . Starting at about day 1000, we imposed an increase in ice discharge and ablation that results in a roughly constant loss of mass to the oceans of 110 Gt a −1 . This rate is comparable with mass balance estimates of the ice sheet between the mid 1990s and mid 2000s ( . The spatial pattern of loss is distributed between increased runoff and discharge as shown in Fig. 3 left.

GRACE (Van den
The spatial and temporal pattern of mass loss imposed for Antarctica is markedly different to Greenland. The variability in SMB is determined by the regional climate model output and is shown in Fig. 7 at annual resolution. There is no secular trend in this component over the time period considered. Ice discharge was increased linearly for the Amundsen Sea Sector of West Antarctica (Pine Island and Thwaites glaciers in Fig. 3 right) from 1995-2000 and then kept constant thereafter. The total increase was 100 Gt over the five year period (Fig. 7). This results in both negative and positive balance years depending on the snowfall in a particular year. The mass loss for the whole ice sheet is somewhat smaller than recent observations suggest (Velicogna, 2009) but the spatial distribution and magnitude of dynamic losses in the Amundsen Sea sector are quite close to observations (Rignot et al., 2008).

Ocean
The global Ocean Model for Circulation and Tides (OMCT) (Thomas, 2002) has been developed by adjusting the original climatological Hamburg Ocean Primitive Equation model (HOPE; Wolff et al., 1996;Drijfhout et al., 1996) to the weather time scale and coupling with an ephemeral tidal model. OMCT has been used to investigate the impact of the general ocean circulation on the Earth's rotation and the gravity field Wünsch et al., 2001;Dobslaw andThomas, 2005, 2007b) and is currently applied to de-alias short-term ocean mass anomalies within the GRACE gravity field processing (Flechtner, 2007).
The model is based on non-linear balance equations for momentum, the continuity equation, and conservation equations for heat and salt. The hydrostatic and the Boussinesq approximations are applied. Water elevations, threedimensional horizontal velocities, potential temperature as well as salinity are calculated prognostically, vertical velocities are determined diagnostically from the incompressibility condition. Implemented is a prognostic thermodynamic seaice model (Hibler III, 1979) that predicts ice-thickness, compactness and drift. Effects of loading and self-attraction of the water column are parameterized by means of a secondary potential proportional to the mass of the local water column . The model uses a time-step of 30 min, a constant horizontal resolution of 1.875 • in longitude and latitude, and 13 layers exist in the vertical.
A quasi steady-state circulation has been obtained from an initial model spun up for 265 yr using climatological wind stresses (Hellerman and Rosenstein, 1983) and mean sea surface temperatures and salinities according to Levitus (1982). Subsequently, OMCT has been forced by 6hourly wind stresses, atmospheric surface pressure, 2 m temperatures and freshwater fluxes due to precipitation, evaporation and runoff obtained from the latest reanalysis project ERA-40 of the European Centre for Medium Range Weather Forecasts (ECMWF), and the Hydrological Discharge Model HDM (Walter, 2007) for the period 1958-1989. The following 5 yr have been simulated using ERA-40 atmospheric data accompanied with climatologically averaged ice-sheet discharges (Sect. 4) and daily continental runoff (Sect. 3) in order to prevent numerical shocks due to sudden changes in forcing conditions. The period relevant for the data-set described here, i.e. 1995 until 2006 has been simulated by applying forcing based on ERA-40 (1995ERA-40 ( -1999 and operational ECMWF atmospheric data (2000)(2001)(2002)(2003)(2004)(2005)(2006) together with daily discharges from continental hydrology (see Sect. 3), daily discharges from the Greenland ice-sheet, and annual ice discharges from Antarctica (Sect. 4).
Total ocean mass variations due to freshwater fluxes are considered up to seasonal time-scales. However, as suggested by Greatbatch (1994), artificial mass variations caused by the applied Boussinesq approximation have been corrected for (cf. Dobslaw and Thomas, 2007a). Oceanic mass fields have been stored every 6-h.
Ocean mass anomalies have been obtained by reducing the mean field averaged over 2000-2005. Due to the shift in atmospheric forcing conditions from ERA-40 towards operational ECMWF analyses in 2000, mean atmosphere and ocean mass fields change at this time substantially, as can be deduced from the differences of mean atmospheric pressure as well as ocean mass for the periods 1995 until 1999 and 2000 until 2005 (Fig. 8). Differences over continental regions are mainly due to the improved spatial resolution in the operational model, allowing a more detailed representation of the orography. The largest deviations over the oceans occur mainly in coastal regions, which are both related to artifacts caused by the required horizontal interpolation towards the OMCT grid as well as to the different spatial resolution. However, large-scale differences with small but still relevant amplitudes exist in the open ocean, e.g. in the southern Pacific, which might lead to systematic biases. In order to remove this bias as much as possible, a mean field for the ERA-40 period has been calculated by minimising the differences between mass anomalies from simulations forced by ERA-40 and the operational ECMWF data in the overlapping years 2000 and 2001, separately for atmospheric and oceanic mass anomalies. The obtained mean field has been removed from all ERA-40 forced data covering 1995-1999. Variability of regional ocean mass anomalies is essentially comparable in the run analysed in Dobslaw and Thomas (2007a). The highest variability occurs in the Southern Oceans as well as in the North Pacific, and, with slightly smaller amplitudes, in the Arctic basin. Small scale signals with significantly higher amplitudes additionally occur in various shelf regions but are generally less reliable due to the 1.875 • horizontal resolution of the OMCT. While variability is fairly comparable for both simulated periods, apparent bottom pressure trends alter significantly when changing from ERA-40 to operational ECMWF forcing fields (Fig. 9). At this point, causes for these apparent changes in the trends are not clear and remain the subject for further studies.
Atmospheric and continental freshwater fluxes affect the total mass contained in the global oceans (Fig. 10). Since atmospheric freshwater fluxes as derived from numerical weather prediction data contain large uncertainties, estimations of total ocean mass by means of numerical models are particularly challenging. ECMWF freshwater fluxes have   been intensively analysed (e.g. Andersson et al., 2005) revealing significant overestimations of precipitation in tropical ocean regions. This leads to a net-influx of water into the ocean, which is obviously unrealistic (see, e.g. Dobslaw, 2007, p. 40ff). Therefore, an additional correction has been applied by requiring that net-fluxes into the ocean must be zero within each year. This condition has been achieved by adding a time-invariant and globally homogeneous layer of fresh water. Thus, variations in total ocean-mass are included on seasonal and sub-seasonal time-scales only.

Solid Earth
Time variable gravity field changes driven by solid-earth processes are primarily caused by glacial isostatic adjustment (GIA) due to continental ice mass changes and concomitant sea-level variations, co-and post-seismic deformation, mantle convection and plate tectonics, and core motions and seismic modes (e.g. Vermeersen, 2005). Specifically GIA and Th. Gruber et al.: Simulation of the time-variable gravity field deformations caused by very large earthquakes (above 8 on the Richter scale) lead to gravity field changes that are observable by space-borne gravimetric missions (e.g. Han et al., 2006;Einarsson et al., 2010), which have a typical foreseen lifespan of a few years to a decade.

Global Isostatic Adjustment
Glacial Isostatic Adjustment (GIA) is the response of the solid Earth to the waxing and waning of Late-Pleistocene Ice. The melting of the ice has left large and persistent geoid and gravity anomalies over regions like Canada and Fennoscandia. Unfortunately, separating the GIA-induced contributions to the geoid from those induced by plate tectonics and mantle dynamics is not always straightforward. For example, it is now widely acknowledged that the deep geoid low above Canada is partly due to non-GIA induced lithosphere and mantle heterogeneities and is only partly attributable to GIA (e.g. Simons and Hager, 1997;Tamisiea et al., 2007). Whereas the geoid above Canada is related to two geodynamical processes, it is thought that secular geoid and gravity anomaly variations are only triggered by post-glacial rebound (e.g. Wahr and Davis, 2002). These secular changes have been included in the temporal gravity field model.
Spherical harmonic coefficients were derived complete to degree order 180 based on a selected Earth and ice model. These coefficients represent the secular change due to GIA. Each spherical harmonic or Stokes coefficient C is represented by: where C 0 is the start value of the Stokes coefficients (from the other contributions to the gravity field model), C GIA is the gravity field linear trend value (change per year), t is time (yr), and t 0 (yr) is the reference time for the static model. The standard elastic Earth model used is PREM (Dziewonski and Anderson, 1981). The thickness of the lithosphere is 120 km, the radial mantle viscosity profile VM-2 (Peltier, 2004) and the standard ice model ICE-5G (Peltier, 2004). Sea-level variations induced by GIA are self-consistently determined by means of solving the sea-level equation (Sabadini and Vermeersen, 2004). Rotational feedback is incorporated in the way described in Chapter 4 of Sabadini and Vermeersen (2004). Details on how lithosphere thickness and mantle viscosity affect geoid heights in terms of spherical harmonics and their detectability by GOCE and GRACE can be found in Schotman and Vermeersen (2005) and Vermeersen and Schotman (2008). The resulting spherical harmonic model is displayed in Fig. 11 in terms of geoid change. It can be observed that the dominant changes indeed occur in Canada and Fennoscandia, but also over Antarctica, with a maximum amplitude of a few mm yr −1 . However, as has already been dealt with in Sect. 4 on ice, for Antarctica also contemporary ice mass variations induce non-negligible geoid signals, which are hard to separate from GIA (e.g. Riva et al., 2009).

Earthquake
Large earthquakes induce local, regional and global gravity field variations, both during and in the days, months, years and even decades after the faulting event. During a faulting event there is an immediate, non-recoverable redistribution of the Earth's mass. This is called co-seismic deformation. Due to the existence of shallow low-viscosity intra-crustal and asthenospheric layers, the redistribution of stress and strain due to the faulting will relax in the days, months, years and decades after the earthquake. This relaxation does not necessarily diminish the co-seismic mass redistribution; the post-seismic deformation can, in fact, enhance the co-seismic component. Co-and post-seismic deformation, and thereby gravity variations, due to large earthquakes depend on the parameters of the earthquake source, such as seismic moment (i.e. the product of the solid-Earth rigidity at the fault, the fault length and relative fault displacement), the type of earthquake (e.g. normal fault, strike-slip fault), the geometry of the faulting event and depth of the earthquake (e.g. Sabadini and Vermeersen, 1997).
One very large earthquake (above 9 on the Richter scale) was the Sumatra-Andaman earthquake, which hit on 26 December 2004. This earthquake has been selected for inclusion in the temporal gravity field model. As with the other mass change fields, a spherical harmonic coefficient model complete to degree and order 180 was derived. The standard elastic Earth model is again PREM (Dziewonski and Anderson, 1981). Sea-level variations induced by co-and post-seismic deformation can be determined selfconsistently. Spherical harmonic expansions were derived for both co-seismic (episodic) and pos-seismic (trend) gravity field changes Vermeersen, 1997 and2004). The co-seismic contribution to the gravity field model is considered to be a step function. The gravity field model changes abruptly after 26 December 2004 according to the co-seismic related 180 × 180 spherical harmonic expansion. The postseismic gravity field change is modeled by a 180 × 180 spherical harmonic expansion as well, where the associated Stokes coefficients represent a linear trend for the first year. After this one-year period, post-seismic relaxation is considered to be negligible.
Thus, Stokes coefficient changes are modeled as: Before 26 where t 0 is the time of the earthquake (yr), C co-Sum represents the co-seismic gravity change, and C post-Sum is the gravity field linear trend value due to post-seismic relaxation (change per year). The associated co-seismic geoid changes are strongly localized close to the earthquake epicenter and are between −30 and +10 mm (Fig. 12, left). The post-seismic relaxation is about a factor of 10 lower for a one-year period (Fig. 12, right). If the spherical harmonic expansion is truncated at degree and order 50, the maximum geoid change is about −10 mm, which is comparable to the values derived from monthly GRACE solutions, which typically provide statistically signals to this truncation degree (Einarsson et al., 2010).

Combined mass fields
As described in the previous sections, estimates of mass variations are provided via the geophysical models for each component of the Earth mass system and in different representations. The different representations include (but are not limited to): pressure, equivalent water height, and gravitational attraction. In addition, sometimes these fields are referenced to specific mean fields. These values are provided on regular or irregular grids or are provided in terms of global spherical harmonic series. In order to compute our "real" time variable Earth gravity field, all mass input fields have to be harmonized. This includes harmonisation of the reference fields, transformation to identical quantities, interpolation to unique regular grids and in some cases interpolation to the 6-hourly time intervals. After this stage, the combined mass fields are converted to 6-hourly gravity potential spherical harmonic series, which are the usual representation for the global gravity field (Heiskanen and Moritz, 1967). In the following sub-chapters we provide a detailed description for the pre-processing of the geophysical mass fields and the conversion to spherical harmonics.

Geophysical model mass fields
The starting point for pre-processing the geophysical model mass fields were the atmospheric pressure fields. As they are globally available on a 1 × 1 • equi-angular grid for every 6 h time interval it was decided to use this as a reference and to convert all other mass fields into this representation. This implies that for the atmospheric mass fields no further processing is required. The global hydrological model (PCR-GLOBWB) provides daily estimates of the total continental water storage on a half degree equiangular grid in terms of volume of water (see Sect. 3). In order to convert this into 6-hourly 1 degree equiangular grids of pressure a number of processing steps had to be performed. First, volume is converted into pressure by the well-known relation: The area size for each cell is provided together with the data set, while density and gravity are assumed to be constant: density = 1000 kg m −3 , gravity = 9.81 m s −2 . In order to determine 6-hourly continental water mass estimates a linear interpolation in time for the daily half degree grids is performed. Then 1-degree block mean grids are computed from the mean of all defined 30 block mean values that are located inside the 1-degree block. Finally, all files are converted to the internally used grid format. After this step, a preliminary check of the grid time series in terms of a trend analysis was performed. At this point, some unrealistic outliers were identified and eliminated. These outliers are located in inland water bodies (lake Titicaca and Caspian sea) and in Greenland. It is well known that the PCR-GLOBWB model does not perform very well in these areas and therefore an elimination of these grid points is recommended. The edited time series then was used for combination with other mass fields.
The continental ice mass model provides ice mass variation estimates for Antarctica and Greenland separately in a polar-stereographic representation on a daily basis (see Sect. 4). Mass variation changes are expressed in terms of equivalent water height change with respect to the previous day. As a first step, the polar stereographic coordinates are converted to geographic latitude and longitude values. Then 1-degree block mean values are computed for each day by a simple mean value operator applied to all data points located inside a 1-degree block. Due to the projection, for a few 1degree grid points in Antarctica there are no data points available. For those, a linear interpolation in longitude was performed. These points are all located close to the pole where very small mass changes occur. Equivalent water height is converted to pressure by applying the following relation: with: EWH Equivalent Water Height in [m]. Like for continental hydrology, density and gravity are assumed to be constant with the values defined above. Daily 1degree grids are then linearly interpolated in time in order to determine 6-hourly ice mass fields. Finally, ice mass changes in terms of pressure are accumulated in order to compute the total variation with respect to the first time step. This means that all ice-mass variation files finally are referred to the 1 January 1995 00:00 UTC. For this reference time, the total ice mass is assumed to be 0 and only variations are available.
Oceanic bottom pressure fields resulting from the OMCT model (see Sect. 5) are available for the same 6-hourly time steps as the atmospheric pressure. They are scaled with a constant factor of 9.7 × 10 −4 , which has to be reversed in order to determine real ocean bottom pressure in [Pa]. As the horizontal resolution of the OMCT model is limited to 1.875 degree a further interpolation to a 1 × 1 • grid was required. For this, a standard procedure based on 2-D linear interpolation was used. As with all other data sets, all ocean bottom pressure files are converted to the internal grid format.
For modelling solid-Earth mass variations, gravity field variations due to post-glacial rebound and due to the Sumatra earthquake were taken into account (see Sect. 6). All together these sources represent three spherical harmonic series that describe the changes of gravity field variations. 6hourly gravity field spherical harmonic series for the solid Earth mass variations are computed with the formulas provided in Sect. 6. As the temporal reference, 1 January 1995 00:00 UTC is again used. Subsequently, these spherical harmonic series coefficients are converted to 6-hourly pressure variations on a 1-degree equiangular grid by spherical harmonic synthesis. During the synthesis, the conversion to pressure is done by applying the necessary factors (see formula in Sect. 7.2). In addition the coefficients are divided by the loading Love numbers, because they have been included in the original series and they will be applied again in the final spherical harmonic analysis. In this way, it is ensured that no double effect of loading is applied. Using this procedure, a complete time series for the solid Earth mass variations in terms of pressure on a 1-degree equi-angular grid is obtained.

Gravity potential spherical harmonic series
Starting from the pre-processed mass variation fields, as explained in the previous section, gravity field variations for the complete time series are computed. All source fields are available as 1-degree equiangular grids in terms of pressure with 6-hourly temporal resolution for the period from 1 January 1995 00:00 UTC to 31 December 2006 18:00 UTC. Gravity potential spherical harmonics are determined by harmonic analysis of the pressure fields applying the following formula. Using this integration formula, harmonic coefficients up to degree and order 180 (corresponding to the 1-degree spatial resolution) are computed for various combinations of geophysical mass fields with respect to a mean field. The mean field for each geophysical model data set is defined as the state on 1 January 1995 00:00 UTC. This implies that starting from this date variations are computed with respect to the mass field at this date. The size of each surface element Figure 13. Mean of monthly signal in terms of spherical harmonic degree variances for different mass variation fields compared to mean monthly error from GRACE time variable gravity field estimates. is taken into account by dS in the equation above. Loading Love numbers are applied in order to take into account the loading of the mass on the Earth's crust and the integration is performed by summation over all surface elements. In order to avoid leakage effects, it is important to apply global data sets. For this reason, the atmospheric mass field always represents the basic quantity, because it is the only global field available. Then step-by-step additional sources are added and the final time series represents the simulated global time variable gravity field. The following data combinations were performed:

Results
In this section, a few examples resulting from this study are shown. Here we only use the resulting spherical harmonic series of the time variable gravity potential in order to characterize them. One way to identify the signal strength of the simulated signal is to compare it with the results for the time variable gravity field signal as determined from the GRACE mission. Sensitivity on a global scale can be analyzed by signal and error degree variances (Heiskanen and Moritz, 1967). They show the mean signal or error per degree of a gravity field spherical harmonic series. When comparing the signal degree variances of the simulated fields with the error degree variances of the gravity field time series as observed by GRACE (in this case the ITG-GRACE2010S time series is used, see Mayer-Gürr et al., 2010) one can identify, if a mission is sensitive to the specific time variable gravity signal or not (see Fig. 13). Figure 13 tells us up to what resolution (in terms of spherical harmonic degree) the time variable signal from individual sources is sensitive for a mission like GRACE. As this representation is a global average one should regard it as a pessimistic estimate, which means that for regions with strong mass variability sensitivity could be even higher. We can identify that hydrology exhibits the strongest signal, while for all other sources the signal has a significantly reduced power. With GRACE, one can estimate monthly changes in the continental water at least up to degree 40 (corresponding to 500 km spatial resolution), while for other sources the error line crosses the time variable signal at about degree 20 (corresponding to 1000 km spatial resolution). When interpreting these results one should not forget that there is still some room to improve the GRACE error curve in future, which would mean higher sensitivity with respect to the individual signals. As this study was designed in preparation for defining future gravity field missions, it is expected that with new gravity field sensors in orbit (laser interferometer) one could reduce the error curve of GRACE by a factor between 20 and 50, which will strongly improve the sensitivity of future missions in terms of spatial resolution. In order to identify the signal versus error for the individual time variable gravity field sources for other than monthly periods, the mean variability was computed for 6 hourly, daily, weekly and semi-annual time intervals (see Fig. 14).
From Fig. 14, one can clearly identify which of the mass fields has higher variability for shorter time intervals and vice versa. Obviously, for the atmosphere-ocean system daily variability is nearly as strong as monthly and seasonal signals, which implies, that for the analysis of gravity field observations from GRACE or any future mission one has to use geophysical models for removing the high-frequency variability during processing (de-aliasing) or one has to observe them directly. Hydrology and ice fields exhibit much slower variability. To estimate hydrological or ice mass changes from space, one needs monthly, or for specific areas with strong signals, weekly observations. Solid Earth processes play a slightly different role, because except for co-seismic events these are very slow processes without seasonal variability. Sudden co-seismic events like earthquakes could have a significant local mass change signal, which is, because of the global mean operator applied for the computation of degree variances, not clearly observable in such representations. In order to have a closer look to the spatial distribution of mass changes, a few monthly maps were computed from the time variable gravity field spherical harmonics (see Fig. 15). Figure 15 identifies, for two representative months, areas where we have monthly mass change with respect to the yearly mean. From the combined mass fields (top row), one clearly observes seasonal hydrological variations in large river basins like the Amazon, Ganges and Zambezi as well as some variability for the large continental ice sheets in Greenland and Antarctica. Due to the high frequency nature of atmosphere and ocean variability, this mostly is cancelled out when computing the monthly mean. Long-term variability from solid Earth processes hardly is visible in this representation, but would become evident when comparing monthly means with a decadal or longer reference field. Comparing the results from the simulated gravity field variations with those observed from GRACE (bottom row), for areas with large signal amplitudes good correlations can be identified (global correlation up to 50 %). Due to the design of the GRACE observing system (polar orbit with one dimensional observations mostly in North-South direction) we have a non-isotropic error behaviour, which becomes visible by the North-South stripes in the GRACE solutions. Applying specific a-posteriori filter techniques (cf. Kusche, 2007) one could reduce these artefacts significantly with the drawback that some signal is reduced as well. As there is a wide range of filters available, for the correlation analysis shown here the originally observed monthly fields were used. This ensures that the signals' amplitudes from the geophysical models are comparable to the results from the GRACE data analysis.
From the kind of comparisons shown in the figures above, one can observe that the simulated time variable gravity field as computed from the coupled geophysical models, to a large extent, contains realistic signal amplitudes and frequencies making it is well suited for use in simulating future gravity field missions.

Conclusions
In this paper, we have detailed the development of a selfconsistent model that realistically describes the motion of mass between the fluid reservoirs of the Earth system. Our model captures mass redistributions between the atmosphere, hydrosphere, oceans, cryosphere and solid Earth and combines them into 6-hourly Stokes Coefficients up to degree and order 180 for the period between 1995-2006. The data have been used to test various mission scenarios, e.g. satellite orbit parameters, allowing us to choose scenarios that optimize the observability and separability of specific mass transport features. In addition, using a known data set such as this as input, has allowed us to determine the effects of different types of satellite error contributions on the gravity solutions.