Combining satellite altimetry and gravimetry data to improve Antarctic mass balance and gia estimates

COMBINING SATELLITE ALTIMETRY AND GRAVIMETRY DATA TO IMPROVE

ANTARCTIC MASS BALANCE AND GIA ESTIMATES

(1) _{Delft University of Technology, Stevinweg 1, 2628 CN Delft, NLD, b.c.gunter@tudelft.nl}

(2) _{Institute for Marine and Atmospheric research Utrecht, Utrecht University, P.O. Box 80.005, 3508 TA Utrecht, NLD} (3) _{School of Civil Engineering and Geosciences, Newcastle University, Newcastle upon Tyne NE1 7RU, UK}

(4) _{School of Geography and Environmental Studies, University of Tasmania, Hobart 7001, AUS}

(5) _{The University of Texas at Austin, Center for Space Research, 3925 West Braker Ln, Ste 200, Austin, TX 78759, USA}

ABSTRACT

This study explores an approach that simultaneously estimates Antarctic mass balance and glacial isostatic adjustment (GIA) through the combination of satellite gravity and altimetry data sets. The results improve upon previous efforts by incorporating reprocessed data sets over a longer period of time, and now include a firn densification model to convert the altimetry volume estimates into mass. When the GIA models created from the combination approach were compared to in-situ GPS ground station displacements, the vertical rates estimated showed good agreement after a systematic bias was removed from the computed GIA models. The new models suggest the potential for GIA uplift in the Amundsen Sea Sector, as well as the possible subsidence in large parts of East Antarctica.

1. INTRODUCTION

The most common approach for estimating present-day Antarctic GIA uplift rates is to utilize a numerical earth and ice history model to reconstruct the surface load time series since the last glacial maximum; however, the ice history and earth model parameters used for the reconstruction can have a wide range of plausible values, introducing uncertainty into the resulting GIA model. An alternative approach is to determine present day GIA rates through the use of space geodetic measurements, such as those from satellite altimetry and gravimetry missions, as well as in-situ GNSS surface displacements and climate data. This data-driven approach was shown to be feasible in an earlier study by Riva et al [1], but there now exist several new and updated data sets that have the potential to improve the accuracy of the GIA estimates produced. This includes recently reprocessed ICESat laser altimetry data (R633), an updated regional atmospheric climate model (RACMO2) and associated firn densification estimates, and an extended time series of Antarctic GNSS site displacements. This study will mostly focus on the contribution of each of these data sets towards the determination of Antarctic GIA, which is directly linked to the associated ice mass change estimates that are simultaneously computed as part of the inversion process.

2. METHODOLOGY

The methodology used to combine the altimetric and gravimetric data sets is adapted from earlier work by Riva et al [1], summarized here for convenience. In short, the technique relies on the fact that satellite altimetry measurements primarily observe surface processes, such as accumulation and ablation, whereas the mass change measurements from satellite gravimetry are sensitive to the mass change of both GIA and surface processes. By exploiting the difference in

density between ice/snow, ρsurf, and the solid earth, ρrock,

the following relationship can be established which

relates the vertical height rates of GIA, ℎ!"#, to the

mass, height, and density values for a given location.

ℎ

=

!!"#$!!!"#$   (1)

A 400 km Gaussian smoothing is applied to ensure the various components in Eq. 1 have the same spatial resolution, but this is only done after elements with equivalent resolution are first combined. For example, the multiplication of the surface density and ICESat height rates is done before applying the smoothing, since these two grids have approximately the same spatial resolution. How the surface and rock densities are treated will be covered in the next section.

3. DATA SETS

Several data sets are used to perform the combination, as well as validate the results. For this study, the mass change estimates were derived from the GRACE and the surface height trends derived from the ICESat mission. The properties of the surface, i.e., surface mass balance (SMB) and firn layer changes, were taken from Antarctic climate and firn densification models. The

solid earth densities were assumed to be 4000 kg/m3 _for

land and 3400 kg/m3 _{under the ice-shelves. The time}

period under investigation covers the entire ICESat mission period, from Feb 2003 to October 2010.

3.1. ICESat

The ICESat mission was the first Earth-orbiting laser

B.C. Gunter(1), O. Didova(1), R.E.M. Riva(1), M.R. van den Broeke(2), S.R.M. Ligtenberg(2), J.T.M. Lenaerts(2)_{, M. King}(3,4)_{, T. Urban}(5)

_____________________________________

Proc. ‘Earth Observation and Cryosphere Science Conf.’ Frascati, Italy, 13-16 November 2012 (ESA SP-712, May 2013)

altimeter and, while no longer operational, it was able to collect valuable information on the long-term surface height change of Antarctica over a period which overlaps the gravity data from GRACE. The surface height change trends used for this study are shown in Fig. 1, and were computed using the latest release of ICESat data (R633). Height changes were computed by comparing aggregate statistics of overlapping footprints within each grid cell (20km equal area blocks), similar to the approach developed by Slobbe et al [2], but using different editing criteria and no slope corrections. The ICESat laser shots are known to have a systematic bias in them that can introduce cm-level errors if neglected [3]. These campaign biases were estimated to be 1.5 cm/yr (with additional periodic variations) using a low-precipitation zone in East Antarctica, in the same line as Gunter et al [4] and Riva et al [1]. The total volume change for the dh/dt values is shown in Fig. 1a, and when integrated only over the grounded ice sheet, is

approximately -59 km3/yr.

3.2. GRACE

The GRACE mission has collected data on the time-variable nature of Earth’s gravity field since its launch in March 2002. A number of research centers produce monthly gravity field models, using different processing methodologies. Two sets of gravity models are examined in this study, one generated by the University of Texas Center for Space Research (CSR), and the Delft Mass Transport (DMT) models computed by TU-Delft [5]. The default resolution of the DMT fields extends to spherical harmonic degree and order 120 (120x120), while the CSR solutions extend to 60x60. To make the comparisons equivalent, the DMT

solutions used were regenerated using a

parameterization matching that of the CSR, i.e., to 60x60. Degree one coefficients were also added to both solutions using values generated by Swenson et al [6],

and the C20 values were replaced with those derived

from satellite laser ranging [7]. For the CSR models, the secular trends that are removed from select zonal coefficients were restored, as these rates are believed to mostly represent the effects of GIA. Both sets of solutions were based on RL04 GRACE data.

For the CSR solutions, a linear trend was estimated through each coefficient across the full time series of unconstrained monthly models. The trend was co-estimated with a bias, annual periodic, S2 (161 dy), and K2 (1362.7 dy) periodic terms. The long-term trend fit was then de-striped using an approach similar to that outlined by Swenson and Wahr [8]. The final result of this processing is shown in Fig. 1c. The long-term change for the DMT solutions were derived using that used by Siemes et al [9], which applies the anisotropic filtering method developed by Klees et al [10] after the long-term coefficient trend is estimated (along with bias, annual, and S2/K2 terms). These trends are shown

in Fig. 1d. Note that, in the final combination with the altimetry data, both GRACE solutions are smoothed with a 400km Gaussian spatial filter.

3.3. Firn densification model

In order to convert the volume changes derived from the ICESat data into mass, the density of the volume change needs to be known. There are many complex processes at work that complicate the determination of these densities, including regional variations in temperature, accumulation, and firn compaction. To account for these, a firn densification model (FDM) developed by Ligtenberg et al [11] is used that is forced at the surface with realistic 6-hourly climate output from the regional atmospheric climate model RACMO2 produced by Lenaerts et al [12]. This FDM model accounts for compaction of the firn over time, and is used in conjunction with the time-varying estimates of the total SMB from RACMO2 to estimate the mass change of the firn layer. Future references to the FDM imply the use of both the time-varying firn and SMB models. Fig. 1b shows the total surface height change as derived from the FDM model over the study period. It is important to note that the height changes of Fig. 1b only represents the surface height change of the firn, and do not reflect changes due to either the solid earth or ice dynamics.

When using the FDM to convert the ICESat volumes into mass, two basic assumptions were made to account for height differences that might exist between the altimetry measurements and the FDM. Because the FDM does not model ice dynamics, any negative differences between the ICESat and FDM surfaces were assumed to be the result of ice dynamics (glacier thinning), and the density assigned to this volume loss

was that of ice (917 kg/m3). Positive height differences

are attributed to an underestimation of SMB by RACMO2, and given a density closer to that of snow using a static density profile similar to that of Kaspers et al [13]. The only exception to this rule was for the region of the blocked Kamb ice stream, where the positive height change is assumed due to a build-up of ice (glacier thickening). These assumptions only deal with potential residual signal observed between ICESat and the FDM. The majority of the surface mass changes come directly from the SMB estimates derived from RACMO2. Mathematically, this is represented by a modification of Eq. 1,

ℎ

=

where !_!"#/!"#$ represents either the density of ice or

the static surface density map, depending on

3.4. GPS ground stations

The output from the combination represented by Eq. 2 is a spatial map of the vertical rates associated with GIA over Antarctica, hereafter referred to as the empirical rates. While this map can be compared against other GIA models in an attempt to assess its accuracy, an alternative approach is to compare the empirical rates with those observed by ground-based GPS stations. Fig. 2 shows the location of the GPS stations used for the comparisons to be shown later. The processing of the GPS displacements follows that of Thomas et al [14], and includes both seasonal and permanent stations. The trends in vertical displacement for these stations have been corrected for (modelled) elastic deformation effects. Furthermore, the time frames for the GPS trend analysis do not necessarily overlap with the GRACE and ICESat data sets. This is not as critical for the GPS rates, since GIA is assumed to evolve at near-constant rates over relatively short geologic time frames (e.g., decades).

Comparisons with the GPS data were done by computing the weighted root-mean-square of the residuals (WRMS) between the vertical empirical or modelled rates and those observed from the GPS stations,

!"#$ = !!⋅ ℎ!"#−ℎ!"#

2 !

!_! , (3)

where the weight,  ! = 1/!_!! _{is the inverse square of}

the uncertainty for GPS station i. These uncertainties ranged from < 0.3 mm (indicated by large circles in Fig. 2), 0.3 to 1.5 mm (medium circles), and > 1.5 mm (small circles).

The GPS data for three locations showed vertical rates that had large differences (> 5mm/yr), or were opposite in sign, to neighbouring GPS sites. The rates from these stations, designated by squares in Fig. 2, were considered as potential outliers and were therefore excluded from the WRMS calculations. In addition, the computed empirical rates showed a systematic bias with respect to the GPS rates that was on the order of 2-3 mm/yr. This bias was removed before the WRMS calculations were computed.

4. RESULTS

The combination approach was used to generate estimates of the GIA and ice mass change from the two GRACE solutions. The resulting empirical rates are shown in Fig. 2, overlaid with the corresponding GPS vertical rates. A summary of the WRMS fit of the solutions to the GPS rates can be found in Tab. 1. For comparison, the same WRMS calculation (including the bias removal and exclusion of potential outlier stations) was performed for three recent GIA models by Peltier (ICE-5G)[15], Ivins and James (IJ05)[16] and

Whitehouse et al (W12a)[17], also shown in Tab. 1. The geographical plots for these other models are also shown in Fig. 2. DMT-1b 1.41 CSR RL04 1.31 IJ05 1.47 ICE-5G 1.61 W12a 1.71

Table 1. The WRMS, in mm/yr, of the difference between the vertical rates of the GPS stations and the empirical

estimates, as well as several recent GIA models. 5. DISCUSSION

Several observations can be made when examining the results of the combinations. First, the empirical rates do not compare well with the GPS rates in the Antarctica Peninsula. This region is particularly dynamic, so there are many factors that may impact the comparison, such as potentially strong elastic effects on the GPS stations, the fact that ICESat is relatively data poor in this region, and the ability of GRACE to resolve the mass change of narrow North-South oriented features. Though not shown, the Antarctic Peninsula also has one the highest levels of uncertainty in the FDM.

In the Amundsen Sea (ASE) sector, the models indicate a level of uplift not typically predicted in this area by other models based on ice history reconstruction. There are several plausible reasons that might explain this signal. First, it could be that the gridded ICESat height change maps may be underestimating the total volume change in the ASE, so the increased surface mass predicted by the FDM gets interpreted as GIA in the combination. Alternatively, the FDM could be overestimating the surface mass balance in the region. It could also be that genuine GIA uplift is taking place in this area, as suggested by a recent study by Groh et al [18]; however, the uplift observed in the study was derived from only two seasonal GPS campaigns, so these results would need to be confirmed by additional long-term GPS measurements.

The widespread negative GIA signal observed in East Antarctica appears to correlate well with the coastal GPS sites that have high precision, and that also predict a small level of subsidence. Similar patterns are also observed in the W12a and IJ05 models.

Looking at the WRMS values in Tab. 1, the bias-corrected GIA solutions are, at a minimum, comparable to the other GIA models derived from more traditional methods, i.e., ICE-5G, IJ05 and W12a. The solution involving the CSR GRACE data had the best fit to the GPS rates, which provides some confidence to the methodology and assumptions used in the processing. Although, there are large parts of Antarctica that currently do not have long-term GPS records, including some areas where the estimated signal is large.

Recently installed permanent stations for some of these areas will help clarify some of these discrepancies, but this will take several years to develop.

6. CONCLUSIONS

This study revisited the approach of estimating present-day GIA and ice mass change using a combination of satellite altimetry and gravimetry. An updated and extended ICESat surface height change map was combined with two different GRACE solutions, along with an advanced regional atmospheric climate model and associated firn densification model. The empirical GIA rates generated from this approach showed good agreement to an independent set of GPS-derived vertical rates when a systematic bias term was removed. The exact cause of this bias offset is being investigated, but

it is likely related to the choice of degree-1 and C20

coefficients, as well as the ICESat campaign biases. The empirically derived GIA rates show some noticeable differences to other recent GIA models derived using the more traditional ice history reconstruction approach, such as in the Amundsen Sea sector, the Antarctic Peninsula, and central East Antarctica. Many of these areas are without long-term GPS records, so the vertical rates cannot currently be verified. Nonetheless, the results from the combination approach are encouraging, and demonstrate that the technique has the potential to reduce the uncertainty surrounding both Antarctic GIA and ice mass change estimates.

7. REFERENCES

1. Riva, R.E.M., Gunter, B.C. Urban, T., Vermeersen, L.L.A., Lindenbergh, R., Helsen, M., Bamber, J., van de Wal, R., van den Broeke, M., and Schutz B. (2009). Glacial Isostatic Adjustment over Antarctica from combined GRACE and ICESat satellite data, Earth and Planetary Science Letters, 288, pp. 516–523.

2. Slobbe, D., Lindenbergh, R., Ditmar, P. (2008). Estimation of volume change rates of Greenland's ice sheet from ICESat data using overlapping footprints. Remote Sensing of the Environment, 112, pp 4204–4213. 3. Urban, T., Schutz B. (2005). ICESat sea level comparisons. Geophysical Research Letters, 32 (L23S10). doi:10.1029/10.1029/2005GL024306

4. Gunter, B.C., Riva, R.E.M., Urban, T., Harpold, R., Schutz, B., Nagel, P., Helsen, M. Evaluation of GRACE and ICESat Mass Change Estimates Over Antarctica. Gravity, Geoid and Earth Observation International Association of Geodesy Symposia Volume 135, 2010, pp 563-569.

5. Liu, X., Ditmar, P., Siemes, C., Slobbe, D. C., Revtova, E., Klees, R., Riva, R., and Zhao, Q. (2010). DEOS mass transport model (DMT-1) based on GRACE satellite data: methodology and validation. Geophysical

Journal International, 181 (2), 769-788.

6. Swenson, S., Chambers, D., and Wahr, J. (2007). Estimating geocenter variations from a combination of GRACE and ocean model output. Journal of Geophysical Research, 113.

7. Cheng, M. and Tapley, B.D. (2004). Variations in the Earth's oblateness during the past 28 years. Journal of Geophysical Research, 109 (B9).

8. Swenson S, Wahr J (2006) Post-processing removal of correlated errors in GRACE data. Journal of Geophysical Researtch, 33(L08402). doi:10.1029/2005GL025285 9. Siemes, C., Ditmar, P., Riva, R.E.M., Slobbe, D.C., Liu,

X., Hashemi Farahani, H. (2012). Estimation of mass change trends in the Earth’s system on the basis of GRACE satellite data, with application to Greenland. Journal of Geodesy, DOI:10.1007/s00190-012-0580-5 10. Klees, R., Revtova, E.A., Gunter, B.C., Ditmar, P.,

Oudman, E., Winsemius, H.C., Savenije, H.H.G. (2008). The design of an optimal filter for monthly GRACE gravity models. Geophysical Journal International, 175 (2), 417-432.

11. Ligtenberg, S.R.M., Helsen, M.M., van den Broeke, M. R. (2011). An improved semi-empirical model for the densification of Antarctic firn. The Cryosphere, 5, 809-819, doi: 10.5194/tc-5-809-2011.

12. Lenaerts, J.T.M., van den Broeke, M.R., van de Berg, W.J. van Meijgaard, E., and Kuipers Munneke, P. (2012). A new, high-resolution surface mass balance map of Antarctica (1979–2010) based on regional atmospheric climate modeling. Geophysical Research Letters 39, L04501, doi:10.1029/2011GL050713.

13. Kaspers, K.A., van de Wal. R., van den Broeke, M.R., Schwander, J., van Lipzig, N.P.M., Brenninkmeijer, C.A.M. (2004). Model calculations of the age of firn air across the Antarctic continent. Atmos Chem Phys 4:1365–1380.

14. Thomas, I.D., King, M.A., Bentley, M.J., Whitehouse, P.L., Penna, N.T., Williams, S.D.P., Riva, R.E.M., Lavallee, D.A., Clarke, P.J., King, E.C., Hindmarsh R.C.A., and Koivula, H. (2011). Widespread low rates of Antarctic glacial isostatic adjustment revealed by GPS observations. Geophysical Research Letters, 38: L22302. doi:10.1029/2011GL049277.

15. Peltier, W.R., 2004. Global glacial isostasy and the surface of the ice-age Earth: the ICE-5G (VM2) Model and GRACE. Annu. Rev. Earth Planet. Sci. 32, 111149. 16. Ivins, E.R., James, T.S., 2005. Antarctic glacial isostatic

adjustment: a new assessment. Ant. Sci. 17 (4), 537549. 17. Whitehouse, P., Bentley, M., Milne, G., King, M.,

Thomas, I. (2012). A new glacial isostatic adjustment model for Antarctica: calibrated and tested using observations of relative sea-level change and present-day uplift rates. Geophysical Journal International. doi: 10.1111/j.1365-246X.2012.05557.x

18. Groh, A., Ewert, H., Scheinert, M., Fritsche, M., Rulke, A., Richter, A., Rosenau, R., Dietrich, R. (2012). An Investigation of Glacial Isostatic Adjustment over the Amundsen Sea Sector, West Antarctica. Global and Planetary Change. doi: 10.1016/j.gloplacha.2012.08.001

a)

b)

c)

d)

Figure 1. Long-term trends for: a) surface height change observed from ICESat, b) surface height change due to surface and firn layer processes as simulated by the FDM, c) the mass change in equivalent water height from the CSR

GRACE solutions, and d) the mass change in equivalent water height from the DMT-1b GRACE solutions

cm/yr

cm/yr cm/yr

a)

b)

c) d) e)

Figure 2. Annual vertical displacement from the GIA models: a) CSR RL04 GRACE and, b) DMT-1b GRACE data, c) ICE-5G, d) IJ05, and e) W12a. Each map is overlaid with the location of the GPS stations used for validation, with the estimated precision of the GPS rates indicated by the symbol size (large is < 0.3 mm/yr, medium is between 0.3 and 1.5

mm/yr, and small is > 1.5 mm/yr). Square symbols indicate potential outlier stations that were excluded from the WRMS calculations.