Introduction
The present-day and future contribution of the Antarctic ice sheet to sea-level rise has considerable uncertainty associated with it (Reference Lemke and SolomonLemke and others, 2007). To reduce this uncertainty, models of the entire ice-sheet–shelf system are needed and the thickness of the floating ice shelves around Antarctica provides an important boundary condition for these and for sub-ice-shelf cavity models. Ice-shelf thickness is also required, close to the grounding line, to calculate the mass balance of the ice sheet using the mass budget approach (Reference RignotRignot and others, 2008). In this method, the velocity of the ice crossing a gate close to the grounding line is multiplied by the ice-shelf thickness at the gate to estimate the flux of ice leaving the ice sheet. This can be compared to mass input, in the form of accumulation, to determine ice-sheet mass balance. In this calculation, the ice thickness at the gate, along with accumulation, are the two dominant uncertainties and account for ∼90% of the estimated error (Reference RignotRignot and others, 2008). Ice-shelf thickness is also required in calculations of sub-ice-shelf mass balance (Reference Joughin and PadmanJoughin and Padman, 2003; Reference Wen, Wang, Wang, Jezek, Liu and AllisonWen and others, 2010).
A warming climate implies a greater vulnerability of the ice shelves to collapse; in recent decades a number of shelves on the Antarctic Peninsula have disintegrated. This can have a significant impact on the flow of glaciers that were previously buttressed by the ice shelf (Reference Rott, Skvarca and NaglerRott and others, 1996; Reference De Angelis and SkvarcaDe Angelis and Skvarca, 2003; Reference Rignot, Casassa, Gogineni, Krabill, Rivera and ThomasRignot and others, 2004; Reference Scambos, Bohlander, Shuman and SkvarcaScambos and others, 2004). There is therefore considerable interest in the stability of these ice shelves, and accurate estimates of their thickness are essential for modelling their behaviour. For example, recent inverse modelling studies have derived ice-shelf rheology from estimates of ice-shelf velocity and thickness (Reference Khazendar, Rignot and LarourKhazendar and others, 2007, Reference Khazendar, Rignot and Larour2009, in press).
Ice-shelf thickness can be measured at point locations in situ or determined from radio-echo sounding (RES). These data, while usually having high accuracy, are generally severely limited by the spatial coverage between tracks, which can be tens to hundreds of kilometres (Reference Lythe and VaughanLythe and others, 2001). Ice-shelf thickness can also be estimated from surface elevation if it is assumed that the ice shelf is freely floating in hydrostatic equilibrium (Reference Bamber and BentleyBamber and Bentley, 1994). The increased accuracy and availability of satellite altimetry in recent decades, newly available, spatially resolved, physically based estimates of firn densities and improved estimates of the geoid for the Southern Ocean make accurate retrieval of ice-shelf thickness from surface elevation a viable prospect. In this paper, we calculate ice-shelf thicknesses for all Antarctic ice shelves with areas of >10 km2. The work builds on the case study of the Larsen C ice shelf by the authors (Reference Griggs and BamberGriggs and Bamber, 2009a), extending our method to the remaining ice shelves of sufficient size. In the limited regions where ice-shelf thickness grids created from densely spaced RES or in situ measurements exist, those data are likely of higher accuracy. In this study, we do not combine the two, but instead use the terrestrially derived observations to provide an independent validation of our thickness estimates. The results presented here supplement the small number of field-based ice thickness estimates available. In particular, our results provide a consistent, continent-wide dataset. We begin by describing the data available and then the methods used and their limitations. Next we describe the gridded surface elevations, their validation and conversion to ice thickness and finally compare the results with airborne data.
Data and Data Processing
Surface elevations are derived mainly from the European Remote-sensing Satellite (ERS-1) radar altimeter (RA). In April 1994, ERS-1 was placed into a geodetic phase of two long repeat cycles of 168 days, equivalent to a single 336 day orbit (phases E and F). This gives an across-track spacing of 2 km at 70° S, an along-track spacing of 335 m with an effective spatial resolution of ∼4 km (Reference Lingle, Brenner and ZwallyLingle and others, 1990). These data are used so that coverage is maximized particularly on the smaller ice shelves which may be missed entirely with a shorter orbit cycle. For example, the operational Ice, Cloud and land Elevation Satellite (ICESat) orbit cycle results in ∼10 times coarser across-track spacing (Reference Schutz, Zwally, Shuman, Hancock and DiMarzioSchutz and others, 2005). High spatial resolution is particularly important close to the grounding line where surface slopes are largest. Studies have shown a slope-dependent bias between the ERS-1 RA data and those from more recent satellite laser altimetry, but this amounts to centimetres over the relatively flat ice shelves (Reference Bamber and Gomez-DansBamber and Gomez-Dans, 2005). The timing of the geodetic phase of ERS-1 means that the ice thicknesses presented here have an effective timestamp of January 1995.
The RA data are the same as those used to derive a widely used 5 km digital elevation model (DEM) of Antarctica (Reference Bamber and BindschadlerBamber and Bindschadler, 1997) and, in combination with ICESat data, to derive a more recent 1 km DEM (Reference Bamber, Gomez-Dans and GriggsBamber and others, 2009; Reference Griggs and BamberGriggs and Bamber, 2009b). The reason that these DEMs were not used directly to derive ice-shelf thickness is primarily related to data coverage and interpolation close to the grounding line (Reference Bamber, Gomez-Dans and GriggsBamber and others, 2009). In this region, there is often a break in slope moving from the floating to grounded ice (Reference Fricker and PadmanFricker and Padman, 2006), which may result in the ERS-1 RA being unable to track the surface (known as loss of lock) and/or off-ranging (e.g. Reference Bamber, Gomez-Dans and GriggsBamber and others, 2009, fig. 12). This can create gaps in coverage on the order of 10 km close to the grounding line (see Fig. 1a, close to Roosevelt Island). Interpolation between the higher grounded ice and the floating ice shelf can result in a positive bias across this data-sparse zone. To avoid this problem, it is better to interpolate across the floating portion only, excluding points inland of the grounding line. The ERS-1 RA data were corrected for slope-induced error and the inverse barometer effect, and data were removed where radar lock was lost (Reference Bamber and BindschadlerBamber and Bindschadler, 1997). No tide model was previously applied to these data.
Elevation (m) of available ERS-1 data points with respect to the geoid on (a) Ross Ice Shelf, (b) Jelbart and Fimbul Ice Shelves and (c) Porpoise Ice Shelf (locations marked by red (a), green (b) and blue (c) boxes on Antarctic map). The elevations are shown plotted over the MOA image of Antarctica, with black indicating open ocean (Reference Scambos, Haran, Fahnestock, Painter and BohlanderScambos and others, 2007). The thick black lines show the position of the MOA grounding and coast lines.
Large parts of the Ross and Filchner–Ronne Ice Shelves are south of the latitudinal limit of the ERS-1 RA. In these regions, we use surface elevations derived from ICESat data (H.J. Zwally and others, http://nsidc.org/data/gla12.html). The Geoscience Laser Altimeter System (GLAS) instrument on board ICESat is a laser altimeter which flew from January 2003 until October 2009, with science data available from February 2003. The lasers are switched on two to three times per year, evenly spaced throughout the year. The ICESat orbit gives an across-track spacing of ∼20 km at 70° S, an along-track spacing of 172 m and a ∼65 m diameter footprint on the surface (Reference Schutz, Zwally, Shuman, Hancock and DiMarzioSchutz and others, 2005). GLAS has a vertical accuracy on the order of 20 cm at the slopes typical of ice shelves (Reference Brenner, DiMarzio and ZwallyBrenner and others, 2007) and does not suffer from loss of lock. It cannot, however, penetrate optically thick cloud, and thus we use multiple laser operation periods giving repeat observations of the same ground tracks, lowering random noise and filling in areas of cloud cover encountered in individual operation periods. The ICESat product used was the level 2 Antarctic and Greenland Ice Sheet Altimetry Data product (GLA12) downloaded from the US National Snow and Ice Data Center (NSIDC), and all data were Release 428 because the later release was not available for all campaigns at the time of processing. Laser campaigns 1a (20 February–21 March 2003) through to 3j (17 February–21 March 2008) were used and were extracted and transformed to the World Geodetic System 1984 (WGS84) ellipsoid using the software provided by the NSIDC. A correction was applied to the data to account for saturation of the laser over ice (Zwally and others, http://nsidc.org/data/gla12.html), and the applied tide model was removed. ICESat data had geophysical quality-assurance filters applied to remove data that may contain residual cloud or other artefacts that affect the elevation estimate (Reference Bamber, Gomez-Dans and GriggsBamber and others, 2009).
To prevent data from grounded ice affecting the interpolation of floating ice, all altimeter returns for the former were removed. The grounding line was determined by combining the Moderate Resolution Imaging Spectroradiometer (MODIS) Mosaic of Antarctica (MOA) (Reference Scambos, Haran, Fahnestock, Painter and BohlanderScambos and others, 2007; T. Haran and others, http://nsidc.org/data/nsidc-0280.html) with estimates determined from interferometric synthetic aperture radar (InSAR) data (personal communication from E. Rignot, 2010). They were combined by resampling the MOA and InSAR grounding lines onto a grid of 100 m resolution. Where differences existed between the two estimates, we chose the values farthest inland to ensure that no floating data were excluded. The ice-shelf edge was determined from the MOA coastline, resampled onto a 100 m resolution grid. Ice rises and grounded islands were also identified using MOA. A number of islands were found to be missing from MOA, including Sherman Island on the Abbot Ice Shelf and Smyley Island and a neighbour on George VI Ice Shelf (Fig. 3, further below). Ice shelves of less than ∼10 km2 were ignored as they are unlikely to contain sufficient data. Elevations of <5 m with respect to the EIGEN-GL04C geoid (Reference FörsteFörste and others, 2008) were also removed to avoid contamination close to the calving front.
The influence of ocean tides was removed from both the laser and radar altimeter data using the TPXO6.2 tide model (Reference Egbert and ErofeevaEgbert and Erofeeva, 2002). This is a global model and represents a least-squares best fit to the Laplacian tide equations and TOPEX/Poseidon and Jason satellite data. The model provided the eight primary tide components as well as two long-period components at 0.25° spatial resolution. Ocean-loading tides and solid-Earth tides were also included at this stage of processing. This model was chosen as it has been shown to have superior accuracy around Antarctica (Reference King, Penna, Clarke and KingKing and others, 2005). It uses a predetermined, inbuilt grounding line and, inland of that, calculates no ocean tide. The tide-model grounding line was significantly different from that determined by MOA/InSAR, in particular for the ice shelves along the coast of Dronning Maud Land (Fig. 3, further below) where it lay roughly at the midpoint of the ice shelf. In such cases, we extrapolated the closest floating point in the tide model to the grounding line we used.
Method
Ice thickness was inferred from surface elevation using the principle of hydrostatic equilibrium. If ice is in hydrostatic equilibrium, its thickness can be determined as
where H i is the equivalent ice thickness, i.e. the thickness that the ice shelf would be if all the ice were at the meteoric ice density, ρ i. Z is the actual ice thickness (i.e. including a layer of variable-density firn), e is the ice-shelf elevation above mean sea level (the freeboard), ρw is the density of the water column under the ice shelf and δ is the firn density correction, i.e. the difference between the actual depth of the firn layer and the depth that the firn would be if it were at the density of meteoric ice.
The assumption of hydrostatic equilibrium is only valid if the ice shelf is freely floating. This is typically only the case at a distance of several ice thicknesses from the grounding line. The first point of hydrostatic equilibrium (defined as point H in Reference Fricker and PadmanFricker and Padman, 2006) is at the outer boundary of the first fringe of tidal displacement from double-difference InSAR (Reference RignotRignot, 1996). The distance between the first point of hydrostatic equilibrium and the grounding line depends on a number of factors in addition to thickness, including ice rheology, dynamics and bedrock topography. Estimates of the distance for two locations on Institute Ice Stream flowing into the Filchner–Ronne Ice Shelf vary between ∼4.2 and 6.4 km (Reference Fricker and PadmanFricker and Padman, 2006) and 9 km on Thwaites Glacier (Reference RignotRignot, 2001). On the Ross Ice Shelf, a mean width of 3.2 km, with a standard deviation of 2.6 km, and two cases of grounding zone widths greater than 10 km have been reported (Reference Brunt, Fricker, Padman, Scambos and O’NeelBrunt and others, 2010). Thus, in general, hydrostatic equilibrium may not be valid up to ∼10 km from the grounding line. Equation (1) employs a single ice density that is constant throughout the column below the firn layer. This assumption is unlikely to be robust in areas of crevassing and also when marine ice is present (discussed in more detail below).
Kriging
We interpolated the data onto a regular polar stereographic grid with a standard parallel of 71° S and a central meridian of 0°. We used a grid spacing of 1 km to maximize the number of gridboxes containing data while preserving smaller-scale features on the ice surface. We used ordinary kriging for the interpolation based on the GSLIB software (Reference Deutsch and JournelDeutsch and Journel, 1998). We chose this approach so that the interpolation was related to the variance in thickness of each ice shelf. The kriging method uses the spatial variability of the surface, mathematically described in the form of a variogram, to determine the interpolated data point, preserving the actual data value at locations where an observation is present. Variogram shape was determined separately for each ice shelf, as data density and surface roughness (i.e. the statistical properties of the input) are not constant between shelves. Variograms were determined from the ERS-1 RA data available on each ice shelf and fitted to an exponential using a nonlinear least-squares algorithm for lags up to 50 km. The parameters of the exponential fit are listed in Table 1. In the case of three ice shelves, Progress, Shackleton and Vanderford, the exponential was fitted up to a lag of 100 km, and a maximum search radius of 100 km was used due to poorer data coverage. In the case of four of the smallest ice shelves, Tucker in Victoria Land, Hull on the Ruppert Coast, Bryan and Loubet on the Peninsula, there were insufficient data over a wide enough area of the ice shelves to calculate representative variograms, so no result is derived for these four shelves. A nugget of 1 m was applied to all variograms. A maximum search radius of 50 km was set and up to 50 nearby data points were considering in the kriging.
Shelf locations and sizes
The Filchner–Ronne and Ross Ice Shelves both extend south of 81.5° S (the latitudinal limit of ERS-1), so elevations from ICESat were needed to produce a thickness estimate for the complete ice shelf. We interpolated the ICESat and ERS-1 datasets separately using the same variogram. The two surfaces were then merged by weighting them using a Hermite basis function over a distance of 20 km at the southern limit of the ERS-1 data. Only in the 20 km region at the southern limit of ERS-1 are both datasets used. Examination of the two datasets close to the southern limit of ERS-1 shows no consistent elevation bias, with differences of 0.29 ± 1.56 m in the 20 km region. The spatial resolution of the resulting merged elevation is lower in the ICESat-only region due to the considerably larger across-track spacing of ICESat compared with the geodetic-phase ERS-1 RA data.
Geoid correction
The surface elevation grids are relative to the WGS84 ellipsoid and need to be transformed to elevations relative to mean sea level. Recent geoids using Gravity Recovery and Climate Experiment (GRACE) satellite data give much improved accuracy compared to older models. We used the EIGEN-GL04C geoid in the mean tide system (Reference FörsteFörste and others, 2008). This geoid is calculated using GRACE satellite data supplemented with surface gravity data from altimetry and gravimetry. Comparisons with independent GPS-derived geoid heights gave a root-mean-square (RMS) difference of up to 0.36 m and an error budget of 0.23 m (Reference FörsteFörste and others, 2008).
Mean dynamic topography
The mean dynamic topography is the height of the mean ocean surface relative to the geoid as a result of non-tidal effects. It is the response of the ocean to planetary rotation, hydrographic variability and wind. It can be estimated by comparing the geoid height with tide-free altimetric measurements of the ocean surface height (Reference Andersen and KnudsenAndersen and Knudsen, 2009). Over ice shelves, this is not possible, and thus there are no estimates of mean dynamic topography for the ice shelves. In the case of smaller ice shelves, the mean dynamic topography should be close to that at the ice front. For the larger ice shelves such as the Ross, Filchner–Ronne and Amery, sub-ice-shelf currents exist, which are substantially different from the general ocean circulation (Reference Nicholls, Makinson and ØsterhusNicholls and others, 2004). Additionally, there is substantial meso-scale variability in mean dynamic topography in the Southern Ocean, and values vary between 0 and −1 m over scales of 100 km. Thus, extrapolating a mean dynamic topography from open ocean across the whole of the larger ice shelves is likely to induce errors of a similar magnitude to, or larger than, the correction. For this reason, we do not apply a correction for the mean dynamic topography to the Filchner–Ronne, Ross and Amery ice shelves. For the smaller ice shelves, we use mean dynamic topography from DTU10MDT (Reference Andersen and KnudsenAndersen and Knudsen, 2009) with respect to the EIGEN-GL04C geoid. We calculate the mean value at the ice-shelf front and extrapolate to the entire ice shelf. Values range between −1.1 and −2.0 m. On the larger ice shelves, we estimate the lack of a correction adds a mean error of 1 m (9 m thickness), with a maximum of 2 m, to the freeboard estimate (Reference Andersen and KnudsenAndersen and Knudsen, 2009).
Firn correction
The firn correction is required to account for the presence of a variable-density layer above the fully compacted meteoric ice, which makes up the majority of the ice-shelf thickness. It is defined by δ in Equation (1) and has units of length, as it represents the difference in thickness of a column of ice with density ρ i and the larger thickness of the actual combined firn–ice column. The value of δ can range between 0, for blue ice areas on, for example, the Amery Ice Shelf (Reference Van de Berg, van den Broeke, Reijmer and van MeijgaardVan de Berg and others, 2005), and >20 m (∼180 m thickness) for areas of convergent flow onto, for example, the Ross Ice Shelf (Reference Bamber and BentleyBamber and Bentley, 1994). The depth and density of the firn layer cannot be detected by satellite, and only a few in situ observations exist. Here we used modelled firn-density corrections as adopted in an earlier mass budget study (Reference RignotRignot and others, 2008). A firn densification model was forced by the output from a regional climate model (Reference Van de Berg, van den Broeke, Reijmer and van MeijgaardVan de Berg and others, 2005, Reference Van de Berg, van den Broeke, Reijmer and van Meijgaard2006), which was, in turn, forced at its lateral boundaries by European Centre for Medium-Range Weather Forecasts re-analysis and operational analysis data. The firn densification model accounts for temperature, accumulation and wind-speed variability (Reference HelsenHelsen and others, 2008; Reference Van den BroekeVan den Broeke, 2008) using 25 year averages (1980–2004) and was run at a resolution of ∼55 km. Spatially, the firn density correction shows large variability, particularly close to the grounding line, but values typically vary from ∼13 to 19 m. The role of katabatic winds on densification rates is important, so including wind-speed variability in the model is key. Previous studies have generally assumed a constant firn-density correction (Reference Lythe and VaughanLythe and others, 2001) due to the absence of information on its spatial variability. Uncertainty in the firn density correction can therefore dominate the error budget in the ice thickness estimation and it increases, in percentage terms, with decreasing ice thickness. The modelled density has been compared to in situ observations from firn cores and agrees well, except in one case in the upstream region of Whillans Ice Stream where horizontal compression, which is not modelled, causes rapid compaction of the firn layer (Reference Van den BroekeVan den Broeke, 2008). Spatial variability of the modelled firn is high, amounting to a standard deviation of 24% over the Larsen C ice shelf.
Densities
We use values of 1027 and 917 kg m−3 for the densities of sea water and meteoric ice respectively. 917 kg m −3 is a standard value for the density of meteoric ice (Reference Benn and EvansBenn and Evans, 1998), and while impurities will cause variability about that value, values only vary by around ±15 kg m−3. Values between 912 and 922 kg m−3 were obtained for the Amery Ice Shelf when using an ice density model constrained by borehole measurements (Reference Fricker, Popov, Allison and YoungFricker and others, 2001). Measurements of average ice densities are not widespread, so the meteoric ice density is used throughout this study, and variability of ±15 kg m−3 is noted as one of the errors. 1027 kg m−3 is the mean global sea-water density. Again, measurements are sparse, particularly for sea-water densities under ice shelves, so assuming a variable density would introduce uncertainty. A density of 1024 kg m−3 has been inferred for the Ross Ice Shelf by using elevations and thickness and assuming that ice density is known (Reference Bamber and BentleyBamber and Bentley, 1994). A value of 1029 kg m−3 has been used for the Amery Ice Shelf from an ocean measurement near the ice front (Reference Fricker, Popov, Allison and YoungFricker and others, 2001). Thus an error of ±15 kg m−3 is assumed.
Other error sources: marine ice, horizontal compression, melt
A number of other factors can affect the ice-shelf thickness, but including them in the calculations induces unknown errors into the estimates. These other factors are discussed in turn, listing their cause and possible magnitude.
Marine ice is widespread beneath the ice shelves and is formed when water at the surface freezing point melts the underside of the ice shelf. The resulting meltwater sets up a thermohaline circulation where cold, buoyant water comes into contact with the ice-shelf bottom and freezes onto it to form porous marine ice with trapped interstitial sea water. Marine ice can be formed in stripes aligned with fast flow from outlet glaciers and has been shown to extend from close to the grounding line right to the ice front on the Larsen C ice shelf (Reference Holland, Corr, Vaughan, Jenkins and SkvarcaHolland and others, 2009). The transition from meteoric ice to denser marine ice can often form the last return recorded by RES, so care must be taken to ensure that airborne ice thickness measurements capture the full ice thickness and not just the meteoric layer. Comparison of RES data with thickness derived from hydrostatic equilibrium on the Amery Ice Shelf has been used to infer marine ice thicknesses (Reference Fricker, Popov, Allison and YoungFricker and others, 2001); however, this requires accurate knowledge of the firn, ice and water densities. Borehole measurements on the Amery Ice Shelf have measured a maximum marine ice density of 938 kg m−3 (Reference Craven, Allison, Fricker and WarnerCraven and others, 2009). If we assume that marine ice accounts for half of the total ice-shelf thickness and that the firn-corrected surface elevation is 90 m, the error induced by not including marine ice would be 5%. This is an upper bound, as, locally, marine ice comprising 50% of the ice-shelf thickness is larger than all studies bar one suggest (Reference Zotikov, Zagorodnov and RaikovskyZotikov and others, 1980; Reference OerterOerter and others, 1992; Reference Craven, Allison, Fricker and WarnerCraven and others, 2009). Marine ice cannot be easily detected from the surface, and its distribution, density and thicknesses, none of which are comprehensively known for ice shelves other than the Amery, would be required for its effects to be included in our calculations. Thus it is considered as an error term of 5% of the thickness derived in the present study.
Two effects are not included in the steady-state model which derives the firn density correction. Firstly, horizontal compaction of the ice is not included but has been shown to be important (Reference Van den BroekeVan den Broeke, 2008). This effect will be most dominant in the onset regions of fast-flowing ice streams (e.g. Reference Bamber and BentleyBamber and Bentley, 1994), but, for the floating ice shelves considered here, may affect the ice shelf in regions of fast-flowing outflow such as where Institute Ice Stream enters the Filchner–Ronne Ice Shelf. The other process which is not modelled is the effect of melt on the firn correction. Melt is widespread on the ice shelves (Reference Van de Berg, van den Broeke, Reijmer and van MeijgaardVan de Berg and others, 2005); however, almost all of that meltwater refreezes within the snowpack. The refreezing process is not modelled in the firn compaction model, so modelled estimates of melt cannot be incorporated into our firn density correction. Attempts to correlate melt from the model with firn density corrections calculated from independent and coincident estimates of ice thickness and surface elevation showed no significant correlation. Melt is most pronounced on the Antarctic Peninsula and here the firn density corrections were lowered by 2.8 m, based on the difference between the average modelled firn-density correction on Larsen C and in situ measurements across West Antarctica (Reference Griggs and BamberGriggs and Bamber, 2009a).
Results and Validation
Table 1 lists the location of each ice shelf ordered by area. The number of ERS-1 data points and percentage of gridboxes containing ERS-1 RA data are also shown. The kriging variogram is defined by the positive variance contribution value, c, and the effective range, a, as listed in the rightmost column of Table 1. Figure 1 shows a typical distribution of ERS-1 surface elevation data from three ice shelves: the Ross Ice Shelf, representative of the largest, most southerly ice shelves (Fig. 1a); the Jelbart and Fimbul Ice Shelves, representative of medium-sized, mid-latitude ice shelves (Fig. 1b); and the Porpoise Ice Shelf, representative of the distributions of ERS-1 data seen on the smallest ice shelves and those farthest north (Fig. 1c). The effect of ERS-1 losing lock in break-in slope regions close to the grounding line is illustrated clearly in Figure 1a, particularly around Roosevelt Island and Siple Dome where data gaps are seen in bands ∼10 km wide aligned with the grounding line. In Figure 1b, a data gap is seen on the Fimbul Ice Shelf around 2190 km north. This is due to the switch from ice mode to ocean mode of the ERS-1 RA (Reference FrancisFrancis and others, 1991).
Comparison of elevation with ICESat
We compare the derived surface elevation grids (DEMs) with independent surface elevation data from ICESat to determine the accuracy of the gridded product. The ICESat data were recorded ∼10 years after the ERS-1 RA data, so differences which may be due to real changes over that time period were accounted for by including an elevation-rate (dH/dt) correction in the comparison. The dH/dt correction was from trends in ERS-1 and ERS-2 data between 1992 and 2001 (Reference ZwallyZwally and others, 2005). We used bilinear interpolation to determine the surface elevation at the exact location of the ICESat observation from the gridded product, and differences are calculated in the sense of ICESat data-point minus interpolated DEM value.
Table 2 lists the mean difference and RMS difference between the two measures of surface elevation for all ice shelves both with and without correction for the dH/dt signal. In most cases, mean differences are negative, implying that the ERS-1 derived elevations are higher than the more recent ICESat elevations. However, dH/dt correction does not fully explain these differences. In some cases such as the Filchner–Ronne, the dH/dt correction makes the mean differences slightly larger.
Mean and RMS difference between ICESat surface elevations and gridded DEM values from ERS-1. All differences are in metres
The comparison with the ICESat data shows that the ERS-1 surfaces are sufficiently accurate. We find mean differences in general below 5 m, and the RMS of those differences generally below 15 m. However, in a number of cases, particularly where the ERS-1 or ICESat data are not evenly distributed across the ice shelf, the observed errors are significantly higher. Figure 2 shows the distribution of differences on three representative ice shelves.
Differences between gridded elevation dataset and ICESat elevations (m) on (a) Filchner–Ronne Ice Shelf, (b) Brunt and Riiser-Larsen Ice Shelves and (c) Getz Ice Shelf. Locations marked by red (a), green (b) and blue (c) boxes on Antarctic map inset in (a). The differences are shown plotted over the MOA image of Antarctica, with black indicating open ocean (Reference Scambos, Haran, Fahnestock, Painter and BohlanderScambos and others, 2007). The thick black lines show the position of the MOA grounding and coast lines, and the green square indicates the colour of zero difference. Note that (a) has a smaller range than (b) and (c).
The Filcher–Ronne Ice Shelf (Fig. 2a) shows results very similar to those on the rest of the largest ice shelves. In general, the differences are close to zero, with larger negative differences (ERS-1 interpolated surface higher than ICESat) close to the grounding line where ERS-1 loses lock. The mean and RMS differences increase by up to an order of magnitude within 10 km of the grounding line, although mean differences still represent typically <10% of the total thickness, being, in general, <10 m (columns 4–8 of Table 2). Large differences are also seen close to the ice front due to the differing ice-front positions between the observation time of ERS-1 and the MOA coastline. Other small departures from approximately zero mean are seen, for example close to 550 km north, 1225 km west, where surface features such as large rifts exist and will have moved in the ∼10 year time period between the two instruments recording. Thus, a component of the RMS difference is due to the time interval between the two satellite datasets.
Figure 2b shows the difference between the ERS-1 gridded dataset and the ICESat data over the Brunt and Riiser-Larsen Ice Shelves. This is typical of the results on ice shelves with areas of no ERS-1 coverage. In the area around 1475 km north, 640 km west, differences are significantly larger than elsewhere and are both positive and negative. This area has patchy coverage from the ERS-1 data, so the interpolation misses real surface features. Again, larger differences exist close to the ice front and grounding line where ERS-1 data are absent. Figure 2c shows the differences over the Getz Ice Shelf. In this case, dH/dt measurements are larger than seen on the Brunt and Riiser-Larsen Ice Shelves, and including the effect of this change significantly improves the differences seen on the northern end of the ice shelf. Very large differences (mean differences in a 30 km × 30 km area of close to 50 m) are seen at a number of places close to grounded islands including close to 1030 km south, 1450 km west, where the ERS-1 data show a large gradient for a small number of data points close to an area of missing data.
Discussion of derived thickness estimates
Ice-shelf thickness is determined for all floating gridboxes using the method outlined above, and the resulting map of ice thicknesses is shown in Figure 3. A small number of gridboxes have firn corrections which are larger than the surface elevation observed, which produces a negative thickness. These gridboxes (1.3% of the total floating gridboxes) were set to zero ice thickness. They occurred mainly in the southern portion of Larsen D, where the firn model is unable to resolve the narrow ice shelf given its 55 km spatial resolution, and on the smaller ice shelves where heavy crevassing and ice melange is present and so the mean ice density is indeterminate. The ice-shelf thickness varies considerably from very thin ice <200 m thick, on some of the smaller ice shelves in Dronning Maud Land and on the Peninsula, to >2000 m close to some of the grounding line of the Filchner–Ronne and Amery Ice Shelves.
Ice thickness (m) for all ice shelves overlaid on the MOA image of Antarctica (Reference Scambos, Haran, Fahnestock, Painter and BohlanderScambos and others, 2007).
Comparison of ice-shelf thickness with airborne radar measurements
To determine the accuracy of the derived ice-shelf thicknesses, we compared them with independent estimates from airborne RES. The five datasets used here cover parts of the Ross Ice Shelf and the Pine Island and Thwaites Ice Shelves and are:
-
1. SOAR CASERTZ (US Support Office for Aerogeophysical Research–Corridor Aerogeophysics of the South East Ross Transect Zone) on the Ross Ice Shelf close to the grounding line, north of Siple Dome
-
2. WMB at the southeastern corner of the Ross Ice Shelf and along the Saunders Coast
-
3. RBG at the outlet of Robb Glacier onto the Ross Ice Shelf close to the grounding line on the western side
-
4. PIG on the Pine Island Glacier ice shelf
-
5. THW on the Thwaites Glacier floating tongue.
The data properties and acronyms are described in the following subsections. Data are compared by bilinearly interpolating the ice-shelf thickness grids to all airborne data points defined to be on floating ice. The mean difference (airborne-gridded satellite), expressing the systematic bias in the comparison, and the RMS difference, expressing the random error with no assumption of a Gaussian distribution, are determined for each dataset. Comparison with other publicly available BEDMAP ice-shelf thickness datasets shows similar results.
SOAR CASERTZ
The SOAR CASERTZ data were recorded by the University of Texas Institute for Geophysics and the US Geological Survey between 1991 and 1997 (Reference Blankenship, Alley and BindschadlerBlankenship and others, 2001). The available ice-shelf data cover the outlet of Bindschadler Ice Stream into the Ross Ice Shelf. The data were collected on a Twin Otter aircraft flying 300–1000 m above the ice surface carrying a pulsed 60 MHz radar capable of penetrating up to 4000 m of ice. Data were sampled every ∼12 m along track, and the instrument had a vertical accuracy of ±11 m for smooth topography such as ice shelves (Reference Blankenship, Alley and BindschadlerBlankenship and others, 2001). No firn correction was applied in the processing of these data, implying that the airborne data may be biased towards thinner ice thicknesses than reality.
Figure 4a shows the differences between the SOAR CASERTZ and the gridded ice-shelf thickness. The mean difference in the entire region is 27.3 ± 33.5 m, reducing to 21.9 ± 32.4 m when considering a zone within 10 km of the grounding line. While the differences are, in general, low, a distinct spatial pattern exists. In some areas close to the grounding line, where the assumption of hydrostatic equilibrium breaks down, differences are larger (up to and over 100 m) and negative, indicating that the airborne observations are lower than the satellite estimates. The mean difference within 2 km of the grounding line is −13.0 ± 52.1 m. In this region, due to the low surface slopes in the grounding zone, ERS-1 does not lose lock, so differences are due to the fact that the ice is not freely floating, rather than poor interpolation. A region of positive differences aligns with a fast-flowing outlet of the Bindschadler Ice Stream. A 100 km2 zone centred on the regions of largest difference shows a mean difference of 53.4 ± 54.8 m. Based on these positive differences, the ice is thicker than is implied from inversion of the satellite data; we suspect this may be due to the presence of a marine ice layer under the fast outlet region, and thus a higher column-averaged ice density than considered in the hydrostatic equilibrium calculation. Areas of basal freezing have been shown to exist in this area (Reference Joughin, Tulaczyk, MacAyeal and EngelhardtJoughin and others, 2004). Errors in this region are <80 m and within 10–15% of the total ice thickness. The overall positive bias across the dataset possibly indicates ice of a higher density than the assumed meteoric ice density or a thin layer of marine ice extending across the entire region. A 100 km2 region representative of this overall bias has a mean difference of 24.7 ± 24.9 m. Assuming no error in any parameter other than meteoric ice density, a bias of ∼25 m in this area indicates a meteoric ice density of 921 kg m −3, which is within our error bounds for meteoric ice density and within observed densities on the Amery Ice Shelf (Reference Fricker, Popov, Allison and YoungFricker and others, 2001). A similar analysis on the region of highest positive difference indicates a column-averaged density of 925 kg m−3, equivalent to a marine ice thickness of 24% of the ice column given a meteoric ice density of 921 kg m −3, or 38% given a meteoric ice density of 917 kg m−3.
Difference between airborne radar ice thickness measurements and inferred satellite thicknesses on the Ross Ice Shelf (m). (a) Comparison with SOAR CASERTZ data north of Siple Dome; (b) comparison with the western Marie Byrd Land survey at the southeastern edge of the Ross; and (c) comparison with Robb Glacier survey at the western grounding line. Locations marked by red (a), green (b) and blue (c) boxes on Ross Ice Shelf map inset in (c). The differences are overlaid on the MOA image of Antarctica (Reference Scambos, Haran, Fahnestock, Painter and BohlanderScambos and others, 2007), the thick black lines show the position of the MOA grounding and coast lines and the green square indicates the colour of zero difference.
WMB
The Western Marie Byrd (WMB) Land Survey (Reference Luyendyk, Wilson and SiddowayLuyendyk and others, 2003; B.P. Luyendyk and D.S. Wilson, http://nsidc.org/data/nsidc-0119.html) was flown in December 1998. The data cover the Sulzberger Ice Shelf on the Saunders Coast and the Shirase Coast region of the southeastern Ross Ice Shelf. Data were collected using the same instruments as SOAR CASERTZ flying at least 300 m above the surface, with the radar producing good images of the ice for thicknesses less than 1000 m but rarely imaging through ice thicker than 1500 m. The track spacing was either 5.3 or 10.6 km in a grid formation over the target region.
Figure 4b shows the difference between the airborne- and satellite-derived data. It is immediately obvious that on the smaller Sulzberger Ice Shelf, agreement is significantly worse than on the Ross Ice Shelf. The mean difference seen is 39.9 ± 144.6 m, rising to −71.8 ± 197.5 m within 10 km of the grounding line (excluding grounding lines with grounded islands). Away from the grounding line on the Ross Ice Shelf, the agreement between the datasets is again better than 90%. Closer to the grounding line, ERS-1 loses lock and, as shown in Figure 1a, areas of ∼10 km width have no data and are purely interpolation in the ice thickness product. This, combined with the breakdown of hydrostatic equilibrium, explains the larger differences (rising to >20%) seen here. The Sulzberger Ice Shelf is likely typical of many small ice shelves. The ice is pinned both at the grounding line and by a number of islands, meaning that relatively little of the ice is truly freely floating. Around 650 km west, there are a number of flight-lines with large differences between the airborne- and satellite-derived ice thicknesses. They are landward of the MOA grounding line plotted in Figure 4b. This is one of the areas where the MOA grounding line was edited to include areas thought to be floating from InSARderived grounding lines. The large differences seen suggest that MOA is probably a more accurate grounding location in this area and the ice was probably grounded.
RBG
The Robb Glacier (RBG) survey was flown by the Institute for Geophysics, University of Texas, in 1999–2000. The instrument set-up was as for the SOAR CASERTZ campaigns, and a single flight-line was flown down Robb Glacier, through the Transantarctic Mountains onto the Ross Ice Shelf on the western side. The mean difference between the datasets is −2.8 ± 48.6 m. In terms of relative difference, this is well below 10% for the third of the profile closest to the centre of the ice shelf. The datasets are also within 10% of each other close to the grounding line. In the rest of the profile, differences are both negative and positive and are of a similar width (∼1 km) to features seen on the surface in the same area in the background MOA image. This suggests that the differences seen are due to small-scale features such as crevasses, which are resolved by the airborne instruments but below the resolution of the ERS-1 RA.
PIG
Two profiles have been flown down the Pine Island ice shelf as described by Reference Bamber and RignotBamber and Rignot (2002). The data were collected in 1998 and use kinematic global positioning for navigation control and should be accurate to ±15 m (Reference Bamber and RignotBamber and Rignot, 2002). Figure 5a shows the comparison between the satellite-derived ice thicknesses and the airborne measurements. The mean difference is −48.2 ± 105.5 m and shows areas with very large differences. While 45% of Pine Island ice-shelf gridboxes contain ERS-1 data, the spatial distribution is patchy, with large areas of no data coverage. The major mismatch between the datasets is seen at the inland side of the ice shelf where there are no ERS-1 data for the first ∼35 km. This region exhibits a large discrepancy between the grounding line derived by MOA and those derived from InSAR measurements. The InSAR grounding line is some 30–35 km seaward of the MOA estimate, and if the comparison is repeated using just the former, the differences improve to −42.7 ± 94.5 m, so the region of largest difference close to the grounding line in the figure is probably part of the grounded ice plain. The remaining differences are most likely due to the complex surface topography on the ice shelf. As Table 2 shows, the errors in the surface elevation are more than large enough to account for the errors seen in the thickness measurements. This highlights the difficulty in determining ice thickness on small ice shelves with complex topography and a highly dynamic nature.
Difference between airborne radar ice thickness measurements and inferred satellite thicknesses (m) on (a) the Pine Island Glacier ice shelf and (b) Thwaites Glacier floating tongue (locations marked by red (a) and green (b) boxes on Amundsen Sea sector map). The differences are overlaid on the MOA image of Antarctica, the thick black lines show the position of the MOA grounding and coast lines (Reference Scambos, Haran, Fahnestock, Painter and BohlanderScambos and others, 2007) and the green square indicates the colour of zero difference.
THW
Two profiles were flown down the Thwaites ice tongue in the 1970s. The data have an accuracy of 1–10% primarily limited by the accuracy of the aircraft navigation, which utilized an inertial navigation system (Reference Bamber and RignotBamber and Rignot, 2002). On Thwaites, the MOA and InSAR grounding lines are in agreement to within ∼5 km, increasing confidence that all the area included is floating. The percentage of gridboxes containing ERS-1 data on Thwaites is 50% and is fairly uniform except for a ∼10 km width data gap near the grounding line and no coverage at the far eastern end of the floating tongue. Figure 5b shows the differences along the two profiles, which have a mean difference of 11.0 ± 70.5 m. The variability along track is significant and does not show any increase towards the grounding line, instead being high throughout the profile. This is probably due to the dynamic nature of the Thwaites ice tongue. The tongue is not a simple ice-shelf surface but is composed of a melange of sea ice, ice-sheet debris and wind-blown snow. Rifting in the tongue is significant, and large icebergs calve along a significant proportion of the grounding line (Reference RignotRignot, 2001). ERS-1 is unlikely to accurately resolve surface features of this type due to its low spatial resolution, which is reflected in the significant differences seen between the ERS-1 derived DEM and those from ICESat in Table 2. On such a dynamic ice tongue, calving significant proportions of the total ice area, the long time period between the airborne and satellite dataset makes meaningful comparison difficult as surface features in the airborne data will be calved by the 1995 timestamp of ERS-1.
Conclusions
Surface elevation measurements from the geodetic phases of ERS-1 recorded in 1994–95 have been used to calculate the ice thickness for all ice shelves with area over 10 km2 in Antarctica. The ERS-1 data were carefully edited to ensure that only data from floating ice were included in the gridding. The data were then gridded using kriging and combined with similar data from ICESat in regions south of the ERS-1 latitudinal limit. Comparison of the gridded ERS-1 and ICESat surface elevation in regions north of the ERS-1 latitudinal limit show high accuracy, particularly on the larger ice shelves. Differences were generally <5 m except in regions of poor ERS-1 data coverage close to the grounding line. The surface elevations were converted to ice thicknesses assuming hydrostatic equilibrium using a modelled firn correction that accounts for spatially variable temperature, accumulation and wind speed. Comparison of the derived ice-shelf thicknesses with independent airborne estimates of ice thickness shows good agreement on the large ice shelves. Differences of >20% are seen in some areas, but these can be explained by areas of missing data and the breakdown of the assumption of hydrostatic equilibrium close to the grounding line. Known errors can be combined and an error estimate on the order of 50 m is produced for Larsen C (Reference Griggs and BamberGriggs and Bamber, 2009a). This may be a significant underestimate in areas of unknown marine ice density and thickness, data gaps and close to the grounding line, which makes the production of an error estimate for the entire dataset problematic as shown by the larger differences seen in the validation in the Pine Island and Thwaites areas. The thickness grids have been shown to be of high quality and have already been shown to improve model predictions when compared to predictions using ice-shelf thickness estimates derived from sparse in situ measurements (Reference Khazendar, Rignot and LarourKhazendar and others, in press). Data from the recently launched CryoSat-2 could help with some of the problems discussed here, as there will be no loss of lock over break-in slope regions and the spatial resolution is significantly higher than that of ERS-1. However, other errors are due to the breakdown of the assumption of hydrostatic equilibrium and uncertainties in the integrated ice column density. Modelling efforts are also under way which will create higher-spatialresolution, time-evolving firn density corrections including the effect of melt which will reduce the uncertainty currently in this estimate. The dataset presented here will be updated in the coming years when these new data become available. The thickness grids can be obtained by contacting the authors (j.griggs@bristol.ac.uk).
Since the manuscript of this paper was submitted, new modelled firn density corrections have become available. These were created by the same model running at a higher spatial resolution than used previously. This has the effect of better resolving features close to the grounding line where there are significant slopes. The final product supplied by the authors will contain these new firn density corrections. The overall effect of this change is small (mean difference between thickness using the two corrections is 0.2 m), whereas regional differences can be significant (up to 100 m), showing the importance of accurate firn-density measurements. The conclusions of the comparison with airborne ice-shelf thicknesses are unchanged, but the values change slightly. The CASERTZ mean difference becomes 19.8 ± 32.8 m, the WMB mean difference becomes −48.4 ± 147.35 m and the RGB mean difference becomes 42.2 ± 63.3 m. The mean differences for PIG and THW become −33.5 ± 112.1 m and −13.2 ± 72.1 m respectively.
Acknowledgements
This work was funded by UK Natural Environment Research Council (NERC) grant NE/E004032/1. We acknowledge the ice2sea project, funded by the European Commission’s 7th Framework Programme through grant No. 226375, ice2sea manuscript No. 27. We thank H.J. Zwally for the elevation change data and M. van den Broeke for providing the firn correction. We also thank all the airborne data providers for allowing us access. The SOAR data are based on work supported by the US National Science Foundation (NSF) Office of Polar Programs (OPP) under grants OPP-9120464, 9319369 and 9319379. RGB data were from the University of Texas Institute for Geophysics collected under grant NSF/OPP award No. 9319379. WMB data were funded by NSF contract NSF OPP 9615281 and obtained from NSIDC. We thank H.A. Fricker and an anonymous reviewer for their helpful comments which substantially improved the manuscript.
