Satellite-derived volume loss rates and glacier speeds for the Cordillera Darwin Icefield

We produce the first icefield-wide volume change rate and glacier velocity estimates for the Cordillera Darwin Icefield (CDI), a 2605 km2 temperate icefield in southern Chile (69.6 W, 54.6 S). Velocities are measured from optical and radar imagery between 2001–2011. Thirty-six digital elevation models (DEMs) from ASTER and the SRTM DEM are stacked and a weighted linear regression is applied to elevations on a pixel-by-pixel basis to estimate volume change rates. The CDI lost mass at an average rate of −3.9± 1.5 Gt yr−1 between 2000 and 2011, equivalent to a sea level rise (SLR) of 0.01± 0.004 mm yr −1 and an area-averaged thinning rate of −1.5± 0.6 m w.e.(water equivalent) yr −1. Thinning is widespread, with concentrations near the front of two northern glaciers (Marinelli, Darwin) and one western (CDI-08) glacier. Thickening is apparent in the south, most notably over the advancing Garibaldi Glacier. The northeastern part of the CDI has an average thinning rate of −1.9± 0.7 m w.e. yr −1, while the southwestern part has an average thinning rate of −1.0± 0.4 m w.e. yr −1. Velocities are obtained over many of the CDI glaciers for the first time. We provide a repeat speed time series at the Marinelli Glacier. There we measure maximum front speeds of 7.5± 0.2 m day−1 in 2001, 9.5± 0.6 m day−1 in 2003 and 10± 0.3 m day−1 in 2011. The maintenance of high front speeds from 2001 to 2011 supports the hypothesis that Marinelli is in the retreat phase of the tidewater cycle, with dynamic thinning governed by the fjord bathymetry.


Introduction
We focus on the Cordillera Darwin Icefield (CDI), the third largest temperate icefield in the Southern Hemisphere (Bown et al., 2013), which along with the Northern Patagonian Icefield (NPI) and Southern Patagonian Icefield (SPI), has experienced a rapid reduction in ice-covered area (Rivera et al., 2007;Masiokas et al., 2009;Lopez et al., 2010;Davies and Glasser, 2012;Willis et al., 2012b).The CDI is located in the southernmost Andes (Fig. 1) in Tierra del Fuego.It is coalesced around two main mountain peaks, Mount Darwin (2469 m a.s.l., Koppes et al., 2009) and the nearby Mount Sarmiento (2300 m a.s.l., Strelin et al., 2008).The icefield covers 2605 km 2 , measured from ice outlines derived from satellite imagery acquired from 2001 to 2004.It extends roughly 200 km west-east from 71.8 • W to 68.5 • W and roughly 50 km south-north from 54.9 • S to 54.2 • S, boardered to the north by the Almirantazgo Fjord and the Beagle Channel in the south.Precipitation during the winter comes predominantly from the south/southwest (Holmlund and Fuenzalida, 1995), and the E-W orientation of the CDI leads to an orographic effect with greater snowfall on southern and western glaciers and drier, warmer conditions on northern and eastern glaciers (Holmlund and Fuenzalida, 1995;Strelin and Iturraspe, 2007;Koppes et al., 2009;Lopez et al., 2010).Previous length-change measurements of CDI glaciers show retreat at northern and eastern glaciers and stable/advancing fronts at southern and western glaciers (Holmlund and Fuenzalida, 1995).The "southern" part or side of the CDI refers to southern and western glaciers, and the Published by Copernicus Publications on behalf of the European Geosciences Union."northern" part or side of the CDI refers to northern and eastern glaciers.The purple line in Fig. 1 shows the divide between the northern (1322 km 2 ) and southern (1283 km 2 ) sides.
There are few studies on the CDI compared to other temperate ice fields (Masiokas et al., 2009;Lopez et al., 2010), such as the Alaskan ice fields, the NPI and the SPI (e.g., Arendt et al., 2002;Rignot et al., 2003;Berthier et al., 2010;Glasser et al., 2011;Ivins et al., 2011;Willis et al., 2012a).Climate and mass balance studies are scarce for southern hemispheric ice bodies outside of Antarctica (Holmlund and Fuenzalida, 1995;Lopez et al., 2010), due to the difficult access and weather.
Studies of observational and reanalysis data indicate the loss of ice at the CDI can be attributed to climatic changes that include 20th century regional decreases in precipitation (Quintana, 2004) coupled with atmospheric warming (Holmlund and Fuenzalida, 1995;Lopez et al., 2010) and dynamic instability at the largest glacier on the icefield, Marinelli (e.g., Holmlund and Fuenzalida, 1995).On a local scale, changes in wind patterns have increased precipitation on the southern side of the CDI (Holmlund and Fuenzalida, 1995;Strelin and Iturraspe, 2007), while decreasing precipitation on the northern side (Holmlund and Fuenzalida, 1995;Koppes et al., 2009).
Temperate ice fields are disproportionately large contributors to SLR (e.g., Arendt et al., 2002;Rignot et al., 2003); Rignot et al. (2003) claim this is particularly true of the Patagonian glaciers, which they say account for 9 % of the nonpolar contribution to SLR.The CDI, along with the NPI and SPI, provides an opportunity to examine the response of different glaciers (e.g., calving vs. noncalving) in different climates (maritime on the southern side versus more continental on the northern side) to regional changes in climate (Holmlund and Fuenzalida, 1995), and unlike the NPI and SPI the contribution of the CDI to SLR has not yet been estimated (Lopez et al., 2010).
The CDI is the closest icefield to the Antarctic Peninsula, a region that has also experienced significant warming.Mass loss at the CDI might be contaminating GRACE measurements of the Antarctic Peninsula, NPI and SPI (Ivins et al., 2011), so our constraints on the mass loss rate occurring at the CDI will help isolate this signal.Thinning and acceleration have been observed on glaciers in the Antarctic Peninsula and the NPI (Pritchard and Vaughan, 2007;Willis et al., 2012a), we assess whether this is the case for any glaciers on the CDI.
In this study we calculate both the elevation change rates ( dh dt ) over the entire CDI and measure glacier velocities using pixel-tracking applied to pairs of optical and radar images.With dh dt and an assumed density of material lost/gained, we can estimate the mass change rate, allowing us to quantify its SLR contribution and compare it with other ice fields.We also use the surface elevation change rates to identify which glaciers are providing the largest contribution to SLR and should be the focus of further study.Additionally, measuring glacier velocities allows an estimate of mass flux out of the glacier if the thickness is known.Increased speed and mass flux through the front of the glacier can cause "dynamic thinning" if it is not balanced by increased mass input.Our results will provide a baseline measurement over many glaciers and areas of the icefield for which ice velocities have not been measured.
The date format used throughout this paper is DD/MM/YYYY.

Data Preparation
The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) has a stereo-imaging capability, enabling DEMs to be generated on-demand by NASA's Land Processes Distributed Active Archive Center (LP DAAC) (Fujisada et al., 2005).In this study, ASTER DEMs (product 14) are used to calculate dh dt , while band 3N images (product 1B) are used for pixel-tracking.Clouds are not masked during data preparation, rather, they are largely removed by excluding elevations based on deviation from the the SRTM DEM (Sect.2.2).NASA's Automatic Registration and Orthorectification Package (AROP; Gao et al., 2009) is used to co-register ASTER images and DEMs to a Landsat GLS image (available from the Global Land Cover Facility) and orthorectify the ASTER L1B images using the Shuttle Radar Topography Mission (SRTM) DEM (acquired in 2000).Landsat GLS images are orthorectified to the SRTM DEM (Tucker et al., 2004); therefore, co-registering the ASTER imagery to the Landsat GLS image effectively coregisters them to the SRTM DEM (see Willis et al., 2012a for details).The SRTM DEM is our default DEM for orthorectifying the ASTER imagery; it is not as adversely affected by clouds, snow and other features that cause errors in optically derived (e.g., ASTER) DEMs (Scherler et al., 2008)

Elevation change rates
Horizontally co-registered ASTER DEMs are subsequently vertically co-registered and a weighted linear regression is applied to calculate dh dt for each pixel.In all, 36 ASTER DEMs (derived from imagery acquired from 2001 to 2011) and the SRTM DEM (acquired in February 2000) are processed, with an average of 4-5 elevations per pixel incorporated into the regression (Fig. 2 shows elevations and regression lines for several points on Marinelli Glacier).Each elevation is weighted by the inverse of the standard deviation of the bedrock elevation differences between its ASTER DEM and the SRTM DEM.This is a common measure of the uncertainty associated with ASTER DEM elevations (Fujisada et al., 2005;Howat et al., 2008).Horizontal misalignment will appear as off-ice ("bedrock") elevation differences when comparing ASTER DEMs to the SRTM DEM, and is therefore included in the standard deviation of the bedrock elevation differences (i.e., the greater the misalignment, the greater the standard deviation, which is the measure of un-certainty that we use for each ASTER DEM).We typically find values from 8 m to 20 m as our uncertainty for ASTER DEMs, similar to other studies (Kääb, 2002;Fujisada et al., 2005;San and Süzen, 2005;Rivera et al., 2007;Howat et al., 2008;Berthier et al., 2010).
In the ablation zone, we exclude elevations from the regression that deviate more than +5/−30 m yr −1 from the first elevation, in the accumulation zone the cutoff is +5/−10 m yr −1 .These cutoffs largely remove the areas of erroneously high elevations due to clouds, as well as outliers caused by shadow, snow and other sources which either obscure the ice or are largely featureless at ASTER resolution.It is important to note that although these cutoffs are based on expected dh dt , they are imposed on each individual map point (see Fig. 2) before applying the weighted linear regression, rather than to already calculated dh dt , which would yield drastically different results.The first elevation in our time series is SRTM for 94 % of the pixels, for the remaining areas we use the first available elevation from ASTER.The average end date is 17/05/2008.The high percentage of points with SRTM elevations makes our cutoff strategy feasible.The radar-derived SRTM DEM is not influenced by clouds like ASTER, making SRTM elevations a reliable starting point.
We limit the maximum negative deviation to −30 m yr −1 in the ablation zone to just capture the maximum thinning and exclude spurious elevations from the regression calculations.Changing this to −60 m yr −1 has no affect on the zones of maximum thinning.Less thinning is expected in the accumulation zone than in the ablation zone; therefore, we allow a maximum negative deviation of only −10 m yr −1 above the ELA.Permitting a greater negative deviation in the accumulation zone only results in incoherent areas of extreme and unrealistic thinning.
In order to determine the ELA and distinguish where to impose our different constraints, we examine late season (January and February) ASTER images.We record the altitude of the transition between bare ice and snow and take this as a rough proxy for the ELA (e.g., Bamber and Rivera, 2007).We are able to measure the ELA for 11 larger glaciers, for the rest of the icefield, we assume a regional ELA of 1090 m (Strelin and Iturraspe, 2007).
The volume change rate at each pixel is the dh dt for the pixel multiplied by the area of the pixel.Summing together the volume change rate of every pixel yields a volume change rate ( dV dt ) for the entire icefield.dV dt is multiplied by the density of glacier ice, set to 900 kg m −3 (e.g., Cuffey and Paterson, 2010) to produce an estimate of the mass change rate.Future ground-based measurements of densities on the CDI will be needed to find more accurate and precise values.
We have a sufficient number of ASTER DEMs to provide dh dt rates for 96 % of the total area of the icefield.Each pixel in the remaining small gaps is filled with the median dh dt value within 1 km.This local method is more accurate than using  the average dh dt of an elevation bin to fill gaps within that bin, or filling in gaps with the mean of all dh dt calculated (e.g., Rignot et al., 2003).For example, the average low-elevation dh dt and overall average dh dt are both negative, these rates would not be suitable for filling a gap near the front of Garibaldi Glacier, which is known to be advancing and where we measure positive dh dt at low elevations.
Penetration of C-band radar into ice and (particularly) snow (e.g., Rignot et al., 2001) is a potential problem when using the SRTM DEM.We use a technique pioneered by Gardelle et al. (2012) and applied by Willis et al. (2012b) to compensate for penetration effects.Willis et al. (2012b) compared X-band SRTM elevations (which have negligible penetration) with C-band SRTM elevations and found approximately 2 m of C-band penetration over the SPI at all elevations.Due to a lack of X-band SRTM coverage over the CDI, we cannot provide a similar analysis here.Instead, we assume the penetration depth is similar to that at SPI and add 2 m to each SRTM elevation over ice, which increases our mass loss rate by about 13 %.We note the CDI is colder than the SPI, which could lead to drier conditions and greater penetration (Rignot et al., 2001), but we do not have adequate X-band data to quantify the difference.
An additional source of mass change that we consider is sub-aqueous mass loss/gain (which does not contribute to sea level rise).We cannot directly measure sub-aqueous ice gains or losses, but make a rough estimate by measuring changes in the areal extent at the front of several glaciers that have undergone relatively large advance or retreat (from image pairs) and multiplying these area changes by the approximate depth below water of the ice.We assume an average depth below water of 150 m for Marinelli Glacier (see Koppes et al., 2009, Fig. 4a and b) and 60 m for the other glaciers, with an uncertainty of ±50 m (the uncertainty on the change in area of each glacier is negligible).Dividing the sub-aqueous volume change by the time interval separating the images used to find area change gives a rate.This calculation shows that subaqueous mass loss, while not well constrained in this study, is an order of magnitude lower than the overall mass change rate calculated from dh dt (< 5 %).

Sources of uncertainty
Below we consider sources of uncertainty on our mass change rate.These include the uncertainty on the elevations incoporated into the regression, uncertainty on the ELA, the effect of varying the maximum deviation allowed from the first elevation, different density scenarios, and uncertainty on the penetration depth of the C-band SRTM DEM.The uncertainty associated with the dh dt for each pixel is calculated from the model covariance matrix (e.g., Aster et al., 2005), which accounts for the uncertainties on the elevations incorporated into the regression.The 95 % confidence interval for the volume change rate uncertainty is calculated using the formula: . U is the total "volume" of uncertainty, calculated by taking the uncertainty at each pixel, multiplying it by the area of the pixel (to determine a "volume" of uncertainty for that pixel), and then adding together "volume" of uncertainty for each pixel where a dh dt is calculated.N is the number of independent pixels (e.g., Howat et al., 2008), which we determine by dividing the total area by the area over which off-ice dh dt are correlated (e.g., Rolstad et al., 2009).We estimate the scale at which the dh dt are independent by finding the area at which the variance of the off-ice dh dt begins to "flatten" (see Rolstad et al., 2009 andWillis et al., 2012a for details on the method), which we estimate to be 1800 m by 1800 m (Fig. 4).This is analogous to the "corner" point on an L curve (e.g., Aster et al., 2005, p. 91, Fig. 5.2) and indicates the lengthscale past which the dh dt are no longer correlated.The total contribution from the uncertainty on individual elevations is ∼ 0.35 Gt yr −1 .
The regional ELA is poorly known.In order to investigate the impact of changing the regional ELA on the mass loss rates, we lower our regional ELA from 1090 to 650 m, an ELA that has been found for several glaciers on the southern and western regions of the CDI (Bown et al., 2013).We also assume an uncertainty of ±200 m on the 11 glaciers where we have estimated the ELA from optical imagery.The difference between the two ELA scenarios is incorporated into our uncertainties (∼ 0.37 Gt yr −1 ).
The deviation allowed from the first elevation has a large impact on the mass loss rate (e.g., Willis et al., 2012b).The rate produced by allowing +10 m yr −1 is contaminated by low clouds but is given as a rough minimum estimate of the mass loss rate.The unsymmetric cutoff (+5 m yr −1 versus −10/−30 m yr −1 ) may bias our results towards thinning but we argue that an unsymmetric cutoff is more physically justified than a symmetric cutoff.Rapid retreat has been independently observed at Marinelli Glacier (Holmlund and Fuenzalida, 1995;Koppes et al., 2009), from which we know that a large amount of thinning must be occurring.The maximum allowed negative deviation from the first elevation of −30 m yr −1 just captures the maximum thinning at Marinelli Glacier (see Fig. 2, point 1).
The cutoff of +5 m yr −1 from the first elevation is based on the limited precipitation data available for this region.The upper limit on precipitation in Patagonia is about 10 m yr −1 (e.g., Holmlund and Fuenzalida, 1995;Rignot et al., 2003;Rasmussen et al., 2007;Koppes et al., 2011), and PRECIS model results show maximum precipitation in the 4000 to 6000 mm yr −1 range for the CDI (Fernandez et al., 2011).One year of accumulated precipitation undergoes densification into firn and eventually ice and provides much less than +10 m yr −1 of elevation change through time, so +5 m yr −1 is chosen as a reasonable upper limit on the maximum thickening expected over the large areas covered in this study.We note that a point measurement may yield a dh dt higher than +5 m yr −1 but we would not expect a single point to be representative of sustained thickening rates over square kilometers.Figure 2 highlights the importance of applying this cutoff, cloud-influenced elevations would seriously degrade the quality of our regressions for each pixel.Tests using a higher positive cutoff produced discontinuous and incoherent "splotches" of extreme positive dh dt that are unrealistic.The only exception is the lower ablation zone of Garibaldi Glacier, where we allow a maximum positive deviation of +10 m yr −1 to accommodate the known advance of Garibaldi (Fig. 3).
The sub-panels in Fig. 2 illustrate the problem of using a positive cutoff of +30 m yr −1 , which would lead to the inclusion of ASTER elevations that are obvious outliers.ASTER elevations in the accumulation zone are generally less reliable due to greater cloud and snow cover (lack of contrast).
www.the-cryosphere.net/7/823/2013/The Cryosphere, 7, 823-839, 2013 Errors due to clouds tend to be positive rather than negative; this is confirmed by comparing the average off-ice dh dt of increasing symmetric maximum allowed deviations.Allowing ±5, ±10 and ±30 m yr −1 results in average dh dt of −0.04 m yr −1 , +0.10 m yr −1 and +0.82 m yr −1 .Allowing a deviation of −30 m yr −1 is necessary to capture maximum thinning, but 330 m of thickening (+30 m yr −1 from 2000 to 2011) is unreasonable in both the ablation and accumulation zone for this region (Koppes et al., 2009;Lopez et al., 2010).
In situ measurements of accumulation rate on the CDI are required to refine our estimates further; the cutoffs we use are the best available based on the literature (e.g., Fernandez et al., 2011).To assess the uncertainty due to the choice of maximum deviation allowed, we find the mode (most frequent occurrence, i.e., peak) of the distribution of elevation differences between all ASTER DEMs and the SRTM DEM (normalized by dividing by the time interval between each ASTER DEM and the SRTM DEM).The mode is not affected by the choice of allowed deviation, and so is an independent measure that we can compare our regression-derived rates to.
We do not present the "mode rate" as representing the best estimate of the average dh dt .
It is approximately what the rate would be if the distribution were Gaussian around the peak dh dt .This implicitly assumes that rates that are equidistant from the mode rate are equally likely (e.g., if the mode rate is −1.5 m yr −1 ) then rates of +18.5 m yr −1 are approximately as common as rates of −20.5 m yr −1 .We consider 1800 m Data (variance of o -ice dh/dt at corresponding length scale on x-axis) Fit Fig. 4. A plot of pixel size (listed as side length of square pixel) vs. variance of bedrock dh dt (calculated using a −10/+10 m yr −1 cutoff), plotted as red dots.A linear fit to log(x)/y shows the trend of the data (blue line).The black circle indicates the point at which the curve is "flattening", i.e., the lengthscale past which the variance no longer changes significantly, indicating that the dh dt are no longer correlated.1800 m is selected as a conservative estimate of the decorrelation length for the dh dt .
this unlikely given the coherent thinning of approximately −25 m yr −1 estimated at the front of several glaciers, and the lack of evidence for a similar degree of coherent thickening anywhere on the icefield.Therefore, the mode rate is probably lower than the the actual thinning rate, and the difference between the mode rate and our regression-derived dh dt is an overestimate of the maximum "bias" due to our choice of dh dt .We add the difference (±1.22 Gt yr −1 ) that results from using the same density and penetration assumptions as our regression rate to our uncertainties to fully account for any possible bias.
We assume all volume change is at a density of 900 kg m −3 , consistent with previous studies (e.g., Rignot et al., 2003;Berthier et al., 2004Berthier et al., , 2010)).However, the density of the lost material is likely variable.We consider two different density scenarios.Assuming that firn (600 kg m −3 ) is lost in the accumulation zone and ice (900 kg m −3 ) is lost in the ablation zone (considered as a possible scenario by Kääb et al., 2012) reduces the mass change rate by 0.3 Gt yr −1 .Assuming an ice density of 900 ± 25 kg m −3 (e.g., Gardner et al., 2012) adds ±0.1 Gt yr −1 to our uncertainty.The effect of density changes are small (∼ 10 %) as most of the mass loss occurs from ice in the ablation zone.We incorporate the 600/900 kg m 3 scenario into our uncertainties, as it adds a further ±0.3Gt yr −1 and is thus a more conservative estimate of uncertainty.
The penetration depth of the C-band SRTM into snow and ice is an additional source of uncertainty.Assuming a range of ±2 for the penetration depth yields a variation of +0.47 Gt yr −1 for 0 m of penetration and −0.48 Gt yr −1 for 4 m.We add ±0.48 Gt yr −1 to our uncertainties to account for penetration depth uncertainty.
Taking the square root of the sum of the squared uncertainties gives us an overall uncertainty of ±1.5 Gt yr −1 .Uncertainties beyond the statistical uncertainty from the regression are added as an average rate to the uncertainties for individual glacier basins in Table 1.

Velocities
In this section we describe how we use data from three different satellites with different resolutions to measure sub-pixel offsets and convert them to glacier velocities.

ASTER
Sub-pixel offsets between ASTER image pairs (pixel resolution of 15 m/pixel) are measured via normalized amplitude cross correlation, with a spatial resolution of 120 m (i.e., a step size between cross correlations of 8 pixels, see Willis et al., 2012a for details).This technique, known as "pixeltracking", has been used to track velocities on many glaciers (e.g., Scambos et al., 1992).
AROP is used to co-register the more recent scene in a pair of orthorectified ASTER images to the earlier scene to minimize misfits.The open source ROI PAC's "ampcor" routine (Rosen et al., 2004) is used to calculate E-W and N-S offsets.The results are post-filtered by excluding offsets with a signal-to-noise ratio (SNR, which is the peak cross correlation value divided by the average) below a manually selected threshold (Willis et al., 2012a).A linear elevation-dependent correction (determined from apparent "bedrock" velocities) is applied to the velocities to correct for the elevation-dependent bias due to imprecise coregistration/orthorectification (Nuth and Kääb, 2011).Figure 5 shows a typical trend for an ASTER pair.We do not know the physical origin of the bias, Ahn and Howat (2011) attribute the systematic elevation-dependent displacement errors over the ice surface to co-registration errors.Applying an elevation-dependent correction based on the displacement of off-ice areas largely removes this (Fig. 6).
Uncertainty for each pair is estimated from motion on ice-adjacent "bedrock" (see Willis et al., 2012a for details), which should be zero.Horizontal misalignment leads to office motions when calculating offsets, so error due to misalignment is included in the uncertainty estimate.

Advanced Land Observing Satellite (ALOS)
We use L-band SAR pixel-tracking (e.g., Rignot, 2008;Strozzi et al., 2008;Rignot et al., 2011)   routine.ALOS interferometry does not yield velocities due to the relatively large motions, long separation between scenes, and changing surface characteristics.

Orthorectification errors from thinning
Pronounced thinning can cause orthorectification errors in optical imagery that do not appear in the off-ice velocities.
The longer the timespan between the acquisition date of imagery used to measure offsets and the DEM used to orthorectify the imagery, the greater the error.Our standard processing uses the SRTM DEM, acquired in February 2000, to orthorectify the ASTER L1B images used for pixel-tracking.
Examining our dh dt results (Fig. 1) reveals this to be a potential problem over the three most rapidly thinning glaciers: Marinelli, CDI-08 and Darwin.For Marinelli Glacier we mitigate this effect by orthorectifying QuickBird 2 imagery from 2011 to a 2007 ASTER DEM rather than the SRTM DEM.Furthermore, the difference in incidence angle between the two QuickBird 2 images is less than 0.2 degrees.This means the base/height (B/H ) ratio (e.g., Fujisada et al., 2005)  For this pair we orthorectify to the SRTM DEM as the time between SRTM acquisition and the ASTER pair is fairly short.The direction of the velocity vectors are consistent with glacier flowlines visible in the ASTER imagery providing confidence in the result.Velocities over CDI-08 are from radar pixel-tracking, orthorectification errors are not a concern because radar images are not orthorectified.

Marinelli Glacier -Flux
We calculate flux along transects perpendicular to glacier flow (as close as possible to the front) for velocities from ).The 25/09/2001 and 13/09/2003 ASTER DEMs do not need to be adjusted because they are coincident with one of the images in the velocity pairs they are being used to orthorectify, and the maximum interval is 18 days (for the 2001 pair).We assume an average glacier depth below water of 150 m (see Koppes et al., 2009, Fig. 4a and b).Adding this to the height gives an approximate thickness.We multiply the glacier thickness by the perpendicular velocity along the transect to calculate flux.Sources of uncertainty that we include in the uncertainty for our flux estimates are the uncertainties on the speed, uncertainty on the depth below water (±50 m), and uncertainty on the DEMs used to obtain elevations.

Elevation change rates
Figure 1 provides a map of dh dt for the CDI.Thinning is concentrated at lower elevations at the front of several tidewater glaciers.Thickening is apparent on the southern side of some of the highest mountains and at the front of the advancing Garibaldi Glacier.We calculate a dV dt of −4.2 ± 1.7 km 3 yr −1 for the CDI (extent shown in Fig. 1).This rate equates to an area averaged dh dt of −1.6 ± 0.7 m yr −1 .Our estimate of the sub-aqueous volume loss rate (−0.12 ± 0.06 km 3 yr −1 ) provides an additional (∼ 3 %) loss, giving a final dV dt of −4.3 ± 1.7 km 3 yr −1 , equivalent to a mass loss rate of −3.9 ± 1.5 Gt yr −1 , assuming surface elevation changes are not due to densification processes and assuming an ice density of 900 kg m 3 .This corresponds to an overall averaged thinning rate of −1.5 ± 0.6 m w.e.yr −1 .The ablation zone, which comprises 60 % of the icefield, accounts for 76 % of the mass loss.The ablation zone has an average thinning rate of −1.8 ± 0.6 m w.e.yr −1 (including sub-aqueous mass change), the accumulation zone has an average rate of −0.9 ± 0.3 m w.e.yr −1 .Three glaciers (Marinelli Glacier, Darwin Glacier, and CDI-08 Glacier) account for 31 % of the mass loss, but cover only 12 % of the icefield area.
We estimate a thinning rate of −1.0 ± 0.4 m w.e.yr −1 for the southern side, which is lower than the northern side (−1.9 ± 0.7 m w.e.yr −1 ).The contrast between north and south is most likely due to warming in the north and changing weather patterns that have increased precipitation on the windward side of the mountains and decreased precipitation on the lee side (e.g., Holmlund and Fuenzalida, 1995), and of course the combination of climate with dynamic instability at glaciers such as Marinelli (e.g., Holmlund and Fuenzalida, 1995;Porter and Santana, 2003;Koppes et al., 2009).Table 1 gives the dV dt for the 16 largest glaciers, with rates for the entire basin, the accumulation zone and the ablation zone.

Velocities
Pixel-tracking provides useful velocities from twenty ASTER image pairs, one QB02 pair and two ALOS pairs with acquisition dates between August 2001 and August 2011 (Fig. 7).ASTER pixel-tracking results are generally better over the northern half of the CDI as there is less cloud cover than to the south (Holmlund and Fuenzalida, 1995;Strelin and Iturraspe, 2007).Composite speed results are shown in Fig. 8, velocities for individual pairs over selected glaciers are shown in Fig. 9. Average front speeds for the tidewater Marinelli Glacier (133 km 2 -the largest glacier of the CDI, e.g., Koppes et al., 2009)  Speeds at the front of the tidewater Darwin Glacier (46 km 2 ) reach a maximum average of 9.7 ± 0.8 m day −1 for the period 25/09/2001 to 02/10/2001.No repeat measurements of motion from ASTER pairs are available for the Darwin Glacier.CDI-08 (127 km 2 ), the furthest west and south of the three most rapidly thinning glaciers, reaches speeds of 2.0 ± 0.5 m day −1 within 1 km of its 15/01/2011 front (04/01/2011 to 19/02/2011 ALOS pair), unfortunately repeat speeds are not available for this glacier.

Individual glaciers
A longitudinal profile of elevation and thinning rates for Marinelli Glacier (Fig. 11) illustrates both terminus retreat and the dramatic thinning occurring at the glacier.Thinning extends from the terminus to the highest parts of the accumulation zone.The longitudinal profile of the glacier surface remains convex as the glacier thins and retreats (Fig. 11), indicating that basal stresses are relatively high and the front remains grounded (e.g., Koppes et al., 2009) as of 13/11/2007.
The Marinelli Glacier retreated ∼ 4 km between 2001 and 2011 (measured from ASTER and QB02 imagery), an average rate of retreat of 0.4 km yr −1 .This is lower than the retreat rate of 1 km yr −1 measured by Koppes et al. (2009) during the late 1990s.Our rate is similar to their average rate of ∼ 0.3 km yr −1 from 1960 to 2005.Koppes et al. (2009) infer volume loss rates at Marinelli Glacier from observed retreat rates.A higher retreat rate provides a higher volume loss estimate.They find that the volume change rate dropped from −0.7 km 3 yr −1 in 1997 to −0.2 km 3 yr −1 by 2005.In contrast, our measurements show that volume loss at Marinelli is sustained at an average rate of −0.7 ± 0.2 km 3 yr −1 until at least 2007.
The accumulation area ratio (AAR) of Marinelli Glacier is 0.38 at an ELA of 1100 m (Fig. 12b).Marinelli underwent a significant reduction in AAR as the ELA moved up to its current altitude, likely contributing to the period of neg-ative mass balance in the mid-1990s preceding the current phase of rapid retreat (e.g., Koppes et al., 2009).The AAR has already moved through the part of the hypsometry curve where small changes in elevation lead to large changes in the AAR, an increase of 200 m from the current ELA of 1100 m would only reduce the AAR from 0.38 to 0.31.The negative mass balance operates in conjunction with dynamic instability; now that large changes in the ELA have relatively little affect on the AAR, the continued recession of Marinelli Glacier will likely be even more dependent on the underlying fjord bathymetry.
CDI-08, the third glacier with a strong thinning signal, is on the southern side of the icefield and faces west.This position and aspect favors increased snowfall in the prevailing climate (Holmlund and Fuenzalida, 1995).Given its location and the direction it is facing, CDI-08 is somewhat anomalous in that almost 2 km of retreat has occurred and the glacier has thinned rapidly between 2001 and 2011.The average thinning rate for CDI-08 is −3.0 ± 0.9 m w.e.yr −1 , compared to −1.0 ± 0.4 m w.e.yr −1 for the southern part of the CDI.
Garibaldi Glacier, on the southern side of the icefield, with a southern aspect, by contrast, has relatively large areas of positive dh dt in its ablation zone.Our dh dt here are based on sparse temporal coverage, and there is only one ASTER image (13/09/2003) that covers the entire front of Garibaldi Glacier.Comparison of a WorldView-1 (optical) image from 09/27/2011 with a 13/09/2003 ASTER image shows that the glacier has advanced by more than 1 km between 2003 and 2011.Figure 3 shows the frontal variation history of Garibaldi over a ∼ 65 yr period from Landsat TM and aerial imagery, and also shows advance in the past decade.Further confirmation of advance in the past decade comes from field reports in the austral summer of 2007, when the glacier was observed destroying trees and frequently calving (Masiokas et al., 2009).This gives us confidence in the positive dh dt we observe for Garibaldi Glacier, despite the limited data.The maximum recorded retreat was reached in 2001, followed by advance to the 2011 front that is within 0.5 km of the 1945 front (the maximum extent Garibaldi during this period), and no more than 2 to 3 km of retreat/advance in the past 65 yr.
The contrast in current behavior between CDI-08 and Garibaldi may be due partly to the lower ELA of CDI-08 and the higher sensitivity of its AAR to changes in climate compared with Garibaldi (Fig. 12a and c show the AAR versus   Post et al. (2011) note that an AAR above 0.8 typically favors an advance in TWG, assuming the ELA of Garibaldi is lower than 900 m, its AAR is greater than 0.8.The behavior of Garibaldi Glacier, with an oscillating, quasi-stable front since at least 1945 is consistent with an AAR that hovers around 0.8 and is relatively insensitive to rising temperatures in the region (e.g., Holmlund and Fuenzalida, 1995;Strelin and Iturraspe, 2007;Lopez et al., 2010).
Our results (positive dh dt and advance in the ablation zone versus relatively more negative dh dt in the accumulation zone) indicate the possibility of surge-like behavior (e.g., Rivera et al., 1997).However, a time series of velocities for Garibaldi Glacier would be necessary to ascertain whether this is the case.It is important to note that the underlying topography may be playing an important role in regulating Garibaldi's behavior as well, e.g., buoyancy forces at the front may be reduced by a submarine terminal moraine and/or a shallower fjord (e.g., Post et al., 2011).

Comparison with other southern hemispheric glaciers
Retreat and thinning is also occurring on other glaciers at similar latitudes in the Southern Hemisphere.Ice masses in the Kerguelen Islands (49 • S, 69 • E), for example, ex-perienced accelerating retreat rates and an average dh dt of −1.4 to −1.7 m yr −1 between 1963 and 2000 (Berthier et al., 2009).Glaciers on the subantarctic island of South Georgia (54.5 • S, 37 • W) have undergone retreat in response to warming since the 1950s, with recession being particularly pronounced from the 1980s to the mid-2000s (Gordon et al., 2008).Widespread retreat is also observed for the glaciers of Heard Island (53.1 • S, 73.5 • E) in the southern Indian Ocean (Thost and Truffer, 2008).
A similar amount of warming has occurred in all these locations, with an increase in average temperature of 0.4 to 1.4 • C south of 46 • S in Patagonia since the beginning of the 20th century (Rosenblüth et al., 1995); 1 • C of warming between 1964 and 1982 for the Kerguelen Islands (e.g., Berthier et al., 2009); approximately 1 • C of warming from 1950 to mid-1990s on South Georgia (e.g., Gordon et al., 2008); and 0.9 • C of warming for Heard Island since the mid-20th century (e.g., Thost and Truffer, 2008).
At the CDI and SPI, the thinning and rapid retreat of Marinelli and Jorge Montt glaciers from topographic pinning points (Koppes et al., 2009;Rivera et al., 2012) anticipates the predicted response of tidewater glaciers in other regions, such as South Georgia, if warming trends continue (e.g., Gordon et al., 2008, Sect. 5.1).

Marinelli Glacier -overview
Below, we estimate flux for Marinelli Glacier using our speeds, then compare our results with Koppes et al. (2009), who estimate the terminus speed and flux of Marinelli from the retreat rate.We find that our results do not agree with Koppes et al. (2009).Whereas they infer a reduction in terminus speed for Marinelli from 2000 to 2005, we find that the front speed in 2003 and 2011 is at least as high as in 2001, and consequently the flux in 2011 is approximately the same as the flux in 2001.We conclude that thinning at Marinelli Glacier is probably dynamic, with bed geometry likely governing velocity and retreat.We then consider Marinelli Glacier as a tidewater-cycle glacier (TWG) in retreat phase (e.g., Meier and Post, 1987;Motyka et al., 2003;Post et al., 2011), and compare it with Jorge Montt Glacier on the SPI.

Marinelli Glacier -flux
We estimate a flux of 0.5 ± 0.2 km 3 yr −1 for the 2001 pair, 0.7 ± 0.2 km 3 yr −1 for the 2003 pair and 0.5 ± 0.2 km 3 yr −1 for the 2011 pair.Flux is highest in 2003 (due to higher speeds than 2001 and a larger front than 2011), but the important point is that the 2011 flux has not dropped relative to 2001 due to speeds at the 2011 front that are higher than 2001 and as high as 2003.Seasonal changes in conditions at Marinelli Glacier could be influencing our results, we note that the two ASTER pairs from which we obtain front speeds are both from September.However, it is possible that the speeds in September 2001 and 2003 are different due to inter-annual variations in the onset of conditions affecting the glacier speed.
The 2011 QuickBird 2 pair is from 30/07/2011 to 16/08/2011.While this is not an entirely different season from the ASTER pairs, it is a month earlier, which is long enough for seasonal variations to possibly play a role in any observed speed differences.We consider it unlikely, however, that a seasonal component of motion is dominant in the 30 % increase in the front speed at Marinelli between September 2001 (7.5 m day −1 ) and August 2011 (10 m day −1 ).Unfortunately, we do not have enough repeat measurements to quantify any seasonal effect on speeds, especially given that our ALOS results over Marinelli from the austral summer do not reach the front.et al. (2009) infer speeds from the retreat rate, with lower retreat rate leading to lower inferred speed.They show Marinelli glacier slowing from a terminus ice speed of 8 m day −1 in 2001 (similar to our 7.5 m day −1 2001 speed) to 5.5 m day −1 in 2003.Our results, however, are evidence the glacier has not slowed from 2001 (7.5 m day −1 ) to 2003 (9.5 m day −1 ).We find that, contrary to the conclusion of Koppes et al. (2009), the front speed and therefore the flux are not decreasing from 2001 to 2003.Furthermore, the front speed in 2011 is higher than 2001, so despite the front being less extensive our estimate of the 2011 flux is as high as for 2001.

Marinelli Glacier -thinning gradient maintains surface slope at the front
The maintenance of relatively high speeds at the front of Marinelli Glacier between 2001 and 2011 can be attributed, in part, to the observed gradient in dh dt .There is rapid thinning and retreat at the front, with slower thinning rates upstream.This imbalance maintains surface slope near the glacier front, sustaining the driving stress according to the equation σ = ρgh sin(α), where σ is the driving stress, ρ is the density, g is gravitational acceleration, h is thickness and α is the slope (e.g., Cuffey and Paterson, 2010, p. 295).

Marinelli Glacier -tidewater cycle
The rapid retreat noted in this study and covered elsewhere (e.g., Koppes et al., 2009;Warren and Aniya, 1999), coupled with the thinning we observe between 2000 and 2011, suggests that Marinelli Glacier is a tidewater-cycle glacier (TWG) in retreat phase (e.g., Meier and Post, 1987;Motyka et al., 2003;Post et al., 2011).Front recession has opened the fjord at Marinelli.Around 1945 retreat began as the glacier receded into deeper water from the arcuate terminal moraine it had been pinned at, with particularly rapid retreat throughout the 1990s and 2000s (Holmlund and Fuenzalida, 1995;Porter and Santana, 2003;Koppes et al., 2009).Fjord bathymetry (and possibly ocean temperatures) likely governs the ice dynamics and calving rate (e.g., Benn et al., 2007;Koppes et al., 2009;Straneo et al., 2010), which probably control thinning at Marinelli (e.g., Koppes et al., 2009).The velocity results for Marinelli Glacier support this hypothesis.The 2011 speeds at Marinelli, which are higher than 2001 and as high as 2003, suggest that the glacier has not yet retreated to fjord depths shallow enough to slow it down.Along with an estimate of the 2011 flux that is as high as 2001, this leads us to predict that Marinelli Glacier will continue to retreat at least until it reaches a new topographic pinning point, which would likely cause the terminus speed to fall and calving rates to drop.There is no evidence of accumulation in the recent history of Marinelli Glacier that would compensate the mass loss due to calving.NCEP-NCAR climate model results show that precipitation has decreased at Marinelli Glacier and temperatures have risen (Koppes et al., 2009).
We do not see similar large thinning signals at adjacent glaciers with similar settings, which is further evidence that dynamic instability, and not melt due to a regional increase in atmospheric temperatures, is the main proximate cause of the current thinning at Marinelli.However, warming and decreased precipitation has likely helped trigger the initial retreat and thinning from the terminal moraine and contributed to the recent history of negative mass balance (Holmlund and Fuenzalida, 1995;Porter and Santana, 2003;Koppes et al., 2009).
Jorge Montt Glacier on the SPI is a somewhat analogous TWG, with a grounded front, thinning occurring at approximately the same rate, and a rapidly receding terminus (Rivera et al., 2012;Willis et al., 2012b).As at Marinelli, terminus retreat opened the fjord at Jorge Montt, with subsequent rapid retreat as the glacier receded into deeper water (Rivera et al., 2012).At Marinelli the retreat is thought to have been initiated by thinning linked to climate changes (e.g., Koppes et al., 2009).Though there is a high degree of variability from glacier to glacier on the Patagonian ice fields (Rivera et al., 2012) the similar behavior of these two glaciers in the same region but on different ice fields suggests that changing climate plays a role in retreat and thinning at Jorge Montt as well, given that the retreat was inferred to have been triggered by warming-induced thinning at the fronts of both glaciers.

Conclusions
We provide the first icefield-wide thinning rates and glacier velocities for the Cordillera Darwin Icefield.We find thinning at an average rate of −1.5 ± 0.6 m w.e.yr −1 for the CDI between 2000 and 2011 (equivalent to a mass loss rate of −3.9 ± 1.5 Gt yr −1 ) with less thinning in the south (−1.0 ± 0.4 m w.e.yr −1 ) compared to the northern part of the icefield (−1.9 ± 0.7 m w.e.yr −1 ).This pattern of thinning is consistent with climate records that show a warming trend in this region from the 1940s through the 1990s, leading to higher temperatures along with decreased precipitation on the northern side of the icefield and an increase in precipitation on the southern side (Holmlund and Fuenzalida, 1995;Strelin and Iturraspe, 2007;Koppes et al., 2009).
Operating in conjunction with local glacier geometry, rising temperatures triggered a period of prolonged retreat at Marinelli Glacier (since 1945) (e.g., Holmlund and Fuenzalida, 1995;Porter and Santana, 2003;Koppes et al., 2009).Recently, Marinelli Glacier retreated approximately 4 km from 2001 to 2011, reflected in our dh dt by the zone of rapid thinning near its front (approximately −25 m yr −1 ).Our repeat velocity measurements support the hypothesis that Marinelli is still in a retreat phase.We predict the retreat of Marinelli Glacier until it reaches a new topographic pinning point.Our dh dt map shows a similar pattern and degree of thinning at Darwin and CDI-08 glaciers, but we do not have the thickness or repeat speed measurements required to estimate flux and decipher whether thinning at these glaciers is dynamic.
All three Patagonian ice fields are losing mass.The NPI, with an average dh dt of −1.0 ± 0.1 m w.e.yr −1 (Willis et al., 2012b) is thinning at a slower rate than the CDI (−1.5 ± 0.6 m w.e.yr −1 ), which is thinning at approximately the same rate as the SPI (−1.6 ± 0.1 m w.e.yr −1 ) (Willis et al., 2012b).Other glaciers at the same latitude (e.g., the Kerguelen Islands, the island of South Georgia, and Heard Island) are thinning and retreating, and have undergone a similar degree of warming to Patagonia (Rosenblüth et al., 1995;Gordon et al., 2008;Thost and Truffer, 2008;Berthier et al., 2009).The CDI contributed 0.01 ± 0.004 mm yr −1 to sea level between 2000 and 2011.The combined NPI, SPI and CDI ice fields contributed 0.078 ± 0.008 mm yr −1 over the same period.This rate compares well with recent GRACE estimates (Jacob et al., 2012), and could be further refined by including, for example, DEMs generated from stereo highresolution imagery.Accurate basin-by-basin knowledge of accumulation rates would also help refine the cutoffs used to exclude spurious elevations, producing better dh dt .

Fig. 1 .
Fig. 1.Map of dh dt for the CDI (area indicated by red box in inset).Three glaciers, Marinelli, Darwin, and CDI-08, stand out with extensive thinning towards their fronts.The purple line indicates the divide between the "northern" and "southern" regions of the CDI.

Fig. 2 .
Fig. 2. Elevation values and dh dt for randomly selected pixels over Marinelli Glacier.The left-most elevation in each graph is the SRTM elevation at that pixel.Blue lines indicate the dh dt calculated for each pixel, elevation points bolded red are excluded from dh dt calculation.The bottom right panel shows the dh dt map, with numbered circles indicating the location corresponding to each graph.

Fig. 3 .
Fig. 3. Frontal variation history of Garibaldi Glacier from Landsat TM, ASTER and aerial photographs.The background is a 15/01/2011 ASTER image.

Fig. 5 .
Fig. 5. Scatter plot of E-W "bedrock" velocities for the 07/09/2001 to 25/09/2001 ASTER pair.The black line shows the linear trend fitted to the blue points, this is removed from the overall E-W velocity results.Red points are excluded when fitting the trend (they are greater than ±σ from the median value).This cutoff makes a negligible difference in the trend for this pair, but for other pairs it can change the trend significantly (e.g., > 1 m day −1 per 1000 m elevation).

Fig. 6 .
Fig. 6.Effect of elevation-dependent velocity correction for 07/09/2001 to 25/09/2001 ASTER image pair covering Marinelli Glacier.Glacier outlines in black, bedrock in gray (correction is not applied to water pixels).(a) Is a map of speeds with no correction applied, (b) is a map of speeds with correction applied.Correction reduces bedrock motion, the mean and standard deviation for the entire pair drops from 0.6 ± 0.4 m day −1 for the uncorrected speeds to 0.4 ± 0.3 m day −1 for the corrected speeds.

Fig. 7 .
Fig. 7. Date intervals for 22 pairs that produce usable pixel-tracking results.The most recent is a QB02 pair (blue), the next most recent is an ALOS pair (green), and the remainder are ASTER pairs (red).

Fig. 11 .
Fig. 11.dh dt and elevations for a longitudinal profile (starting towards the front) on Marinelli Glacier.A map of dh dt is shown in (a), the track used in (b) and (c) is plotted in green.(b) Gives elevation profiles for different dates along the green track in (a).The colorscale for the elevation profiles indicates the relative time of acquisition.Dark blue is the SRTM elevation (22/02/2000), the red track extending to the front is an ASTER DEM from 13/11/2007.From the elevation profiles it is clear that the front has retreated between 2-3 km between 2000 and 2007.(c) Shows the dh dt profile for the green track in (a).The colorscale for (a) is the same as (c).

Fig. 12 .
Fig. 12. Hypsometry for Garibaldi (a), Marinelli (b) and CDI-08 (c) glaciers from the SRTM DEM.The curve indicates what the AAR would be for the corresponding elevation.The black dot indicates the current ELA, the red dot indicates a 200 m upward shft.This would reduce the AAR for Garibaldi Glacier from 0.89 to 0.82, for Marinelli Glacier form 0.39 to 0.31, and for CDI-08 Glacier the AAR would drop from 0.71 to 0.46.
Fig. 13.dh dt and elevations for a longitudinal profile (starting towards the front) on Darwin Glacier.A map of dh dt is shown in (a), the track used in (b) and (c) is plotted in green.(b) Gives elevation profiles for different dates along the green track in (a).The colorscale for the elevation profiles indicates the relative time of acquisition.Dark blue is the SRTM elevation (22/02/2000), the red track extending to the front is an ASTER DEM from 15/01/2011.(c) Shows the dh dt profile for the green track in (a).The colorscale for (a) is the same as (c).

Table 1 .
Volume change rates and average dh dt for the 16 largest outlet glaciers on the CDI, with "northern", "southern" and whole icefield totals in the last three rows.Includes estimates of sub-aqueous volume change for Marinelli, CDI-08, Garibaldi and Darwin.
The ALOS pairs fail to capture the high speeds near the fronts of several glaciers due to decorrelation caused by strain and possibly melting.Raw ALOS SAR data are processed using ROI PAC and offsets are produced by "ampcor".The results are SNRfiltered and run through the elevation-dependent correction www.the-cryosphere.net/7/823/2013/The Cryosphere, 7, 823-839, 2013