Sensitivity of Goce Gradients on Greenland Mass Variation And Changes in Ice Topography

Abstract The Gravity field and steady state Ocean Circulation Explorer (GOCE) maps variations in the gravity field by observing second order derivatives (gradients) of the Earth gravitational potential. Flying in the low altitude of 255 km and having a spatially dense data distribution of short wavelengths of the gravity field, GOCE may be used to enhance the time varying gravity signal coming fromthe GRACE satellites. The GOCE gradients may potentially be used for the determination of residual masses in local regions. This can be done using Least-Squares Collocation (LSC) or the Reduced Point Mass (RPM) method. In this study, different gravity field solutions are calculated by the use of RPM, LSC and GOCE gradients, respectively. Gravity field time series are created and presented for the six consecutive months of GOCE gradient observations, data being acquired between November 2009 and June 2010. Corresponding gravity anomaly results are used for the calculation of ice mass changes by the use of theRPMmethod. The results are then compared with the computed topographic effect of the ice by the use of a modified topographic correction and the Gravsoft TC program. The maximal gravity changes at the ground predicted from GOCE gradients are between 2 and 4 mGal for the period considered. The gravity anomaly estimation error arising from the GOCE gradient data using only Tzz with an associated error of 20 mE is 11 mGal. This analysis shows the potential of using GOCE data for observations of ice mass changes although the GOCE dataset is limited to only six months. We expect four years of GOCE gradient observations to be available by mid-2014. This will increase the accuracy and spatial resolution of the GOCE measurements, which may lead to an accuracy necessary for observing ice mass changes.


Introduction
The Gravity eld and steady state Ocean Circulation Explorer (GOCE) is the most advanced gravity space mission to date, and since 2009 it has mapped global variations in the gravity eld with remarkable detail and accuracy.The success of the GOCE mission depends on adequate methodologies being developed for extracting the gravity eld from its observations and for combining the gravity eld with information from other sources.In such a way, GOCE derived data has a possibility to add the insight needed in order to further our understanding of the physical processes occurring inside the Earth and on its surface.This can be done through a broad range of application, e.g. in the elds of geodesy, solid Earth physics, oceanography and cryospheric research.
In order to get a better picture of the e ect of the present-day climate changes in polar regions, research has been carried out focusing on the regions where changes in the Greenland ice sheet have been observed, e.g.Krabill et al. (2004), Velicogna and Wahr (2005), Joughin et al. (2010).Most of the largest outlet glaciers have experienced accelerated retreats, causing an increased ice discharge (Sørensen, 2010).The largest accelerations have been observed by the Helheim glacier, Kangerdlugssuaq and Jakobshavn Isbrae (JI) (Howat et al., 2011).Lowering at rates of 30-35 m/yr (Levinsen et al., 2013), JI has the potential to in uence sea level rise more than any other single feature in the Northern Hemisphere.
A recent GRACE study by Velicogna and Wahr (2013) shows that using Release 5.0 (RL05) GRACE elds for January 2003 through November 2012, a mass change of − ± Gt/y for Greenland is found, and a loss that migrated clockwise around the ice sheet margin to progressively a ect the entire periphery.Comparison of di erent studies for recent mass loss in Greenland (Svendsen et al., 2013) shows local disagreement between the data products.However, all models agree that the acceleration in mass loss is largely con ned to the west-northwestern part of Greenland.
Establishing the presence of an acceleration on the order of magnitude found in the Greenland Ice Sheet re-quires more than 5 years of data, and we nd that the GRACE time series available are now long enough to establish the presence of such an acceleration.Maximally four years of reprocessed GOCE gradient observations will to be available by mid-2014, which may add the supplement GRACE derived solution for mass changes.
The GRACE satellites are 450 km above the Earth's surface, which makes them relatively insensitive to shortscale terms in the gravity eld.This means that the processing centers remove short-scale terms by truncating their solutions to a nite set of low-degree harmonics 60 for Release 5.0 (RL05) elds calculated by CSR (the Center for Space Research at the University of Texas) and 90 for solution calculated by GFZ (GeoForshungsZentrum in Potsdam), see Velicogna and Wahr (2013).This results in spatial resolution of approximately 330 km.Gradient data from the GOCE satellite mission, having the spatial resolution of approximately 100km, may improve the resolution of observed ice mass changes.
The full range of possibilities for the use of GOCE gradients are still not fully recognized.Even though ocean studies already bene t from the GOCE gradients (Herceg, 2012), their contribution to the solid Earth applications is still to follow.
The aim of this study is to present the methodology for using GOCE gradients to possibly extract the signal that corresponds to the signal of change in JI ice mass.In order to present GOCE capabilities of detecting ice mass change by using second order derivatives of gravitational potential observed by GOCE, the used time interval spans from November 2009 to June 2010.

GOCE T zz gravity gradients
As the result of the processing of satellite data, the GOCE mission provides many products.One of these products are gravity gradients given in the Terrestrial Reference Frame (TRF).The available GOCE level 2 TRF data products, which are o ered to the GOCE user community, are packed in monthly and sub-monthly packages (Table 1) in the Earth Explorers File Format Standards (EEF) format (Floberghagen et al., 2010).
The GOCE gradient dataset contains 6 elements of full gravity gradient matrix with associated errors.In this study only Tzz GOCE gradients were used, since adding the Tyy GOCE gradients component when using the Least-Squares Collocation (LSC) method only provides marginal improvements to the results obtained when solely using Tzz (Tscherning and Arabelos, 2011).
In order to produce the gravity anomaly residuals, the long wavelength part of the gravity eld must be therefore subtracted.In this study, the gravity eld contribution up to harmonic degree and order 36 is subtracted from the GOCE gradients (Fig. 1).GOCE Tzz gravity gradients in the Greenland region range between − .and + .Eötvös.Strong gravity gradients (around + .Eötvös) are present in the southern part of Greenland and along the central-eastern parts of the ice sheet (Fig. 1a).This indicates excess of masses compared to the reference model (gradients computed from normal the potential of the reference ellipsoid), an excess, whose source may come from the Earth's interior or from the Greenland Ice Sheet.Near-zero and weak positive Tzz gravity gradients (up to 0.2 Eötvös) are present in most parts of Greenland, while weak negative gravity gradients are only present in the JI region.
Removing the contribution up to spherical harmonic degree and order 36 subtracts the contribution of the gravity eld signal with wavelengths longer than 500km.This allows calculations in local regions, since the results are not contaminated by the gravity eld outside of the study area.For the truncation of the GOCE gravity gradients, the EGM2008 global geopotential model is used (Fig. 1b and Fig. 2).The EGM2008 up to degree 2190 has been developed by Pavlis et al. (2012), and it incorporates shipborne, airborne, and satellite altimetry derived gravity anomalies.It also has bene ted from the Gravity Recovery and Climate Experiment (GRACE) based satellite solutions.
Truncation of the GOCE Tzz gradients, yields the corresponding gradient anomalies (Fig. 1b for combined months and Fig. 2 for monthly solutions).Calculated gradient anomalies for the Greenland region range between − .and + .Eötvös, which indicates a gradient signal of around + .Eötvös, which was removed in truncation.Comparing the full-spectrum GOCE Tzz gradients and the truncated result, it is clear that truncation has the strongest e ect in the southern part of Greenland as well as along the central-eastern coast, i.e. places where the full spectrum gradient signal was strongest.
Monthly solutions for the GOCE Tzz gradient anomaly show a very similar gravity signal (Fig. 2).However, some months have data gaps that leaves the part of Greenland without data coverage (Fig. 2d for February and Fig. 2e for March).
The spatial resolution of GOCE data is around 100 km (this can be observed on Fig. 2), while that of GRACE barely exceeds 330 km.The high spatial resolution of GOCE observations may lead to more detailed mapping of ice mass changes.
High Tzz gradient anomalies (approximately + .Eötvös) are still present in the southern part of Greenland (Fig. 2), while the area along the central-eastern coast shows near-zero and weak negative gradient anomalies.Very signi cant is that strong negative gradient anomalies are present in the area near two of Greenland's largest out-let glaciers, i.e.JI and the Helheim glaciers.This may reect low-density masses in these regions or rather a lack of masses compared to the reference model and the removed contribution up to spherical harmonic degree and order 36.

Prediction of the gravity anomaly using GOCE T zz gradients
For regional gravity anomaly prediction in this study, two methods are used, one of them being LSC (developed by C. C. Tscherning, 1974) and the other the Reduced Point Mass (RPM) method (developed by Herceg, 2012).Both methods show good agreement in prediction of the gravity anomaly residuals, which will be demonstrated in the comparison section.Compared to the LSC method, the RPM method needs less time for calculation since the number of equations needed to be solved in LSC depends on the number of observations.With RPM, the number of equations is dependent on the number of gridded point masses.

. The Least-Squares Collocation method
LSC enables the use of many types of observables (Krarup, 1969) for the estimation of gravity eld quantities and their errors.
The basic observation equation for LSC is where e i is the error contribution.The estimate of T LSC is obtained by   where TLSC (P) is the estimated quantity, C Pi is a vector which contains the covariance values between the estimated quantity and observations, and C = C ij + σ ij , where C ij and σ ij are matrices containing the signal covariance and the error covariance values, respectively, between the observations y (Tscherning, 1976).
Analytic covariance models are used in LSC for regional applications.The optimal covariance function for a region is selected by tting an analytic model to empirically determined values.Here, the empirically determined values are obtained by using the program EMPCOV (Forsberg and Tscherning, 2008) on data in a certain region.
To t empirical covariance functions to isotropic analytic models, the Fortran program cov t.f from GRAVSOFT package was used (Knudsen, 1987).The parameters obtained by tting the covariance function were used in the collocation prediction method as well as the reduced point mass prediction.
Point-mass functions are harmonic functions, which may be used to represent T either globally or locally, and they can be expressed by either closed expressions or as sums of Legendre series.In both cases at least the two rst terms must be removed since they are not present in T, i.e. zero-degree (l=0) harmonic terms cancel each other out and the rst degree harmonics (l=1) are equal to zero since the origin of the coordinate system is chosen to coincide with the geocentre.For local applications the e ect of a global gravity model is generally removed and later on restored.
The anomalous gravity eld, T, at computational point Q is modeled by a set of base functions, each obtained as the anomalous gravity potential from each point mass m i located at the position P i on the Bjerhammer sphere with radius R M estimated with using LSC method.This radius of point masses, R M , is smaller than that of the Earth, R E .
where V and U are the Earth's gravity and normal potential, respectively.
Normally, a gravitational potential function is related to the position in terms of spherical coordinates (ϕ, λ, r).
If the potential is considered, V = V(λ, φ, r), the gradient of V is de ned in the geocentric (x, y, z) or local rectangular (η, ξ , ζ ) Cartesian coordinate system in the normalized basis as (Herceg, 2012): Here the local coordinate system is described by the basis ep with vectors in triad (e , e , e ), which are not unit vectors.
When using point masses, the closed expression for the vertical second order derivatives of the potential in the (η, ξ , ζ ) system are given as (Herceg, 2012): where t is the cosine of the spherical distance ψ between two points (P and Q), and the partial di erential equation of second order satis es the Laplacian di erential equation, ∆V = .
Reduced gravity potential, V R is represented through reduced point masses, which are point mass functions where the rst n harmonics have been removed.For the expressions of gravity gradients, when using reduced point masses, the derivative of the sum of a nite Legendre series is used: Here, m is degree of the Legendre polynomial P, s contains the parameters or the anomaly degree variance model, t is the cosine of the spherical distance ψ between two points, and a l represents the multiplication factor, which depends on the derivative with respect to r Q .For these, we have series expansions similar to Eq. ( 6) where the terms for the rst derivative are multiplied with a = −(i + )/r Q and for the second derivative with a = (i + ) (i + )/r Q .The derivative of the sum of a nite Legendre series can be computed easily using a recursion algorithm (Tscherning and Rapp, 1974).For  Thus, the expression for the vertical second order derivative of the gravitational potential in the ζ direction is (Herceg, 2012): Here, the partial derivatives of the spherical distance ψ with respect to the latitude and longitude are expressed through t, as in Herceg (2012).
This approach is able to tailor the algorithm for point mass depth and grid spacing relations.The method provides the calculation of both full and reduced gravity eld quantities using either full or reduced point masses, respectively.Gravity anomaly determination by means of the reduced point masses can be used as an alternative method to the conventional ones for local gravity eld determination.
The prediction of reduced point masses be used for the estimation of absolute masses in local regions, but may supplement the determination of relative mass changes.The gravity gradient anomalies can then be used in the RPM method for the calculation of ice mass change in a given time interval.

. Comparison of the LSC and RPM results
The predicted gravity anomaly residuals found using collocation have a standard deviation (std) of .mGal at the surface, while the std of the RMP results is 35.32 mGal.Both methods presented here show good agreement in the prediction of gravity anomaly residuals (Fig. 3).Strong gravity anomaly signals are present mostly in the southern and north-eastern parts, where it exceeds 100 mGal.Negative gravity anomaly residuals are present mostly near In the following section, the RPM method will be used for determination of the gravity anomaly change occurring from November 2009 -June 2010 and corresponding relative mass changes in Greenland.
The results of the gravity anomaly changes will then be compared to the gravity anomaly changes resulting from changes in ice mass topography calculated from a Digital Elevation Model (DEM) for the JI drainage basin.
The method and data used for the prediction of gravity anomaly changes by in ice mass topography is described in the following section.

Prediction of mass changes using GOCE T zz gradients and the RPM method -Jakobshavn Isbrae drainage basin
Even though an observed mass loss in Greenland has been substantial over the last decade (Fig. 4), the largest lowering of the ice surface is observed over the largest outlet glacier on Greenland's west coast, namely Jakobshavn Isbrae (JI) (Fig. 5).This glacier loses the most mass and has doubled its contribution to global sea level rise over the last decade (Howat et al. (2011), Khan et al. (2010), Joughin et al. (2014)).Greenland mass loss observed using gravity gradients may supplement the methods for observing surface elevation changes in determining whether the ice surface is thinning or thickening.In order to analyze the GOCE grav- ity gradient sensitivity on ice mass changes, a study on JI drainage basin ice mass change will be presented.
Gravity anomalies arising from changes in the surface topography can be approximated using a simple Bouguer correction.When applying it, the masses are approximated by a at, homogeneous and in nite plate, a socalled Bouguer plate, which has the density ρ and the height h (Sansó et al., 2013): If we assume that the surface consists of only ice and that this has a density of .g cm , the Bouguer correction can be expressed as The gravity anomaly changes (Fig. 6b) resulting from changes in ice mass topography (Fig. 6a), calculated using Eq. ( 14), show variations of ± mGal.This corresponds to an ice mass topography change of ± m.Here, the JI drainage basin is studied using Digital Elevation Models (DEMs).The spatial resolution of the data is .
• , i.e. approximately m in the north-south direction and around m in the west-east direction (Levinsen et al., 2013).These were developed by co-registering stereoscopic imagery from SPOT-5 to airborne and spaceborne laser data agreeing in time and space.The latter data are acquired with the National Aeronautics and Space Administration (NASA) laser scanner on-board the Airborne Topographic Mapper (ATM), and the Geoscience Laser Altimeter System (GLAS) on-board the Ice, Cloud and land Elevation Satellite (ICESat).The imagery has a high spatial resolution of approximately × m however elevation errors on the order of − m; the laser data has errors one order of magnitude smaller, however are constrained to satellite orbits and ight paths.Thus, by minimizing planimetric and elevation-dependent biases in be- tween the two datasets, any o sets between them were reduced, and high-resolution DEMs with low elevation errors were produced (Levinsen et al., 2013).A similar calculation is made with the Gravsoft TC program (Tscherning et al., 1992).This can calculate the direct topographic e ects of all masses above a reference level, assuming the density to be constant.The computation is based on two digital elevation models, a detailed and a coarse, which are used in the inner and outer zones, respectively.The two grids are assumed to have common boundaries, which is the case if the coarse grid has been constructed from the detailed grid by averaging.The integration of the terrain e ects is performed using the formulas for the gravitational e ects of a homogeneous rectangular prism.Depending on the geometry and accuracy, either exact formulas, spherical harmonic expansion, or a centered point mass approximation is used.The e ect of the ice topography change on the gravity anomaly is calculated using the TC program (Fig. 7).The results from using the Bouguer correction (Fig. 6b) and the TC program (Fig. 7), respectively, show a very good agreement in the prediction of gravity anomaly changes from a change in topography.However, more details in the gravity anomaly changes are present when using TC.The maximum gravity changes observed by the surface of the JI drainage basin range from − mGal for the period considered.The error when estimating gravity anomalies from GOCE gradient data using only the Tzz gradient component is mGal when the associated data error is mE.The gravity anomaly change occurring from November 2009 -June 2010 shows variations of ± mGal when using GOCE T ZZ gradients (Fig. 8).However, the small signal in the results is probably contaminated by the error in the GOCE T ZZ gradient dataset.The corresponding ice sheet residual mass change (Fig. 9) shows a variation of ± Mg.The largest changes occur in the center of the ice sheet, and not around the ice margin, thereby ignoring the rapidly changing outlet glaciers.
To further inspect the sensitivity of GOCE Tzz gradients on ice cover thickness, the TC program was used to compute the topographic e ects of the gravity gradient.The gravity gradient topographic e ects were calculated for two time intervals, 04.08.2007 (Fig. 10a) and 02.08.2008 (Fig. 10b).The maximum vertical gradient change at satellite altitude for the Jakobshavn area from August 2007 -August 2008 is .mE (Fig. 10c).
The gradients observed with GOCE have a minimum error of mE for the along-track component, and a combination of all four available components could potentially lower this to mE.Furthermore, we would need a longer period of GOCE Tzz gradient observations; a quali ed guess is ve years in order to observe this small signal.We expect maximally four years of observations to be available by mid-2014.Furthermore, starting from August 2012, the GOCE mission control team initiated the lowering of GOCE at a rate of approximately 300 meters per day.The objective is to lower the satellite from its nominal orbit of 255 km to a recommended height of 235.This will increase the accuracy and spatial resolution of the GOCE measurements, which may lead to the accuracy necessary to observe ice mass changes by the methods presented here.

Conclusions
The estimated gravity anomalies on the ground are very similar for the two methods used.However, both methods are a ected by the lack of several tracks in the monthly datasets.An idea for the calculation of the quantities year   reprocessed GOCE observations will be available by mid 2014.
The maximum gravity changes at the ground are between and mGal for the period considered.The gravity anomaly estimation error arising from the GOCE gradient data using only Tzz with an associated error of mE is mGal.Using more gradient components in the Gradiometer Reference Frame (GRF) would certainly reduce this error, presumably to − mGal.Thus, change detection may be made possible over two to three year periods.
There is hope for future GOCE data acquisitions to be able to further our understanding of mass changes using gravity gradients and gradient anomalies.We expect maximally four years of GOCE gradient observations to be available at the middle of 2014.Furthermore, lowering of the GOCE satellite down to an orbit of 235 km will increase the accuracy and spatial resolution of the measurements.This may then lead to the accuracy in order to observe ice mass changes using the GOCE gravity gradients.A new GOCE-type mission, with improved accuracy, will undoubtedly provide the possibility to detect areas of mass gain or losses over short time periods.
(a) GOCE Tzz gradients.(b) GOCE Tzz gradient anomalies when the contribution up to degree 36 is subtracted.

Fig. 1 .
Fig. 1.Reprocessed GOCE Tzz gradients (a) and GOCE Tzz gradient anomalies (b) when the contribution from Earth's Gravitational Model 2008 (EGM2008) up to spherical harmonic degree and order 36 is subtracted.The dataset consists of data from November 2009 to June 2010.The given units are Eötvös.

Fig. 2 .
Fig. 2. Monthly solutions for GOCE Tzz gradients when the contribution from the global model EGM2008 up to spherical harmonic degree and order 36 is subtracted.The given units are Eötvös.

( a )Fig. 3 .
Fig. 3. Comparison of predicted gravity anomaly residuals found using LSC (a) and the RPM (b) respectively.Area of Jakobshavn Isbrae glacier is outlined red in the di erence (c).The observations used are acquired in December 2009.

( a )
Surface elevation changes derived from a combination of airborne and spaceborne laser data with stereoscopic imagery.(b) Gravity changes.

Fig. 6 .
Fig. 6.Surface elevation changes over JI drainage basin (a), and the corresponding gravity changes (b) calculated using Eq.(14).The observation period is from 04.08.2007 to 02.08.2008.The shown corresponds to the region of Jakobshavn Isbrae glacier shown in Fig. 5.

Fig. 7 .
Fig. 7. Gravity changes calculated using the Gravsoft TC program from 04.08.2007 to 02.08.2008.The area shown corresponds to the region of Jakobshavn Isbrae glacier shown in Fig. 5.

Fig. 9 .
Fig. 9. Ice sheet residual mass change between November 2009 and June 2010 observed using GOCE T ZZ gradients.

Fig. 10 .
Fig. 10.JI Tzz gravity gradients at satellite altitude for the two di erent time intervals, (a) and (b), and the change between the two (c) due to the mass change.The given units are Eötvös.
by year and season by season was also tested, however did not lead to a higher signal in ice mass change.The maximum vertical gradient change at satellite altitude for the Jakobshavn Isbrae drainage basin area from August 2007 -August 2008 is .mE.The gradients observed by GOCE have a minimum error of mE for the along-track component, and a combination of all four available components could potentially lower this error to mE.Consequently, a 5-year observation period is necessary in order to observe such a small signal.Four years of Bereitgestellt von | Bibliothek des Wissenschaftsparks Albert Einstein Angemeldet Heruntergeladen am | 26.03.15 09:33

Table 1 .
GOCE TRF gradient data packages and related time interval covered.