Introduction
Basal melting of Antarctica’s floating ice shelves accounts for 15–35% of the total mass loss from the ice sheet, with iceberg calving removing the overwhelming majority of the remainder (Reference Jacobs, Hellmer, Doake, Jenkins and FrolichJacobs and others, 1992, Reference Jacobs, Hellmer and Jenkins1996). Despite the generally frigid oceanic conditions close to the ice shelves, water temperatures are ~1°C above the freezing point at depth (Reference Jacobs, Hellmer, Doake, Jenkins and FrolichJacobs and others, 1992). These temperatures are maintained year-round and, owing to the large specific heat capacity of water, represent an ample reservoir of heat available for melting. The addition of meltwater to the sea water that circulates beneath the ice shelves cools and freshens it, and one product of this process, Ice Shelf Water, is an ingredient of Antarctic Bottom Water (Reference JacobsJacobs, 1986; Reference Foldvik and GammelsrFoldvik and Gammelsrød, 1988). This latter water mass is observed at depth throughout most of the world’s oceans, and its northward flow is a major component of the global meridional overturning circulation. Despite this pivotal role in ice-sheet–ocean interactions, there are only a handful of direct observations of melting at the base of ice shelves. Almost all published figures are of steady state melt rates; that is, the melt rate required to maintain the ice shelf in a state of equilibrium. The steady-state melt rate is deduced from measurements of the other terms in the mass-balance equation and an assumption that the thickness of the ice shelf at a given point in space is temporally constant. Such calculations have obvious limitations, such as the impossibility of determining the role of basal melting in driving ice-shelf thinning or retreat, or of diagnosing any change in basal melting that might be a response to ocean warming or cooling.
The grounding line separating a marine ice sheet from its surrounding ice shelves has long been a focus of interest among glaciologists. Early theories suggested that the location of a grounding line could be inherently unstable (Reference WeertmanWeertman, 1974; Reference Thomas and BentleyThomas and Bentley, 1978), making the marine ice sheet of West Antarctica vulnerable to complete collapse (Reference MercerMercer, 1978; Reference Thomas, Sanderson and RoseThomas and others, 1979). More recently, this conclusion has been disputed (Reference Hindmarsh and PeltierHindmarsh, 1993, Reference Hindmarsh1996; Reference Le Meur and HindmarshLe Meur and Hindmarsh, 2001), although the processes that determine the location of the grounding line are still subject to uncertainty (Reference Le Meur and HindmarshHindmarsh and Le Meur, 2001). One process that is rarely considered in analyses of ice flow in the vicinity of a grounding line is the change in the basal mass balance of the ice as it goes afloat. Melting beneath grounded ice is driven by geothermal heat flux and frictional dissipation and is generally of the order of centimetres per year or less. The thickest floating ice is found immediately downstream of the grounding line. Since the depression of the freezing point of sea water with depth leads to high melt rates beneath thick ice, melt rates near grounding lines are typically of the order of metres per year (Reference Jacobs, Hellmer, Doake, Jenkins and FrolichJacobs and others, 1992). This rapid change in the basal mass balance as ice flows across the grounding line has interesting implications. Thinning of grounded ice would allow sea water beneath it, the basal melt rate would rise, increasing the thinning rate, and the ice would be unlikely to reground. Conversely, if thickening of floating ice were to exclude the underlying sea water, the melt rate would drop, causing the ice to thicken further. In order to understand the processes that lead to a stable grounding-line position, we must understand not only the complexities of ice flow in the transition zone from ice sheet to ice shelf but also what controls the basal melt rate near the inland margins of an ice shelf. The work reported here represents a first step towards gaining an understanding of the latter.
Over the period from early January to late December 2001 we conducted a series of experiments to measure the basal melt rate near the grounding line of Rutford Ice Stream, in the southwestern corner of the Filchner–Ronne Ice Shelf (Fig. 1). The key to our technique is a precise measurement of the ice-shelf thinning rate by phasesensitive radar. The thinning rate can be partitioned between vertical strain and melting, given contemporaneous measurements of the vertical strain rate. The advantages of our technique are its high spatial and temporal resolution and the fact that we measure the melt rate relatively directly and need make no assumption about the state of equilibrium of the ice shelf. Reference Corr, Jenkins, Nicholls and DoakeCorr and others (2002) described the technique and its application to a site on the George VI Ice Shelf. Here we briefly review that description and provide in the Appendix a detailed derivation of the principal equation we use to describe the thinning of the ice shelf.
Fig. 1. (a) Map of the Weddell Sea sector of Antarctica, indicating the study area near the grounding line of Rutford Ice Stream. (b) Enlargement of the study area showing survey network (+), radar sites (×) and the location of a 43 day tidal record obtained by Reference StephensonStephenson (1984) (o). The black dotted line shows where GPS and seismic measurements were obtained by Reference VaughanVaughan (1995) and Reference SmithSmith (1996), respectively. Shading indicates surface elevation (m), contoured from airborne survey data obtained at locations marked by white dots. The bold black line indicates the grounding line as inferred by Reference Doake, Alley and BindschadlerDoake and others (2001) from satellite radar interferometry, and ice flow is from left to right. (c) Contours of ice thickness generated from the same airborne survey data. Vertical bars at the numbered radar sites indicate the derived basal melt rate.
Theory
Most glaciological observations made from the surface of a glacier or ice sheet are inherently Lagrangian in nature, being made at a series of moving points marked by stakes planted in the surface. We exploit this aspect of our phasesensitive radar measurements by considering the mass balance of the moving column of ice beneath our survey markers. We use the radar to determine how the thickness of the column between two radar horizons, typically a strong internal reflection, assumed to originate from a material surface within the ice, and the ice-shelf base, evolves in time. Since the internal reflector may lie within the firn in the upper part of the column, we cannot assume that the entire column is incompressible. We use a more general form of the mass-balance equation:
derived from the assumption that within the firn, vertical compaction and the viscous response to deviatoric stresses are independent processes, the effects of which are additive. In deriving the compaction rate, we have made the conventional assumption that the density of the firn is a function only of the time elapsed since its deposition at the surface and that the accumulation rate has been constant during the period over which the present firn column was deposited. The Appendix contains a more detailed presentation of these assumptions and we discuss their validity in a later section of the paper. In the above expression, u is a three-dimensional velocity vector with components (u, v, w) parallel to orthogonal (x, y horizontal and z vertical, positive upwards) axes, while 
is the surface mass flux, defined as negative for a downwards flux, and ρ is the firn/ice density.
Integrating the above equation between the upper radar reflector and the ice-shelf base, we obtain an expression for the conservation of mass in the moving column of ice:
where H
e is the effective ice thickness derived from conversion of the two-way travel time of the radar pulse transmitted through the firn/ice to a distance using a constant scale factor equal to half the radar velocity in solid ice. The unknown basal mass flux is denoted 
, while t denotes time and n refractive index, which we have assumed to be directly proportional to density. The subscript i indicates values for solid ice, and overbars indicate depth-averaging. The term ρ(hu
) refers to the density at the depth of the upper reflecting horizon. If this upper reference horizon is deep enough to lie within solid ice, the third term on the lefthand side is zero, and the change in effective thickness is simply the sum of the contributions from horizontal convergence/divergence of the ice flow and basal melting/freezing. Otherwise the size of the third term, which physically results from compaction beneath the reference depth, must be estimated. A complete derivation of this expression appears in the Appendix. Note that we have dropped two terms associated with the covariance of density and horizontal velocity over the depth of the ice column. These are zero either if the horizontal velocity is constant with depth, or if the upper reflector lies within solid ice.
We can apply this same equation to the thickness of ice and/or firn, Hei , measured between two internal reflectors. In this case, the last term on the lefthand side is zero, since both upper and lower interfaces are material surfaces, and the rate of change of reflector separation depends only on strain and compaction. If the internal reflectors are located within the solid-ice part of the column below the firn layer, the compaction term drops from the equation describing temporal changes in their separation. We then arrive at a simple expression for the horizontal divergence term that appears above:
Making use of the above two equations, we can derive the basal melt rate of an ice shelf directly from repeat radar sounding. If suitably deep internal reflectors are not available, the horizontal divergence can be measured separately by the repeat survey of a pattern of surface markers. Knowledge of the horizontal divergence also enables us to extract information on firn compaction from the relative motion of internal reflectors within the upper layers of the ice shelf.
Method
The practical application of the technique is illustrated in Figure 2. At the start of the experiment a column of ice was marked and the thickness between a prominent internal radar reflector and the ice-shelf base was measured. Each thickness measurement actually comprised six individual soundings, spaced at short intervals, to ensure the correct identification of stable internal reflectors. The radar was a step-frequency system consisting of a vector network analyzer (HP8751A) and a pair of identical broadband antennas, which were separated by a distance of 5 m (Reference Corr, Jenkins, Nicholls and DoakeCorr and others, 2002). We used 1601 frequency steps with an interval of 32 kHz, giving a response equivalent to a pulsed radar system having a centre frequency of 295.6 MHz and a bandwidth of 51.2 MHz. The broad bandwidth ensured good spatial resolution, and since both the amplitude and phase of the radar signal were recorded, changes in the range of reflectors could be measured to a small fraction (~1%) of the wavelength (56.8 cm), giving centimetre-level accuracy in the ice-thickness changes.
Fig. 2. Schematic picture of a marked ice column at times t 1 and t 2. The distance h u–h b is measured at both times using phase-sensitive radar. The difference is the result of horizontal divergence and basal melting. The former is determined from the relative motion of either survey poles planted in the ice-shelf surface or internal reflectors identified on the radar record, leaving the latter as the only unknown.
Survey poles were set in the snow surface around the site of the radar sounding (Fig. 1b), and inter-stake distances were measured using global positioning system (GPS) techniques. The stake pattern was triangular to provide measurements of extension or compression in three directions. Trimble 4000 series GPS receivers were used to record dual-frequency, carrier-phase data, sampled at 5 s intervals. FastStatic processing techniques then yielded inter-stake distances to ~1 cm or better with observation times as short as 8–10 min. For our experiments the absolute motion of the ice shelf was of secondary importance, so we adopted a slightly unusual survey technique that used no fixed reference stations. Instead we used multiple roving receivers and designed the survey such that all the required inter-stake distances were directly observed GPS baselines, ensuring no degradation of accuracy through the need to infer distances from the measured coordinates of the (moving) survey markers.
On subsequent visits all the measurements were repeated. The initial and repeat radar soundings provided the ice thickness and material derivative of thickness that appear in Equations (1) and (2). Usable internal radar reflectors appeared from near the surface to at least 600 m depth at all the sites studied, so horizontal divergence and hence basal melting could be derived from the radar data alone. In this paper, we also make use of the results of the GPS surveys, as these provide improved temporal resolution of the horizontal divergence term, albeit an areal average over a region surrounding the radar site. For calculating the melt rate, we favoured the value derived from the radar measurements over a full year, this being a point measurement on the ice column of interest taken over a timescale long enough to average out short-term variability.
The sites discussed here (Fig. 1) were established in the austral summer of 2000/01, with repeat visits being made that summer over periods ranging from a few days up to 2 weeks. These observations provided information on shortterm variability in the ice deformation, dominated by tidal forcing, and summertime melt rates. All sites were revisited in 2001/02, giving estimates of long-term strain, compaction and year-round melting.
Study Area
Rutford Ice Stream is an active outlet glacier that discharges 14–18 Gt a−1 of ice into the southwestern corner of the Filchner–Ronne Ice Shelf (Reference CorrCorr and others, 1996; Reference Rignot and JacobsRignot and Jacobs, 2002; Reference Joughin and PadmanJoughin and Padman, 2003). This is just under 10% of the total inflow from the grounded ice sheet to the ice shelf (Reference Joughin and PadmanJoughin and Padman, 2003). Reference Doake, Alley and BindschadlerDoake and others (2001) give an overview of the ice stream, its catchment basin and its grounding zone. The latter is complex (Fig. 1), with a sinuous grounding line and ice thickness ranging from >1600 to <2000 m. Ice velocities peak at around 400m a−1 near the grounding line, and decrease steadily downstream (Reference Stephenson and DoakeStephenson and Doake, 1982; Reference Jenkins and DoakeJenkins and Doake, 1991). In the centre of the ice stream a ridge of grounded ice extends 15–20 km downstream. The bed knoll that produces this feature carves out a channel of thinner ice, so that a surface trough appears downstream on the ice shelf. Either side are tongues of thicker ice that produce ridges on the ice shelf. The ridge on the western side contains the thickest floating ice, but smaller channels of thinner ice, which are presumably also formed by bed irregularities upstream of the grounding line, are apparent within it. Lateral thickness gradients are large and the ice is not locally in isostatic equilibrium with the underlying ocean so close to the grounding line. The surface troughs are more subdued than they would be if they were in equilibrium, with the channels supported above their equilibrium level by the surrounding ridges, which are consequently depressed below their equilibrium level. The radio-echo sounding profiles used to generate the maps in Figure 1 are not sufficiently close to define all the details of this complex structure. In particular, the downstream end of the grounding-line promontory has a steeper surface than is apparent in the contour maps, and the region immediately downstream is depressed below its equilibrium level (Reference StephensonStephenson, 1984). The process of adjustment towards equilibrium over these regions must be ongoing, and this has implications for some of the observations we discuss later. The grounding zone of Rutford Ice Stream has been the subject of much study since the initial survey of Reference Stephenson and DoakeStephenson and Doake (1982). Crucially for the work presented here, a 43 day tidal record was obtained by Reference StephensonStephenson (1984) using a tiltmeter sited close to the centre of our new network of survey stakes (Fig. 1b).
Results
GPS surveys
Since all the radar sites were situated close to the grounding line, in the zone of tidal flexure, we anticipated that there might be some signature of this tidal forcing in the horizontal strain field. We therefore adopted a survey strategy in which we observed baselines for a number of hours, then processed individual 10 min segments as static observations to produce short time series of inter-stake distances. Relative motion of the survey poles was small over the timescales of these individual baseline observations, and we found that our strategy of static processing accurately reproduced the baseline lengths we calculated from kinematic processing of the data (Fig. 3a). We predicted the tidal height (h T) at the time of each 10 min observation using tidal constituents calculated from analysis of Stephenson’s (1984) data. The records of baseline length were then fitted to a model of the form:
with the coefficients A 0–2 for each leg of the stake network being found by least squares. Figure 3 shows results for one leg over a variety of timescales. Clearly the tidal signal dominates if observations are made only over one summer, and even in the year-round measurements it could introduce significant error if it were not correctly removed. From the least-squares fits for each of the legs, we were able to calculate the long-term, viscous strain rate:
and the tidal, probably elastic (Reference VaughanVaughan, 1995), strain per metre of tidal elevation:
Fig. 3. Variations in inter-stake distance measured at the upstream end of the eastern stake network. (a) Continuous record obtained over 13 hours, showing results derived from kinematic processing (black line) and static processing of 10 min data segments (red dots). (b) Intermittent record obtained over 12 days (red dots) and prediction (black line) derived by fitting a linear trend plus a constant times tidal elevation to the measurements. The box near day 21 indicates the segment plotted in (a). (c) Record obtained over 1 year (red dots) and prediction (black line). The box near the start indicates the segment plotted in (b).
In both cases, principal components were calculated for each combination of three legs making up a triangle of the survey network. Results are plotted in Figure 4.
Fig. 4. (a) Long-term viscous strain rates and (b) short-term (probably elastic) strains per metre of tidal elevation, determined from relative displacement of survey poles (+).
At all points the steady, viscous strain rates (Fig. 4a) show longitudinal compression, approximately in the direction of flow, associated with the deceleration of the ice shelf. The thicker ice to the west spreads laterally, although at the upstream end of the stake lines it appears to be constrained by the grounding-line promontory to the east. The thinner ice, immediately downstream of the grounding-line promontory, experiences mainly lateral compression, presumably due to the spreading of the thicker ice either side. Note that only at the downstream end of the western line does the sum of the principal strain rates give a significant net divergence, although divergence is what we would normally anticipate on an ice shelf. Elsewhere the sum is either very small or significantly compressive. An obvious interpretation of this is that, if it were not for basal melting, the ice shelf would be thickening in the downstream direction (on the eastern line it is anyway, even with melting). However, we should emphasize here, in the light of a later discussion, that we have measured and plotted the horizontal strain rates at the ice-shelf surface.
The tidal strains indicate a broad zone of stretching perpendicular to the grounding line on the high tide. We also infer that there is a possibly narrower zone of compression that we have only sampled on the upstream triangle of the eastern stake line. At low tide the strain changes sign. This pattern of longitudinal extension and compression is consistent with bending of the ice shelf between a region of no vertical motion where the ice is grounded and a region of uniform rise and fall on the tide at some distance from the grounding line. A similar picture of surface strain can be inferred from the vertical deflections measured by Reference VaughanVaughan (1995) using repeat kinematic GPS surveys at low and high tide. Line AA′ of Vaughan’s (1995) study runs through our eastern stake network (Fig. 1b). Parallel to the grounding line we see compression everywhere on the high tide, for which we do not have an explanation. Presumably it is some effect of the complex two-dimensional shape of the grounding line. Once again, we emphasize that these are observations of strain at the surface. If the ice shelf is bending as a thin beam, we would anticipate that these horizontal strains decay linearly with depth, passing through zero at a neutral surface within the ice shelf, then growing to be approximately equal in magnitude but opposite in sign at the ice-shelf base. We investigate whether this is the case later in the paper.
Phase-sensitive radar: long-term observations
Derivation of basal melt rate from displacement of internal and basal reflectors
Figure 5a shows a typical radar record from site 1. The received signal fades away after the direct breakthrough from the transmitter, and we see progressively weaker internal reflectors, which eventually disappear into the noise at around 600–700 m depth. At 1569 m depth, a strong reflection from the ice-shelf base appears. Our subsequent analysis assumes that the internal reflectors are material surfaces, so that the motion of the reflectors precisely mirrors the motion of the ice itself. The processing of the radar data was described by Reference Corr, Jenkins, Nicholls and DoakeCorr and others (2002). From the processed record obtained on the initial visit, we select reflectors with relatively high amplitude and constant phase. This typically yields in the region of 60–100 suitable internal horizons and one or more good reflections from the ice-shelf base. The data obtained from repeat soundings are then scanned and the phase change at these discrete points recovered (Fig. 6). Where the reflectors are closely spaced, it is a simple matter to infer the integer number of wavelengths that must be added to the phase difference to obtain the overall change in range. However, once the reflectors vanish into the noise, we lose track of the phase, and the integer adjustment for the basal reflector is indeterminate. The best solution is to obtain a preliminary estimate of the change in range through cross-correlation of the amplitude of the basal reflections and use this to estimate the integer ambiguity. The wide bandwidth of our radar means that the rise time of a pulse reflected from a near-specular reflector is around 30–40 ns, equivalent to a few metres in ice (Fig. 7). The relatively sharp edges of the reflections enable the initial and repeat soundings to be aligned with a precision better than plus or minus one wavelength. Thus the cross-correlation removes the integer ambiguity, and the total change in range is simply the sum of the correlation delay and the phase difference. In a few cases, the shape of the return changed over the course of the year-long observations such that amplitude correlation was itself ambiguous. In these cases, we estimated the integer number of wavelengths of movement by extrapolating the changes observed over the shorter (1–2 week) time intervals, during which the shape of the reflector was sufficiently stable and the total movement was less than one wavelength anyway. The total year-round movement was once again obtained by adding in the observed phase difference. We discuss the extrapolation of short observations in more detail later.
Fig. 5. Radar record (upper sub-panels) and displacement (lower sub-panels) of internal (cyan dots) and basal (red dots) reflectors measured over periods of around 330 days at sites 1–14 (a–n). For site locations see Figure 1. In the lower sub-panels the black line indicates the vertical strain that would result from the measured horizontal divergence. The (six) blue lines represent linear fits to the displacement of internal layers between the two vertical dashed lines on the left, while the (six) green lines represent quadratic fits to the same data. The boxes in (a) indicate the areas shown enlarged for clarity above.
Fig. 5 (cont.). For caption, see facing page.
Fig. 6. Phase differences determined over a period of 333 days at prominent reflectors on one of the six radar records from site 1. (a) Raw data incorporating whole wavelength ambiguities. (b) Unwrapped data with the integer ambiguity determined where there are sufficient good-quality data points. The error estimates are derived from the signal-to-noise ratio of the individual reflectors, assuming a background noise level of –120 dB. For clarity the error bars have been omitted from other figures.
Fig. 7. Basal reflections obtained during initial (black) and repeat (red) soundings made 333 days apart at one of the six closely spaced points comprising site 1. (a) Raw data showing the increase in range over the observation period. (b) Processed data, with a time/range shift obtained by cross-correlation of the two records applied to the repeat sounding. Dashed red lines indicate the impact of adding (subtracting) one wavelength to (from) the calculated shift.
In the lower sub-panel of Figure 5a we plot as a function of depth the displacement of the prominent reflectors, measured as described above, that appear on repeat radar records obtained over a period of 333 days. The displacement is referenced to a single, strong reflection, in this case at 107m depth, and data from all six of the closely spaced soundings are shown. At all points below the reference level, the range increased over the observation period, implying that the ice column is thickening at this site. This is consistent with the observed horizontal convergence noted above. Nearer the surface the effect of horizontal convergence is offset by firn compaction, completely so for the near-surface layers where the range remains approximately constant. In principle, we could calculate the vertical strain between any two reflectors as ΔH ei/H ei0, but as the original reflector spacing becomes smaller so the errors in the reflector displacement are magnified. Instead we assume that below the firn layer, where compaction is negligible, the vertical strain rate is a simple function of depth and fit an appropriate model to the measured layer displacement by least squares. The strain at any depth can then be calculated from the gradient of this line.
At site 1, below 100 m depth, the displacement of the internal reflectors appears to be a linear function of depth (Fig. 5a), as we would anticipate if the vertical strain rate were constant with depth. Extrapolation of the best-fit, linear trend to the depth of the ice-shelf base indicates the displacement of the basal reflector that would result from strain of the ice column alone. The fit is a weighted least-squares estimate based on the data between the limits indicated in the figure and is done separately for each of the six repeat soundings. The observed increase in range of the reflector is less than the value obtained from the extrapolation, implying that some mass has been lost from the base. The difference (1.17 m over 333 days) gives us a melt rate of 1.28±0.05 m a−1 (Table 1). This value is an average of the six individual measurements, and the error reflects both the scatter in the estimates of vertical strain as well as the absolute phase error on the basal reflector. If we were to use the measured surface convergence to estimate the vertical strain, a slightly lower melt rate of 1.04±0.04 m a−1 would result. We would generally attribute such differences in the implied vertical strain rate to the existence of spatial gradients in the strain rate over the area sampled by the triangle of stakes surrounding the radar site, and we would therefore favour the result derived from the direct radar observation of strain rate on the column of ice under study, even though the formal error estimate is slightly higher.
Table 1. Estimates of basal melt rate on Rutford Ice Stream
Figure 5b–f show results from the remaining five sites on the eastern network. All show a similar pattern in the strength of the internal reflectors, such that they vanish into the noise at around 600–750 m depth. In all cases (except that of site 2), an estimate of the vertical strain rate based on a linear fit to the measured layer displacement as a function of depth is consistent with that inferred from the horizontal convergence within the errors of the respective measurements. This agreement implies either that horizontal gradients in the strain rate are negligible at these sites, or that they are nearly linear, such that the spatial averages sampled by the stake triangles are equal to the point values at the triangle centres. At site 2 there is a sharp change in the vertical strain rate at about 150 m depth, approximately coincident with a minimum in the amplitude of the internal reflections. We do not have an explanation for this phenomenon. Closure of crevasses would manifest itself in this way, in that not all of the horizontal compression would be translated into vertical extension. However, this would imply very deep-penetrating crevasses, and the measured surface compression is consistent with the near-surface, rather than the deep, vertical strain.
Despite the relatively unambiguous results for strain thickening of the ice at these sites, determination of the melt rates is complicated by the nature of the radar reflection from the ice-shelf base. Rather than a single, sharp reflection, we see a complex, often double-peaked return, so that determining the depth of the nadir point on the ice-shelf base is difficult. We would normally pick the leading edge of the reflection, but this is a less obvious choice in those cases where later parts of the return are often significantly stronger. Worse still, at sites 2 and 3 the shape of the return changes in time, such that after an interval of 1 year it is no longer possible to match the various parts of the return. In these cases, the basal movement measured over short intervals of 1–2 weeks has been extrapolated over the longer interval. This procedure is justified in a later discussion of the short-term records. At site 2, using the measured horizontal convergence to determine the vertical strain rate would give us a melt rate of –0.24±0.15 m a−1 (i.e. basal freezing), but the linear fit to internal layer movement below the break point described above yields melting of 2.08±0.17 m a−1. At sites 3–6 the basal returns have two distinct peaks and we have estimated the motion of both. Calculated melt rates range from –0.51 ± 0.18 m a−1 to 1.42 ± 0.05 m a−1 (Table 1). Because of the good agreement between the two estimates of vertical strain at these latter sites, it makes little difference which we use.
Where there are two identifiable peaks in the basal return, the apparent differential motion is up to ~1 m a−1, with the more distant peak always showing less motion, implying a higher melt rate. There are several possible explanations for this. It could simply reflect an additional uncertainty of around plus or minus two wavelengths in our measurement of range changes. However, in the case of site 4, where we see maximum differential motion, both peaks were identifiable throughout the observation period, so we used the more reliable amplitude correlation to estimate the integer number of wavelengths moved by each reflector over a full year. Another possibility, that the first reflectors could be internal rather than basal, is suggested by the fact that the derived melt rates are generally close to zero. Such a strong internal reflector so deep in the ice column would be surprising, but Reference SmithSmith (1996) identified seismic reflectors, possibly englacial debris, 100–250 m above the base of the ice shelf along a 4 km line running approximately parallel to ice flow within ~100 m of sites 1 and 2 (Fig. 1b). A final possibility is that we have simply measured differential melting of a rough ice-shelf base, with protrusions melting more rapidly than hollows.
Results for the eight sites on the western stake line are shown in Figure 5g–n and are summarized in Table 1. In most cases, the integer ambiguity was removed by the more accurate method of amplitude correlation of the basal reflections obtained 1 year apart, but at sites 9 and 12 we had to rely on extrapolation of the short-term measurements, as discussed earlier. We find generally larger discrepancies than on the eastern lines between the vertical strain rates inferred from the survey of the stake network and those determined from a linear fit to the displacement of the internal radar reflections. Once again we assume these discrepancies arise from spatial gradients in the strain rate, not surprising given the high lateral thickness gradients (Fig. 1c), and we therefore put more trust in the values derived from the radar measurements alone. Most of the derived melt rates lie between 0 and 1 m a−1, but there are extremes of below –1 m a−1 and nearly 2 m a−1.
One aspect of the above results that is puzzling is the apparent occurrence of basal freezing at sites 3, 7, 13 and 14. We would not anticipate basal freezing so close to areas of basal melting, unless there were significant basal topography, which we do not observe. Further, if salty ice were accumulating at the base of the ice shelf, we would not expect to see a basal reflection at this range. At sites 3, 7 and 13 the apparent freezing rates are close to zero, and, although the formal error estimates suggest that they are in most cases significantly negative, a slightly larger error, which could result from a miscalculation of the integer number of wavelengths over which the base has moved, would permit a zero or slightly positive melt rate. However, the result at site 14 cannot be dismissed in this manner and indicates a problem with the analysis we have discussed so far.
Depth variation of vertical strain rate
A solution to the freezing conundrum is indicated by the motion of the deeper internal layers. Along the western stake line, we see a more gradual decline in the amplitude of the internal reflections and we can determine their motion down to depths of 800–900 m. Since we use a weighted least-squares technique to estimate the straight-line fit to these data, the slope of the line is determined primarily by the stronger near-surface layers, where the phase errors are smaller. Lower in the ice column, the measured reflector motion often appears to drift away systematically from the straight-line fit (particularly at site 14), suggesting that the vertical strain rate may not be constant with depth. If it is not constant, our estimate of the basal motion due to strain, and hence our estimate of basal melting, is incorrect.
A possible reason for depth dependence in the vertical strain rate is apparent in Figure 1. The sites are close enough to the grounding line that the ice shelf is not floating in isostatic equilibrium everywhere. Over the more pronounced variations in ice thickness, presumably generated by bed irregularities beneath the ice stream, the surface is still in the process of adjustment towards flotation. Vertical shear stresses must therefore be present in the ice shelf, as illustrated in Figure 8, to balance the gravitational forces acting to push the thin ice down and the thick ice up. The net effect of the vertical shear stresses is to induce relative horizontal compression at the top of the ice column and relative extension at the bottom where the ice is supported above its equilibrium level and is thus sinking (Reference Casassa and WhillansCasassa and Whillans, 1994). The converse is true where the ice is depressed below its equilibrium level and is thus rising. If we assumed the ice to be uniform and the shear stress to be constant with depth, the horizontal divergence would be a linear function of depth. The relative displacement of two internal layers is simply the integral of the vertical strain between them. Thus, with the vertical strain varying linearly with depth, we would expect to see layer displacements that are a quadratic function of depth. Adding a quadratic term to the least-squares fits to the data in Figure 5 also gives us a means of investigating the significance of the apparent nonlinearity while introducing minimal complexity to our model of layer movement.
Fig. 8. Deformation of neighbouring columns of ice that are not freely floating in isostatic equilibrium with the underlying sea water. Buoyancy forces are balanced by horizontal gradients in the vertical shear stress, which give rise to differential horizontal compression and extension over the depth of the ice column. If the horizontal divergence varies linearly with depth, the effect is the same as when the ice shelf is bent about a neutral surface.
The quadratic fits to the data are shown in Figure 5 and the implied melt rates are given in Table 1, along with the estimated uncertainties based on the standard errors of the regression coefficients and the scatter of the six independent lines fitted through the individual repeat soundings at each site. Although the quadratic terms are generally small, introducing them can alter the derived melt rate by >1 m a−1, because of the great thickness of the ice. At the sites where the internal reflections remain strong enough to determine phase changes down to 800–900 m depth, the quadratic lines give visually better fits to the deeper data even beyond the depth range used for the curve fitting. However, given the importance of the quadratic terms for our results, we also investigated their significance using an F test. For any model that we fit to our data we can evaluate the sum of squares of the deviations explained by the model:
and the sum of squares of the residual, unexplained deviations:
where the yi
are our individual measurements, 
is their mean, 
are the corresponding estimates from the fitted model, and n is the number of data points used in the fit. Assuming any particular model is correct, the unexplained deviations provide us with our best estimate of the true variance in the data arising from observational errors:
where m is the number of model parameters. We test for the significance of the higher-order term in the quadratic model by comparing the increase in the value of S e resulting from the inclusion of the extra model parameter with the remaining variance about the new fit:
If the value of F is significantly greater than 1, we can be confident that the higher-order model does a better job at explaining the total variance of the data about their mean value. Since the numerator has only one degree of freedom, we are effectively comparing the increase in the explained variance with our estimate of the true variance. For a denominator with 60–100 degrees of freedom, as we have here, F of around 7 indicates 99% confidence that the extra variance explained by the more complex model is significantly larger than the residual variance due to observational errors and any remaining unmodelled variance. F values for all sites are given in Table 1.
The F test gives us a very high level of confidence that at 10 out of our 14 sites the apparent non-linearity in the layer displacements is a genuine feature of the data that is well modelled by a second-order polynomial. In the remaining four cases, the non-linearity is weak, such that the difference between the melt rates derived from the two models is barely significant anyway. We find that allowing for a linear variation in the vertical strain rate with depth eliminates the apparent freezing rates described above at all but one site, and here the uncertainty in our model parameters now permits zero freezing or a small melt rate. The quadratic fits themselves also appear physically reasonable, in that the curvature (Fig. 5) generally fits with what we would expect from our simple picture of the shear stresses (Fig. 8) and the expected sign of the isostatic anomalies at each site.
Vertical shear stresses in the ice shelf
As a final check on this initially unexpected result, we estimate the likely magnitude of the vertical shear stresses induced by deviations from isostatic equilibrium and their impact on the vertical strain rate. For simplicity we assume that the thickness variations that give rise to the shear stresses and the associated bending of the ice are purely twodimensional. If the deviations of the surface elevation from equilibrium, δh, are entirely supported by vertical shear stresses, we can express the horizontal stress gradient as:
where the z axis is vertical and the x axis is chosen to lie in the bending plane. The shear stress is related to the shear strain rate by the flow law:
The rate factor, A ~ 10−18 Pa−3 a−1 (Reference Jenkins and DoakeJenkins and Doake, 1991), the exponent n = 3, and the effective strain rate is ~10−3 a−1 (Fig. 4a). Deviations from the equilibrium elevation are typically up to ~10 m (Reference Doake, Frolich, Mantripp, Smith and VaughanDoake and others, 1987), and the length scale of the anomalies is ~1–10 km. Thus we would expect to see horizontal gradients in the vertical shear strain rate of up to 10−7−10−6 m−1 a−1.
Rewriting Equation (2), such that we retain the depth-varying horizontal velocities, we can express the layer movement over time interval Δt as:
the vertical gradient of which gives us our estimate of the horizontal convergence or divergence as a function of depth:
Taking the second derivative, and using the assumption that the bending that gives rise to velocity variations with depth is purely in the x direction, we arrive at:
We have described the layer displacement as a quadratic function of depth:
so the quadratic coefficient gives us an order-of-magnitude estimate of the horizontal gradient of the vertical shear strain rate:
The quadratic fits to our year-round (Δt ~ 1 year) observations on internal layer displacement (Fig. 5) give values of A 2 lying in the range 6 × 10−9 to 9 × 10−7 m−1, and are thus entirely consistent with our simple theory of their origin. Although it was a surprise to find a depth variation in the vertical strain rate on an ice shelf, in hindsight, we should perhaps have anticipated it at this particular location.
Final results
The third column of melt-rate values given in Table 1, including the quadratic term, thus represents our best estimate of the true rates (also shown graphically in Fig. 1). We obtain values that range from near 0 to 2.5 m a−1. Melting is generally higher on the eastern stake line, where the ice is thinner. Along both lines the melt rates generally decline with distance from the grounding line. The numbers are consistent with several earlier estimates of spatially averaged steady-state melt rates (Table 2), but they appear inconsistent with the more recent estimates. The remote-sensing techniques employed by Reference Rignot and JacobsRignot and Jacobs (2002) and Reference Joughin and PadmanJoughin and Padman (2003) afford more uniform sampling of the ice velocity field than allowed by the earlier ground-based measurements, so the later results might be considered more robust. However, if melt rates of ~10 m a−1 are required for steady state, the ice shelf must be thickening rapidly. We discuss the state of equilibrium of the ice shelf in a later section.
Table 2. Estimates of melt rate required for equilibrium on Rutford Ice Stream
Phase-sensitive radar: short-term observations
Temporal variability of basal reflector motion
At all sites, we made repeated measurements over the course of the initial 2 week visit in January 2001, then again over the course of the 1 week revisit in December 2001. The aim was to investigate whether we could detect any short-term variability in the melt rate. The nearest oceanographic data have been obtained at a site that is some 300 km distant from the Rutford Ice Stream grounding line, on the western coast of Korff Ice Rise (Fig. 1). Despite the fact that this site is over 500 km from the open water, the data show a strong seasonality (Reference Nicholls, Makinson and WeissNicholls and Makinson, 1998), and our current understanding of the circulation within the cavity beneath the Filchner–Ronne Ice Shelf is that the inflowing waters that carry this seasonal signal to Korff Ice Rise are the primary source of heat for melting at the Rutford grounding line (Reference Nicholls, Makinson and JohnsonNicholls and others, 1997; Reference Jenkins and HollandJenkins and Holland, 2002a, b). We might therefore expect that year-long measurements of melting would differ from those obtained over a short timescale, unless that time was coincidentally one with melting near the average of the seasonal extremes. At Korff Ice Rise a rapid warming occurs in September–October each year, with relatively high temperatures persisting until early winter (Reference Jenkins, Holland, Nicholls, Schröder and ØsterhusJenkins and others, 2004), and we would not expect this signature to be more than a couple of months delayed in its progress to Rutford Ice Stream. Our January measurements should therefore pick up the maximum of any seasonal cycle. Established theory suggests that tidal mixing is likely to play a role in regulating heat transfer to an ice-shelf base in the vicinity of a grounding line (Reference MacAyealMacAyeal, 1984). We therefore also anticipated the possibility of melt rates changing over the fortnightly spring/neap tidal cycle.
The most intensively studied of our sites was site 1, where we made a series of ten measurements over a 15 day period in January 2001. The depth of the basal reflector is plotted in Figure 9a as a function of time. The thickening trend of about 13 cm over the 15 days is clearly visible and a simple linear fit to the data yields a rate of 2.95 ± 0.02 m a−1. This is lower than the 3.09±0.02 m a−1 apparent from the year-long measurement, which would suggest a slightly higher-than-average melt rate during the summer. Although the six measurements made at each time appear to be well clustered, the scatter of the groups about the straight line indicates significant deviations from the linear trend. A clue as to the cause of these deviations is given by the data from the eighth and ninth observations, which were made close to consecutive low and high tides on day 29. The implication is that the ice thickness undergoes a tidal oscillation, presumably as a result of net horizontal compression/extension induced by tidal bending of the ice shelf. We were also interested to see if these data can tell us anything about changes in melt rate over the spring/neap tidal cycle. Our hypothesis was that the melt rate should scale with the square of the tidal amplitude. We therefore fitted a model of the form:
to our thickness data and used an F test to determine the significance of the last two terms. We found that introducing the third term (F of 208) accounted for most of the observed variance about the straight line, while adding the fourth term (F of 1) provided no demonstrable improvement in the fit. The fit obtained using the first three terms only is shown in Figure 9a along with the underlying linear trend. Our best estimate of A1 gives us a thickening rate of 3.05 ± 0.02 m a−1, which is indistinguishable from the year-long trend, given the magnitude of the errors in our measurements. We therefore conclude that, at the level of certainty we can measure, the melt rate is steady, with no significant seasonal or fortnightly variation.
Fig. 9. Variations in ice thickness over a 2 week period at sites 1 (a) and 9 (c). Black crosses indicate measurements, while red lines indicate predictions derived by fitting a linear trend plus a constant times tidal elevation (b) to the measurements. The blue lines in (a) and (c) indicate the underlying linear trends in thickness.
This finding enables us to use the short-term observations of thickness change to determine the year-round changes by simple extrapolation, as discussed in the previous section. The advantage of the short-term measurements is that we can unambiguously track the basal motion. The disadvantage is that on these timescales that motion is subject to considerable tidal contamination.
At site 9 we made five measurements over a period of 10 days in January 2001 (Fig. 9c), so we can undertake a similar analysis to that described above for site 1. The overall thinning we observed is only about 3.5 cm, and a simple linear fit to the thickness data yields a trend of –1.23 ± 0.05 m a−1. Once again there is scatter about the straight line, and the inclusion of a term that is proportional to tidal elevation in our model improves the fit dramatically (F of 15, compared with a value of around 8 for 99% confidence with fewer degrees of freedom in this case). The underlying linear trend in the more complex model is –1.46 ± 0.10 m a−1. This knowledge makes it possible to identify the weak return in the December 2001 data and apply an appropriate integer wavelength correction to each of the six shots. Then adding the phase difference measured between the initial and repeat soundings yields an overall thickness change of –1.43 ±0.04 m a−1 (Fig. 5i), and this is the value we use in deriving the melt rate.
At sites 2, 3 and 12 we have only the minimal set of two measurements of thickness in January 2001 and two measurements in December 2001. However, with two independent determinations of thickness change it is in principle possible to separate the component caused by steady thickening/thinning from that due to the difference in tidal elevation. We consider this procedure to be less precise and we have incorporated an additional half-wavelength into the error estimates. The same technique has been applied to the secondary basal reflections at sites 5 and 6.
Ice-shelf deformation at tidal frequencies
The series of observations made at sites 1 and 9 provide, in addition to the results discussed above, a wealth of information on the response of the ice shelf to tidal bending. We have already shown the effect of net compression at site 1 and net extension at site 9 on the high tide (Fig. 9). In the case of site 1, the net horizontal strain is of the same sign as that measured at the surface (Fig. 4), while at site 9 it is of the opposite sign. This could suggest that the neutral surface is found at different levels in the two ice columns, but it is more likely to be a result of uncertainty in the sign of the much smaller surface divergence at site 9. We have seen that the response to a steady bending force is a vertical strain rate that is a linear function of depth (Fig. 8), so we might anticipate seeing a signature of this in the internal layer motion over the tidal cycle. Figure 10 shows layer motion as a function of depth between visits 8 and 9 and visits 9 and 10 at site 1. We have removed the steady motion of the internal layers as a result of the long-term strain, and that of the basal reflector as a result of both strain and melting. The residual displacements of the internal reflectors show the expected quadratic dependence on depth. Although noise in the internal layer displacements means that the quadratic terms are not well constrained by these data alone, the higherorder terms are nevertheless highly significant, with F values of 33 and 21 respectively. However, assuming that the motion we see in Figure 10 is due entirely to tidal bending, we can fit the basal reflector motion to the same model, and this procedure yields the lines shown in the figure.
Fig. 10. Displacement of internal (cyan dots) and basal (red dots) reflectors measured over opposite halves of the tidal cycle at site 1. The (six) blue lines indicate quadratic fits to the individual records. Black lines indicate the mean of the individual fits (solid) normalized by tidal elevation difference (dashed).
When normalized by tidal height difference, the two sets of observations give a consistent picture of the response to bending (Fig. 10). The neutral surface is at about 1000 m depth. Above that level, the ice column extends by about 10 mm per metre of tidal elevation, while beneath that level it contracts by about 4 mm. The tidal range at the Rutford Ice Stream grounding line is nearly 6 m (Fig. 9b), so that overall the tidally induced thickness change can be as much as 35 mm. Note that the full amplitude of the thickness oscillations is not apparent in Figure 9, because the depth of the basal reflector is plotted relative to the reference internal at 107 m, whereas about 25% of the extension above the neutral surface occurs in the upper 100 m of the ice column (Fig. 10). The vertical strain can be obtained from the gradient of the curves in Figure 10, and is plotted, normalized by tidal elevation difference, in Figure 11a. Included in Figure 11a are results from all nine time intervals between the ten visits to site 1. The uncertainty in these results rises as the tidal elevation difference falls, but there is a high degree of consistency between the observations. The only real outlier is that calculated over the time interval between visits 4 and 5. Visit 5 was made at mid-tide when the elevation was changing rapidly, so the most likely cause of the discrepancy is timing errors. Since our primary aim was the measurement of weekly to annual melt rates, we simply did not anticipate the need to time the observations with sufficient accuracy to place them precisely within the tidal cycle.
Fig. 11. Vertical strain, normalized by tidal elevation difference, as a function of depth at sites 1 (a) and 9 (b). Colour coding indicates separate observations with tidal elevation differences as indicated. Solid lines show the best fit, and dashed lines the 95% confidence intervals.
A weighted mean value (with the outlier removed) of the vertical strain at the ice-shelf surface is (2.04 ± 0.01) × 10−5 m−1 of tidal elevation. This compares remarkably well with a horizontal divergence of (–2.14 ± 0.15) × 10−5 m−1 determined from the survey network (Fig. 4). An ice shelf’s response to tidal loading is normally considered to be elastic, and the generally accepted value of Poisson’s ratio is 0.3 (Reference VaughanVaughan, 1995). Using this value and the measured vertical strain, we would infer a horizontal divergence of (–4.77 ± 0.03) × 10−5 m−1. Our measurements show that, at least at the ice-shelf surface and when the loading is applied at tidal frequencies, Poisson’s ratio is close to 0.5. From its surface value the vertical strain falls linearly to a value of zero (the neutral surface) at a depth of 960 ± 40 m, or 61 ± 3% of the way down the ice column. It is more than likely that the horizontal strain, which is the direct product of the bending, falls linearly with depth, implying that Poisson’s ratio must be 0.5 everywhere. From the gradient of the strain with depth we can determine the radius of curvature of the bend imposed on the ice shelf by the tide. On the extremes of the tidal excursion the radius reaches a minimum of about 15 × 106 m, four orders of magnitude greater than the ice thickness, so the conventional assumption that the ice shelf can be treated as a thin beam is clearly valid.
A similar analysis of the data from site 9 yields the vertical strains shown in Figure 11b. Once again there is one outlier, which is the result of a clearly noisy radar record. With this outlier removed, the mean strain at the surface of the ice shelf is (–5.8 ± 0.4) × 10−6 m−1 of tidal elevation. The measured horizontal divergence is (-1 ± 5) × 10−6 m−1. Given that we found discrepancies along the western line in our measurements of the long-term strain rate between spatially averaged values derived from the survey network and point measurements derived from the radar records, we consider these results to be consistent with the Poisson’s ratio of 0.5 determined at site 1. In this case, the neutral surface lies at 1140± 70 m depth, slightly deeper than at site 1, but at a remarkably similar relative depth of 59 ± 4% of the total ice thickness, while the minimum bending radius is about 65 × 106 m.
In principle, we could undertake a similar analysis on any of the short-term records. However, individual observations can clearly suffer from excessive noise if the difference in tidal height is small and/or the measurements were made at a time of rapidly changing tidal elevation. A dedicated pair of measurements timed to coincide with high and low tides would be sufficient, provided the steady long-term reflector motion were known.
Phase-sensitive radar: near-surface observations
Thus far we have considered the behaviour of the ice column at depth, where the density can be considered constant and the compaction term in Equation (1) is zero. To complete our discussion of the radar data, we now return to the long-term observations (Fig. 5) to investigate what the measured displacement of the near-surface layers can tell us about the compaction process. In Figure 12a the displacement is plotted relative to that at 150 m for all sites. The different trends reflect the differing horizontal convergences and divergences at each site, and this is the dominant signature through much of the firn layer. Clearly the firn layer does expand and contract in thickness in response to the horizontal compression and extension of the ice shelf. Our assumption in deriving Equation (1) (see Appendix) was that this process was volume-conserving, with firn densification proceeding independently as a function only of the rate of burial by new snowfall. Under this assumption, the residual layer displacement, after removal of the effect of horizontal divergence, gives us the rate of firn compaction (third term in Equation (1)), which should be similar at all sites. This is indeed the case (Fig. 12b), with the striking exception of site 2 where we have already noted that the vertical strain in the upper 150 m does not fit well with that inferred lower in the ice column.
Fig. 12. Reflector displacement measured in the upper 175 m at all 14 radar sites. (a) Total displacement referenced to 150 m depth. (b) Residual displacement after removal of the vertical strain determined from the relative motion of deeper reflectors (+). Orange circles indicate the residual after removal of vertical strain derived from measured surface convergence at site 2. The solid black line indicates the theoretical displacement derived from Equation (1) using the surface accumulation and density data of Reference Jenkins and DoakeJenkins and Doake (1991), with 95% confidence intervals given by the dashed lines.
We can estimate the theoretical size of the compaction term from knowledge of the surface accumulation rate and density profile. Reference Jenkins and DoakeJenkins and Doake (1991) measured an accumulation rate of 450 ± 65 kg m−2 a−1 at the location of site 1 and estimated a density–depth profile based on measurements of firn density in the upper 2–5 m at seven sites on the Ronne Ice Shelf. The compaction rate inferred from these data is also shown in Figure 12b. The agreement between these entirely independent determinations of compaction is remarkably good. Clearly, had we been unable to identify a reference layer deep enough to lie within solid ice, using a shallower reflector and an estimate of compaction beneath it based on the earlier observations of Reference Jenkins and DoakeJenkins and Doake (1991), i.e. the solid black line in Figure 11b, would have introduced additional errors of no more than a few centimetres per year to our estimates of basal melting. From this we conclude that the density profile appears to be in near equilibrium with a relatively steady surface accumulation rate.
The good agreement between observation and theory at sites where the net horizontal divergence ranges from –2.8 × 10−3 to 1.1 × 10−3a−1 supports our assumption that horizontal strain has no impact on the firn density, at least at the depths of our measurements. This result can be explained by the fact that firn compaction is driven by the isotropic stress:
which is dominated by the hydrostatic pressure. We would only expect densification to be influenced by other stresses if they were large enough to significantly alter the isotropic stress. Within the solid ice, longitudinal deviatoric stresses 
are typically ~105 Pa or less, but they are smaller in the firn layer because the effective shear viscosity, 
, is a strong function of density (Reference Ambach, Huber, Eisner and SchneiderAmbach and others, 1995). In the near-surface layers, where the density is around 500 kg m−3 or less (Reference Jenkins and DoakeJenkins and Doake, 1991), the effective viscosity, and hence the deviatoric stress, is likely to be an order of magnitude smaller than in the solid ice. Thus at a few metres depth the hydrostatic pressure will already have exceeded the longitudinal deviatoric stresses, and beyond that the latter will make an ever smaller contribution to the isotropic stress. At the depth of our measurements it is therefore safe to assume that horizontal strain is volume-conserving, while firn compaction is driven by the hydrostatic pressure.
Of course the hydrostatic pressure is itself modified by the vertical strain of the firn column. If the maximum strain rates we have observed were sustained over the entire period of the snow–firn–ice transition (around 150years), the pressure exerted on firn of any particular density could have been altered by up to 25% as a result of the cumulative strain. Assuming that the compaction rate is a cubic function of the isotropic stress (Reference Ambach, Huber, Eisner and SchneiderAmbach and others, 1995), this extreme case would lead to almost a factor 2 change in the compaction rate. However, these numbers are appropriate to ~100 m depth where the compaction rate is near zero anyway. At 50 m depth the maximum observed strain rate (acting over about 75 years in this case) would give a 10% change in the hydrostatic pressure and hence a 30% change in the compaction rate, while at 20 m (firn age around 25 years) the limits are 3% and <10% respectively. Thus even these unrealistic upper limits on the changes in compaction rate caused by the net strain of the firn column lie well within the errors of our observations (Fig. 12).
The curious behaviour at site 2 defies simple explanation. If we assume that the horizontal convergence is constant in the upper 125 m and equal to that measured at the surface, the residual layer motion is consistent with that at the other sites. This would appear to discount the closure of deep crevasses as an explanation of the anomalous behaviour in the upper 150 m, since we see no volume change other than that due to vertical compaction. However, this leaves us with an observation of a rapid change in the horizontal convergence at mid-depth in an apparently solid ice column. One possibility is that we are seeing the impact of crevasses that lie outside the stake network, but it is difficult to see why such a far-field influence should not similarly affect other sites.
We note that within the firn column, the last term in Equation (A22) of the Appendix, involving the covariance of density and horizontal velocity, is formally non-zero, because of the depth variation in the velocity we discussed earlier. However, we found that the vertical gradient of the horizontal divergence was at most ~10−6 m−1 a−1, so variations over the 100 m firn column are at least an order of magnitude, and are typically two orders of magnitude, lower than the mean divergence. Thus the covariance terms, while non-zero, are generally negligible.
State of Equilibrium of the Ice Shelf
Our analysis thus far has exploited the inherently Lagrangian nature of our observations, and we have therefore avoided the need to make any assumptions about the state of equilibrium of the ice shelf. Taking a more conventional, Eulerian view of our observations, we can investigate to what extent an assumption of steady state would be valid. The survey lines were set up approximately parallel to ice flow, so choosing a local x axis along this line allows a simple decomposition of the measured total derivative of ice thickness:
If the ice shelf is in a steady state, the first term on the righthand side is zero, and the steady-state thickness profile along a flowline can be calculated from
and compared with the measured thicknesses.
The ice velocity is available as a by-product of our GPS survey, even though we have not needed it in our previous analyses. We used no fixed reference station for the survey, but stand-alone point positioning, despite contamination from the vertical motion of the GPS antennas on the tide (Reference KingKing, 2004), is accurate to a few tens of centimetres. Our year-long observations of the displacement of the survey markers therefore give us ice velocity to a few per cent, and we have interpolated these data onto the lines of our radar sites (Fig. 13a). The derived steady-state thicknesses are compared with the observed thicknesses in Figure 13b.
Fig. 13. (a) Velocity (interpolated from measurements made at the survey poles shown in Fig. 1) as a function of distance along the four lines of radar sites. (b) Measured (+) and calculated (lines) thickness as a function of distance along the lines of the radar sites. The calculated thicknesses were derived from interpolations of velocity (a) and thinning/thickening rate (Fig. 5) measurements and an assumption of steady state.
On the eastern lines we see very good agreement between observation and calculation, implying that the ice-shelf thickness distribution is in equilibrium even on these relatively short spatial scales. On the western lines the measured thicknesses tend to be smaller than the steady-state values. Discrepancies are as large as 150 m, but this may not signify a departure from steady state. There are pronounced thickness gradients perpendicular to flow, clearly shown up by the difference in thickness between the two western lines, which are <1 km apart. It appears that the lines may have drifted slightly east of the true flowlines, such that the further west line has overlapped the nearer one, and the nearer one has drifted towards the narrow channel of thinner ice to its east (Fig. 1). Recall that the site furthest downstream that appears to be too thin in Figure 13b was the one most affected by vertical shear stresses, suggesting close proximity to a thickness anomaly.
We therefore find no strong evidence for a significant departure from steady state, a finding that is consistent with the observed steadiness of the grounding-line position (Reference RignotRignot, 1998). Derived steady-state melt rates should therefore be consistent with our directly measured melt rates. This places the numbers in Table 2 in a different light. While the earlier estimates, based on ground-based and airborne measurements, are entirely consistent with our new data, the more recent estimates, based primarily on satellite remote sensing, appear to be anomalously high.
Summary and Conclusions
We have described the application of an experimental technique, whereby the melt rate at the base of an ice shelf can be measured at high spatial and temporal resolution, to the region immediately downstream of the grounding line on Rutford Ice Stream. The study area is around 800 km from the open ocean and contains some of the thickest floating ice to be found anywhere on the planet. We find melt rates ranging from –0.11 ±0.31 to 2.51 ±0.10 m a−1, with an overall mean of 0.85 m a−1 and a standard deviation of 0.69 m a−1. There is no obvious global correlation with ice thickness or with distance from the grounding line, although there is a general downstream decrease in melting along each of our survey lines. We found no significant variation of the melt rate with season or with tidal range, suggesting that melting might be driven by relatively steady upwelling from a pool of dense water that is topped up, rather than completely replenished, by the seasonal inflow of new High Salinity Shelf Water from the ice front. Seismic soundings (Reference Smith and DoakeSmith and Doake, 1994) reveal a deep trough to the south of the survey area that presumably acts as a channel for any dense inflow. Within this trough the water column is about 500 m thick, so the vertical heat flux there is likely to be small. Heat is presumably delivered to the ice–ocean boundary layer closer to the grounding line, where the water column pinches out, and is gradually used up as the water within the boundary layer flows away from the source of upwelling. Although tidal forcing undoubtedly provides energy for vertical mixing, we see no evidence at this location that it influences the rate at which heat is delivered to the ice–ocean boundary layer. Perhaps the input of kinetic energy is always sufficient to keep the water column well mixed. Another possibility is that tidal mixing maintains a variable mixed-layer thickness, but there is steady divergence of the mixed-layer flow such that a near-constant entrainment rate is required to maintain that thickness.
Our observations of melting compare well with several earlier estimates of the melt rate required for steady state, and this close agreement is consistent with our observation that the ice-shelf thickness distribution is presently close to a steady state. However, our numbers are significantly lower than some more recent estimates of melting that are also based on the assumption of steady state. We consider it unlikely that the discrepancy is simply a result of comparing point measurements (the radius of the first Fresnel zone of our radar at the range of the ice-shelf base is 21–24 m) with spatial averages. Our measurements are fairly well spread, sample a wide range of ice thicknesses, and cover distances from the grounding line that vary from <3 to around 20 km, measured along flow. All our results are of a consistent magnitude and show no signs of under-sampled spatial variability, which we would not expect on theoretical grounds. The earlier estimates, with which our results compare well, are also averaged over a wide range of spatial scales (Table 2). We therefore conclude that the discrepancy results from errors in the more recent estimates.
The advantages of using satellite-derived data in the estimation of melt rates are the continent-wide coverage and uniform sampling. The major disadvantage is that ice thickness cannot, as yet, be directly measured from space. On ice shelves the thickness is commonly derived from surface elevation, measured by satellite altimetry, and an assumption that the ice is freely floating. However, altimeters do not perform well at the ice-shelf margins, where the slope of the ice surface introduces an error (Reference Brenner, Bindschadler, Thomas and ZwallyBrenner and others, 1983) and, as we have already discussed, the assumption of free flotation is not uniformly valid. Figure 14 shows how the observed thickness (Fig. 1) along the Rutford Ice Stream grounding line compares with that derived from three pan-Antarctic digital elevation models, based primarily on satellite altimetry (Reference Bamber and BindschadlerBamber and Bindschadler, 1997; Reference Liu, Jezek, Li and ZhaoLiu and others, 1999, Reference Ambach, Huber, Eisner and Schneider2001). There are considerable errors, with mean deviations of up to 25%, in the derived ice thickness, and we note that even if the measured elevations were perfect, the thickness would be overestimated by nearly 10% because the ice at the grounding line is far from freely floating. Downstream of the grounding line, the ice is closer to flotation and its surface is flatter, so errors in the derived thickness will be smaller. Melt rate estimates are based on the change in ice flux over a given distance, so such non-systematic errors in the thickness data give large errors in the calculated melt rate, and could be a generic problem in the vicinity of a grounding line. Another possible source of error is the tidal modulation of ice-shelf motion that could introduce bias into measurements of displacement made over a short period of time.
Fig. 14. Surface and basal elevation of Rutford Ice Stream along the grounding line indicated in Figure 1. The distance axis runs from east to west (top to bottom of Fig. 1). Note the break in scale at sea level. Dark shading indicates the surface and base observed by airborne radar sounding. Lighter shading indicates the base that would be inferred from the observed surface elevation and an assumption that the ice is freely floating. Surface elevation extracted from three commonly used Antarctic digital elevation models (Reference Bamber and BindschadlerBamber and Bindschadler, 1997; Reference Liu, Jezek, Li and ZhaoLiu and others, 1999, Reference Liu, Jezek, Li and Zhao2001) is also shown along with the inferred basal elevation associated with each. ERS-1: European Remote-sensing Satellite 1; RAMP: RADARSAT Antarctic Mapping Project.
Almost as a by-product of our melt rate measurements, we have gained unprecedented insight into the response of the ice shelf to tidal forcing. As far as we can tell, the strain of the ice is directly proportional to the tidal elevation, and thus presumably the applied stress, suggesting an elastic response. The bending radius is very large, implying that the application of thin-beam theory to the problem is valid, but we also see evidence for the complex, two-dimensional nature of the bending. Our measurements have allowed us to evaluate a Poisson’s ratio of 0.5, to identify a neutral surface at a depth equivalent to 60% of the total ice-shelf thickness and to observe the resulting tidally induced thickness oscillations. We have also made high-resolution measurements of the impacts of horizontal strain and compaction on the firn column and shown that the two are essentially independent, at least below 10–20 m depth (firn density above 650–700 kg m−3) and at the level of accuracy of our observations (a few cm a−1 of layer motion). Finally, the precision and detail of our radar measurements have enabled us to observe directly the impact of vertical shear stresses that are present within an ice shelf that is not freely floating. We have found variations of the vertical strain rate of up to ~10−3a−1 over the depth of the ice shelf, enough to have a significant impact on our calculations of the melt rate.
The experimental technique we have discussed is the first, to our knowledge, that provides a direct measure of the melt rate at the ice-shelf base without the need to penetrate the ice shelf. As a result, the method offers unprecedented spatial and temporal resolution. Compared with traditional techniques used to estimate steady-state melt rates, such observations can provide entirely new insight into the processes operating in the ocean beneath Antarctica’s ice shelves.
Epilogue
In January 2001 we set up what we anticipated would be a simple experiment to measure the basal melt rate of the Ronne Ice Shelf near the grounding line of Rutford Ice Stream. Our pilot study on the George VI Ice Shelf (Reference Corr, Jenkins, Nicholls and DoakeCorr and others, 2002) had yielded a model dataset, from which we could easily extract the basal melt rate. However, on Rutford Ice Stream the melt rates proved to be much more elusive. To get at them we had to measure not just the oscillations in vertical strain rate and ice thickness caused by the rise and fall of the ice shelf on the tide, but also the depth variation of the long-term strain rate caused by deviations of the ice-shelf surface from the level of isostatic equilibrium. When we were designing the experiment we had no idea that these effects would be of an order of magnitude that would bother us, or that we would have to measure them. The fact that we have been able to measure them is a testament to the immense potential of the phase-sensitive radar as a tool in glaciology.
Although we have described Lagrangian measurements of thickness change, modern GPS techniques allow precise reoccupation of a geographically fixed point and make Eulerian measurements of thickness change equally feasible. Hence, the original application of phase-sensitive radar envisaged by Reference Nye, Berry and WalfordNye and others (1972), Reference NyeNye (1975) and Reference Walford, Holdorf and OakbergWalford and others (1977) to measure the thickness change of an ice sheet is readily achievable with off-the-shelf instrumentation. Although the geographical coverage attainable with such a technique will never rival that provided by satellite altimeters, as a complementary measurement it would be immensely powerful. Not only is it a direct measure of the thickness change, rather than the change in surface elevation, but it is also possible to measure independently the long-term thickness evolution below the firn layer as well as the shorter-term effects of firn compaction. In addition to thickness changes, phasesensitive radar observations provide a wealth of detail on the internal deformation of an ice mass. Our identification of depth dependence in the strain rate of an ice shelf, in particular the curious behaviour observed at site 2, provides a hint of the new insight that could result if such observations were routinely made on glaciers and ice sheets.
Acknowledgements
We gratefully acknowledge the invaluable assistance provided by C. Day and D. Routledge in collecting the data presented above. A lively and informative discussion with S. Evans and J. Nye helped us to clarify the assumptions we made in deriving Equation (1).
Appendix: True and Observed Motion of Radar Horizons in an Ice Body
The observations we make with our phase-sensitive radar system yield a precise record of the temporal evolution of the two-way travel time to a series of reflecting horizons in the ice. We wish to relate these observations to changes in the mass of the ice body. To do this we make three basic assumptions:
-
1. The internal reflecting horizons we see in the radar record are material surfaces.
-
2. The rate of surface accumulation has been steady for at least as long as it has taken to form the present firn column, and firn compaction rates are such that the density–depth profile in any particular column would be steady if there were no divergence or convergence in the horizontal flow.
-
3. The volume, and hence the density, of any material element of ice or firn is unaffected by divergence or convergence in the horizontal flow, so the otherwise steady depth–density profile of a firn column is subject to vertical strain in the presence of a non-uniform horizontal velocity field.
The first of these will be true so long as the reflections we see arise from discrete discontinuities in the ice. We maximize the chance of this being the case through careful selection of the most prominent internal reflections. The latter two assumptions will not be true in detail, but are convenient approximations that greatly simplify the problem. Assumption (2) is equivalent to ‘Sorge's law’ of densification (Reference BaderBader, 1954), while assumption (3) allows us to apply Sorge's law to an ice mass experiencing divergent or convergent horizontal flow. Making assumption (2) restricts our analysis to dry accumulation zones, but these encompass most of the Antarctic ice sheet, including the majority of the ice shelves. Under assumptions (2) and (3) combined, the viscous response of the firn to horizontal stresses and its compaction under the weight of the overlying material are independent processes, the effects of which can simply be added together to give the net vertical strain of the firn. If the strain rate is steady, a steady depth–density profile will still result at any particular geographical point, but its precise form will depend on the strain rate. In the paper, we justify assumption (3) on theoretical grounds and show that our data are consistent with its implications.
We consider the evolution of a column of ice and firn experiencing arbitrary surface accumulation, horizontal flow and basal melt. Its location is defined with respect to a space-fixed Cartesian coordinate system (x, y, z), with the (x, y) axes horizontal and the z axis vertical and positive upwards. Particle motion with respect to these axes is denoted by velocity components (u, v, w) respectively. We identify three types of surface within the column (Fig. 15): the upper and lower boundaries, z = h s and h b; isochrones, z = ht , which link all material particles that were on the upper surface at the same instant in time; and isopycnals, z = h ρ , which link all material particles having the same density. We can only identify discrete isopycnal surfaces within the firn. The solid ice beneath is assumed to be of uniform density, ρi throughout.
Fig. 15. Motion of surfaces within a column of ice and firn subject to arbitrary surface accumulation, basal melt and vertical strain. Mass is added to or removed from the column through the upper and lower surfaces, h s and h b. Isochrones, h t, are material surfaces connecting particles that were on the upper surface at the same instant in time, t. They sink lower in the column as mass is added at the upper surface and come closer together as the firn density increases. Vertical strain modifies the rate at which isochrones approach each other, but does not alter the rate at which the firn density increases. The isopycnal surfaces, hρ , also converge or diverge as the ice column strains. We assume that this is the only process that affects the depth–density profile. In particular, firn compaction causes no change in the profile, thus leaving the depth of the isopycnal surfaces unaltered, and this implies that the mass flux across all isopycnals is the same.
For an arbitrary surface, z = h, we can write a kinematic boundary condition:
which expresses the rate at which material particles cross the surface (righthand side) as the difference between the vertical velocity of the surface (first three terms on the lefthand side) and the vertical velocity of the material particles that are instantaneously on the surface. The vertical component of the mass flux across the surface is 
and it has the same sign as w. For the upper and lower boundaries of the ice column z = h
s and h
b the mass flux across the surface is simply (minus) the surface accumulation rate and (minus) the basal melt rate respectively. For the isochrones the mass flux across the surface is zero. The vertical velocity of the isochrones is equal to that of the material particles that lie on them, and hence those same particles lie on the isochrones for all time. This is our definition of a material surface.
The isopycnal surfaces require special consideration. Any process that leaves the density of material particles unaltered will by definition leave those material particles on the same isopycnal surface. The isopycnal will then be a material surface. In contrast, any process that alters the density of material particles causes them to move from one isopycnal surface to another. There will then be a mass flux across the isopycnal surfaces, which will no longer be material surfaces. Assumption (2) allows us to write a simple expression for the mass flux across an isopycnal surface. If the density–depth profile is to be unchanged by accumulation, compaction or indeed any other process that alters the density of material particles, then the physical separation of two isopycnal surfaces must not be altered by flow across them. The physical separation is simply the mass between the surfaces divided by the density. Since the latter is fixed by definition, the total mass between any two isopycnals must remain unaltered by the mass fluxes across them. The only way to ensure this is for the mass flux across every isopycnal surface to be equal. The upper surface of the ice mass, z = h
s, is an isopycnal surface, with density ρs
, and the mass flux across this surface is . Therefore, the mass flux across all isopycnal surfaces must equal ms, and we can write down the kinematic boundary condition that must hold at each and every one:
Under assumption (3), there must be no additional mass flux across the isopycnal surfaces in response to convergence or divergence in the horizontal velocity field. Equation (A2) has the desired properties in that it places no restriction on the vertical motion of the isopycnal surface, which is free to move anywhere in (x, y, z) space with the ice flow, but ensures that the mass flux across it is always consistent with assumption (2). Further, when the surface accumulation rate is zero, the isopycnal surface is a material surface. Isopycnals and isochrones would then remain parallel and would move apart or come together in unison under the influence of the vertical strain rate. Hence, Equation (A2) is a general statement of assumptions (2) and (3).
Note that there is no requirement that the vertical mass flux be constant with depth in the ice column. All we have stated is that the mass flux across isopycnal surfaces, which are free to move relative to the upper boundary of the ice mass, is constant. However, Equation (A2) does imply that the motion of the isopycnal surfaces relative to the material particles (and thus relative to isochrones) is a function of density and accumulation rate only. Thus, if the accumulation rate is constant, each material particle undergoes the same history of densification, so firn of any given density will always be of the same age, although that density can occur at different depths below the upper surface at different horizontal locations.
The spatial or temporal gradient of a constant-density (isopycnal) surface is related to the density gradient along a geopotential (constant z) surface as follows:
Multiplying Equation (A2) by −(∂ρ/∂z) and substituting from Equation (A3), the expression for the kinematic boundary condition can be rewritten as:
This is our generalized version of Sorge’s law that incorporates assumptions (2) and (3).
We now use the principle of mass conservation,
to derive an expression for the true relative motion of two arbitrary radar horizons, z = h u and h l. In general, these could be either internal reflectors or the air–ice and ice–water boundaries at the top and bottom of the ice shelf. Using Equation (A4) we can rewrite the mass conservation equation (A5) as:
which simplifies to:
We integrate Equation (A7) between our arbitrary upper (z = h u) and lower (z = h l) surfaces to obtain:
Note that, although our primary focus in this paper is the study of ice shelves, we have made no assumption about the constancy of the horizontal velocity with depth. Our equations are thus generally applicable to any ice body. We now use the kinematic boundary conditions at the upper and lower surfaces:
to rewrite Equation (A8) as:
We next decompose the horizontal velocities into depthmean and depth-varying components, such that:
and the depth-varying (primed) components, by definition, integrate to zero between the upper and lower surfaces. This yields the desired expression for the temporal evolution of the true vertical separation of two surfaces in a moving column of ice:
The second and third terms on the lefthand side represent the contributions from vertical strain and densification respectively, while the final terms give the contributions resulting from the net gain or loss of mass.
In order to make use of Equation (A12) we must relate the true separation of the upper and lower surfaces (h u − h l) to the two-way radar travel time we observe. The travel time is twice the integral of the radar velocity over the distance between the reflectors:
where n is the refractive index of the firn and ice the radar signal passes through and c is the speed of light in vacuo. In dry firn we assume that the effective refractive index is simply related to the porosity, ϕ, by:
where n a and n i are the refractive indices of air and solid ice respectively (Reference Kovacs, Gow and MoreyKovacs and others, 1995). Since the density of the firn is given by an analogous expression,
and to a good approximation we can take n a = 1 and ρ a = 0, we can combine Equations (A14) and (A15) to give:
Now if we introduce a decomposition of the density into depth-mean and depth-varying components, analogous to the treatment of the velocities above, then substitute the expression (A16) for the refractive index into Equation (A13) and perform the integration, we arrive at:
This expression gives us the observed two-way travel time of a radar pulse travelling between two reflectors in terms of the true separation of the reflectors plus a correction term related to the total mass of ice and firn between them. The quantities that we measure with the phase-sensitive radar system are the travel time, T, and its temporal evolution, DT/Dt.
We now need an expression for the temporal evolution of the total mass between the upper and lower surfaces. To obtain this, we return to the mass conservation equation in its original form (A5) and depth-integrate:
Applying the kinematic boundary conditions, (A9) and (A10), we can rewrite this expression as:
and decomposing the velocities and density into their depth-mean and depth-varying components we arrive at:
Using (A17) to combine (A12) and (A20), we can form an equation for the mass balance of a column of ice in terms of the quantities that we measure with the radar plus some additional terms related to the variation of density with depth:
It is common practice to scale the measured travel time by half the radar velocity in solid ice to yield an equivalent ice thickness, H e = Tc/2n i. Introducing this scaling, we arrive at our final expression for the temporal evolution of the observed thickness:
This expression is valid for a completely arbitrary choice of upper and lower surfaces and arbitrary density–depth and velocity–depth profiles. For an ice column of constant density the only terms that would survive would be the first two (total derivative of thickness and contribution from horizontal divergence) and the fourth and fifth (net gain or loss of mass at the upper and lower surfaces). For the general case considered here, there are corrections to the observed thickness change for the case when the mass gained or lost at the upper and lower surfaces is not fully compacted solid ice (first two terms in the braces), densification of the firn between the reflectors by compaction (last term in the braces), and a correction for the case when depth variation of the horizontal velocity field leads to uneven expansion or contraction of parts of the ice column having different density, thus causing a change in the mean density of the column (last term on the lefthand side).
In most situations of practical importance, the above expression simplifies somewhat. For example, if the horizontal velocity is constant with depth, as is usually the case on an ice shelf, the last term on the lefthand side vanishes. The lower surface will normally be deep enough that the density there is that of solid ice, in which case the second term in the braces is zero. Conventional determinations of the mass balance use the top surface of the ice shelf as the upper reference surface, in which case the other two terms in the braces cancel. No correction is then required for the varying radar velocity in the firn, because the mean density of the ice and firn column remains unchanged. If the upper reference surface is an internal radar reflector, the mass flux across it is zero, and the first term in the braces and the fourth term (immediately following the braces) vanish. Then if the bottom surface lies in solid ice and there are no variations of velocity with depth, the equation is the same as that presented by Reference Corr, Jenkins, Nicholls and DoakeCorr and others (2002), except for the small correction factor 1/n i in front of the compaction term. If both upper and lower surfaces are internal horizons, the first two terms in the braces and the two terms following the braces are all zero. Temporal evolution of the observed thickness then only depends on vertical strain and compaction, provided there are no variations of velocity with depth. If all the internal reflectors used are deep enough that they lie in solid ice, all correction terms for varying density and varying velocity vanish, and we arrive at the simple version of the equation generally used in this paper.
