Ablation from calving and surface melt at lake-terminating Bridge Glacier , British Columbia , 1984 – 2013

Bridge Glacier is a lake-calving glacier in the Coast Mountains of British Columbia and has retreated over 3.55 km since 1972. The majority of this retreat has occurred since 1991. This retreat is substantially greater than what has been inferred from regional climate indices, suggesting that it has been driven primarily by calving as the glacier retreated across an overdeepened basin. In order to better understand the primary drivers of ablation, surface melt (below the equilibrium line altitude, ELA) and calving were quantified during the 2013 melt season using a distributed energy balance model (DEBM) and time-lapse imagery. Calving, estimated using areal change, velocity measurements, and assuming flotation were responsible for 23 % of the glacier’s ablation below the ELA during the 2013 melt season and were limited by modest flow speeds and a small terminus cross-section. Calving and surface melt estimates from 1984 to 2013 suggest that calving was consistently a smaller contributor of ablation. Although calving was estimated to be responsible for up to 49 % of the glacier’s ablation for individual seasons, averaged over multiple summers it accounted between 10 and 25 %. Calving was enhanced primarily by buoyancy and water depths, and fluxes were greatest between 2005 and 2010 as the glacier retreated over the deepest part of Bridge Lake. The recent rapid rate of calving is part of a transient stage in the glacier’s retreat and is expected to diminish within 10 years as the terminus recedes into shallower water at the proximal end of the lake. These findings are in line with observations from other lake-calving glacier studies across the globe and suggest a common large-scale pattern in calving-induced retreat in lake-terminating alpine glaciers. Despite enhancing glacial retreat, calving remains a relatively small component of ablation and is expected to decrease in importance in the future. Hence, surface melt remains the primary driver of ablation at Bridge Glacier and thus projections of future retreat should be more closely tied to climate.


Introduction
Since the end of the Little Ice Age, glaciers across the globe have been shrinking at an accelerated rate (e.g.Dyurgerov and Meier, 2005;Radić and Hock, 2011;Gardner et al., 2013;Zemp et al., 2015).Although this retreat has been irregular, a general trend of 20th century retreat is pervasive and well correlated with an increase in global mean temperatures (Oerlemans, 2005).The reduction in ice cover in mountainous regions has raised concern about potential changes in the timing, volume, and duration of summer streamflow (e.g.Stahl et al., 2008;Marshall et al., 2011).These changes have major implications for hydroelectric projects, agriculture, aquatic habitat, water quality, and eustatic sea level rise (Barry, 2006;Radić and Hock, 2011;Gardner et al., 2013).While recent glacier retreat is well documented (e.g.Kaser et al., 2006), the projection of future retreat is critical to the management of water resources and understanding the evolution of riparian and aquatic habitats (Milner and Bailey, 1989;Cowie et al., 2014).
Due to their sensitivity to air temperatures and precipitation, variations in glacial size and volume serve as important high-altitude climate change indicators (Oerlemans, 2005;Kaser et al., 2006).However, glaciers that terminate in bodies of water have been shown to exhibit changes in mass balance that are at least partially independent of climate on decadal timescales (Warren and Kirkbride, 2003;Post et al., Published by Copernicus Publications on behalf of the European Geosciences Union.

2011
).This blurring of the climate-glacier signal is due to calving, which can be an important additional source of ablation (Benn et al., 2007a), and makes predictions of future retreat more difficult (Van der Veen, 2002;Motyka et al., 2002).However, the potential for calving glaciers to lose large volumes of ice over single seasons (even during years of positive mass balance) suggests that they can contribute disproportionately to eustatic sea level rise (Meier and Post, 1987;Dyurgerov and Meier, 2005), highlighting their important role in glacier response to climate.
Recently, there has been an increase in the number of studies examining the response of freshwater-calving glaciers to climate change.Most of the research exploring the dynamics of lake-calving glacier systems has focused on a few major regions: Alaskan glaciers Mendenhall and Yakutat (Motyka et al., 2002;Boyce et al., 2007;Trüssel et al., 2013), Tasman Glacier in the Southern Alps of New Zealand (Warren and Kirkbride, 2003;Dykes and Brook, 2010;Dykes et al., 2011), and several glaciers along the Patagonian Hielo Sur, most notably Perito Moreno, Nef, and Upsala glaciers (Warren et al., 2001;Stuefer et al., 2007;Sakakibara et al., 2013).Here we present new data from Bridge Glacier, a laketerminating outlet glacier of the Lillooet Ice field in the Coast Mountains of British Columbia, Canada.
The long-term retreat of calving glaciers has been found to follow a step-like pattern in which periods of stability are followed by a dramatic retreat, often coinciding with terminus flotation (Warren and Kirkbride, 2003;Boyce et al., 2007;Dykes et al., 2011).In many cases, flotation is achieved through thinning near the terminus due to successive years of high melt rates.Flotation can also be achieved by frontal retreat into deeper parts of a proglacial lake or fjord.At Mendenhall Glacier, climate-induced thinning led to increased instability and propensity to calve (Motyka et al., 2002) and eventually to the collapse of the terminus and retreat into shallower waters (Boyce et al., 2007).Similar findings have been made at Tasman Glacier in New Zealand (Warren and Kirkbride, 2003;Dykes and Brook, 2010) and in Patagonia (Warren and Sugden, 1993;Warren and Aniya, 1999;Skvarca et al., 2002), suggesting that retreat due to climatic warming may enhance calving rates over decadal timescales.Additionally, flotation can cause thinning due to an increase in terminus flow speeds (Rivera et al., 2012;Sakakibara et al., 2013), creating a positive feedback loop enhancing calving, and accelerating retreat rates.
This study investigates ablation due to calving and surface melt at lake-terminating Bridge Glacier.Here we define "ablation" as the process by which ice is lost from the glacier, both by calving and surface melt below the equilibrium line altitude (ELA) (Cogley et al., 2011), and do not include snow and firn losses.We use "surface melt" to refer to all net ablation of glacial ice through melting at the surface below the ELA; we assume ablation of snow is not significant, and it is not counted."Calving flux" is used throughout to refer to the ablation of glacial ice via frontal melting and ice-berg discharge at the terminus.Surface melt and the calving flux are estimated for the 2013 melt season from field measurements and a distributed energy balance model (DEBM).These results are then used to calibrate a mass balance model and calving model, which are applied to reconstruct calving fluxes and surface melt from 1984 to 2013.Calving rates and the relative contribution of calving to ablation from Bridge Glacier are then compared with findings from other laketerminating glaciers in Alaska, New Zealand, and Patagonia.Commonalities in the nature and timing of the calving flux and surface melt allow for a broad understanding of the pattern of calving losses over the transient calving phase of a retreating alpine lake-terminating glacier.

Study area
Bridge Glacier (50 • 48 11 N, 123 • 38 40 W), an outlet of the Lillooet Ice field, is located in the Pacific Ranges of the Coast Mountains of southwestern British Columbia, Canada, roughly 175 km north of Vancouver (see Fig. 1).The glacier had an area of 83 km 2 as of September 2013, extending from an elevation of over 2900 m at Bridge Peak, to 1390 m, where it terminates in a proglacial lake, locally known as Bridge Lake.Seventy-one percent of the glacier's area lies above 2100 m, which was approximately the average end-of-season snow line from 1985 to 2013.Bridge Glacier lies on the lee side of the humid coastal Pacific Ranges and terminates in a valley in the drier interior Chilcotin Ranges.Synoptic air flow is predominantly from the west, generating heavy snowfall on the highest elevation, most westerly areas, while the eastern flank of the glacier is drier, with a mean 1 May SWE of 600 mm (BC Ministry of Environment, 2014).
Bridge Lake has grown from under 2 km 2 in 1972 to over 6 km 2 in 2013 as the glacier retreated across an overdeepened basin (see Fig. 2).The distal (east) end of the lake traps numerous large (several hundred square metre surface area) tabular icebergs which are pressed along a submerged terminal moraine by persistent katabatic winds and have been present, in most cases, for several years.
Daily streamflow is measured by the Water Survey of Canada site "Bridge River (South Branch) Below Bridge Glacier" (Water Survey of Canada, 2015) and is available from 1978 to present.The hydrometric site is located less than 2 km downstream of the distal (east) end of Bridge Lake, and 60 % of its catchment area (144 km 2 ) is occupied by Bridge Glacier.Temperature and precipitation for the region are obtained from Environment Canada climate station Vancouver International Airport, BC (49 • 12 N, 123 • 11 W; elevation is 4 m; ID no.1108447) (Environment Canada, 2015).Air temperature at the Vancouver climate station is a significant predictor of both mean annual flow at the Bridge River gauge (r 2 = 0.65, p < 0.001) and of Bridge Glacier ELAs (r 2 = 0.32, p = 0.001), suggesting it is an adequate broadscale climatic proxy.

Weather data
Three automatic weather stations (AWSs) collected data from 20 June to 12 September 2013 to provide input data for a DEBM (see Fig. 1).One weather station was installed on the glacier (glacier AWS), collected air temperature, humidity, wind speed, and direction, and reflected shortwave radiation at 10 min intervals.A second weather station (ridge AWS), installed on a ridge ∼ 250 m above the glacier toe and hence shielded from strong, persistent katabatic flow, collected ambient temperature and solar radiation.A third weather station, located along the shore of Bridge Lake (lake AWS) approximately 3 km from the terminus, on a partially submerged end moraine, measured incoming longwave radiation, air temperature, humidity, wind speed, and rainfall.Rainfall was also measured at an exposed nunatak north of the main arm of the glacier (Nunatak TLC on Fig. 1) to estimate the precipitation gradient over the glacier tongue.
Incoming shortwave and longwave radiation was collected off-glacier due to our inability to ensure the sensor remained level at glacier AWS.
In order to evaluate surface melt derived from melt modelling, 3 m long ablation stakes were installed at six locations in the ablation area between 1500 and 1600 m.Due to logistical challenges, and to obtain results that could also be used to evaluate velocity estimates, the stakes were located within 2 km of the terminus (Fig. 1).The stakes were installed on 18 June and were resurveyed and re-drilled on 19 July and 13 September 2013.

Bathymetry
Bathymetric data were collected using a Lowrance HDS Gen2 depth sounder (Lowrance, 2011), with a depth range of 500 m and horizontal GPS accuracy of ±5 m.Due to the presence of large, unstable icebergs throughout the lake, depth measurements were taken at 893 discrete points in an irregular grid.Access to the terminus and the middle part of the lake was hindered by the presence of icebergs, necessitating the inclusion of an additional 74 points which were added by linear interpolation using known depths along east-west transects.The bathymetric data were processed using the gstat package in R (R Core Team, 2013;Pebesma, 2004), and interpolated onto a 10 m grid using inverse distance weighting.Water depth for the 2013 calving front was estimated from a cross-section parallel to, and roughly 500 m from, the June 2013 terminus position.

Flow speed
The terminus flow velocity was measured by tracking features on images taken by two time-lapse cameras (TLC), at Nunatak TLC and lake TLC, set up to capture the floating terminus and the glacier surface roughly 1 km up-glacier.
Points were tracked manually using Tracker video analysis and modelling tool (Brown, 2014).Raw pixel displacement was converted into distances using known camera angles and several ground control points following Harrison et al. (1992) and Eiken and Sund (2012) (see Chapter 4 in Chernos (2014) for further details).Eight points in close proximity on the glacier surface (< 200 m) were tracked from each camera throughout the study period using daily noontime images.
Filtering routines discarded roughly 10 % of the tracked data points due to a negative measured displacement or loss of target.Daily surface velocities were generated by averaging the daily displacements for each tracked point, and the average summer velocity was calculated by averaging the total displacement for each tracked point throughout the study period.Study-period time-lapse velocity measurements were complemented with an end-of-summer survey of ablation stakes; results were found to agree within the error of our Garmin eTrex GPS (±5 m).

Satellite imagery and elevation data
The change in terminus area during the 2013 study period was computed from Landsat images on 23 June and 11 September 2013.Shapefiles for both scenes were generated by manually delineating the terminus in Google Earth.The change in area was then calculated using the rgeos package in R (R Core Team, 2013).Annual terminus positions and ELAs from 1984 to 2012 were reconstructed from Landsat imagery.All Landsat images were taken between 12 September and 24 October to represent end-of-season snow lines.Annual terminus retreat rates (m a −1 ) were calculated by measuring the areal retreat, averaging it by the terminus cross-section (width), and correcting for a full calendar year.
Because the intent was to model only ice melt, and not the melting of snow cover, the DEBM was constrained to the area below the snow line.Daily snow line elevations were determined by loess smoothing of snow line elevations estimated from nine Landsat images obtained from the Landsat-Look Viewer (US Geological Survey, 2014)  We applied a distributed energy balance model driven by data from the three AWSs and a digital elevation model of the glacier surface from 2006.As our purpose was to calculate surface melt below the ELA, we only consider ice melt (not snow or firn melt), and hence we only modelled surface melt for the area of exposed glacial ice below the snow line at each time step.
Surface melt (M), in m (we) day −1 , is calculated as where Q M is the sum of available energy at the surface (W m −2 ), L f is the latent heat of fusion (3.34 × 10 6 J kg −1 ), and ρ i is the density of ice (917 kg m −3 ).Energy supplied to the glacier surface is positive, while energy flux away from the surface is negative.The available energy for melt is calculated as where Q * is the net radiation, Q H and Q E are the sensible and latent heat flux, and Q R is sensible heat of rain.All energy fluxes are in W m −2 .We assume that all energy fluxes occur at the ice surface (Oerlemans, 2010;Munro, 2001); subsurface and subglacial melt is neglected.

Net radiation
Net radiation (Q * ) is calculated as the sum of incoming (↓) and outgoing (↑) shortwave (K) and longwave (L) radiation as follows: where S and D are the direct and diffuse components respectively of incident shortwave radiation, and α is the albedo of ice.
Reflected shortwave radiation was measured on-glacier over bare ice in the ablation area, throughout the melt season.Incoming shortwave radiation was measured from the off-glacier ridge AWS.Differences in shading between the two sites were found to be negligible.To minimise the effects of small discrepancies in shading, uneven cloud patterns, and low solar angle errors (Oerlemans, 2010), the daily ice albedo (α) is assumed constant throughout the day, and is calculated as where the integrals are calculated over the period of daylight each day.Albedo was only estimated from glacier AWS and was kept constant across the glacier.Although this limits the model's representativeness over the whole glacier, given the model is only applied over exposed glacial ice, this simplification is not expected to have an appreciable impact on the volume of melt modelled.Direct shortwave radiation (W m −2 ) for each grid point on the glacier surface is calculated as where K ex i,j is the potential direct solar radiation at grid point (i, j ) and K ex is the potential direct solar radiation at glacier AWS.Measured global radiation was separated into direct and diffuse components based on the ratio of observed to potential shortwave radiation following Collares-Pereira and Rabl (1979) and Hock and Holmgren (2005).Potential direct radiation was corrected for slope geometry and diffuse shortwave radiation is calculated for all cells when K ex > 0 (Hock and Holmgren, 2005;MacDougall and Flowers, 2011) as where D o is the global diffuse radiation, and φ i,j is the sky view factor at each grid point (i, j ).Due to the difficult logistics (and likely spatially variable results) involved in measuring the albedo for the surrounding non-glaciated terrain (α terrain ), a constant value of 0.17 was assumed, which is typical of dark, rocky surfaces (Oke, 2000).Uncertainties associated with this assumption should be minor in practice, given that sky view factors for the glacier are high (∼ 0.95).The sky view factor was calculated using SAGA GIS software following Oke (2000) and a 25 m lidar DEM.The algorithm integrates the maximum horizon angles (H ) for each grid cell, for each azimuth angle (1 • interval).A maximum 10 × 10 km search window was implemented to reduce computation time.
In order to spatially distribute incoming shortwave radiation, each grid point is modelled as either shaded or sunlit.A shading algorithm was implemented that calculates the maximum horizon angle for each grid point within a 10 × 10 km window, using 10 • azimuth bins.At each time step, if the horizon angle is greater than the elevation angle (Z), the grid point is shaded and only receives diffuse radiation.For times when the horizon angle is smaller than elevation angle, the grid point receives both direct and diffuse radiation.
www.the-cryosphere.net/10/87/2016/The Cryosphere, 10, 87-102, 2016 Incoming longwave radiation was measured directly at the lake AWS and was computed at each grid point as follows: where L terrain is the longwave radiation emitted by surrounding terrain.Longwave radiation emitted by the terrain was computed using the Stefan-Boltzmann law with a terrain emissivity of 0.95 (Oke, 2000) and the assumption that terrain temperature is equal to air temperature.Although atmospheric longwave radiation over the glacier and at an offglacier site could be expected to differ due to the effects of katabatic flow on near-surface air temperature and humidity, the difference in humidity between glacier and lake AWSs was less than 10 %, while air temperatures at lake AWS are 1.6 • C warmer.Furthermore, Shea (2010) measured incident longwave radiation at on-glacier and off-glacier sites at the same elevation at Place Glacier and found little systematic difference over all sky conditions.Longwave radiation emitted by the ice surface was computed from the Stefan-Boltzmann law using an emissivity of 0.98 (Oke, 2000).The surface temperature was set to 273.15 K.This assumption of a continuously melting ice surface is reasonable considering that on-glacier air temperature was always above 0 • C during the study period and only below 2 • C for 3 h.

Turbulent heat fluxes
Sensible and latent heat fluxes are calculated using the bulk transfer approach: where c air is the specific heat capacity of air (1006 J kg −1 K −1 ), u is the wind speed (m s −1 ), T g is the on-glacier air temperature, T s is the glacier surface temperature (held constant at 273.15 K), L v is the latent heat of vaporisation (2.50 × 10 6 J kg −1 ), e g and e s are the vapour pressures (hPa) of air and glacier surface (held constant at 6.11 hPa, assuming the glacier surface is at the melting point), and P is the atmospheric pressure (hPa) at glacier AWS.The turbulent transfer coefficient C (unitless) is calculated using stability corrections based on the bulk Richardson number, using a roughness length for momentum of 2.5 mm for ice (Munro, 1989(Munro, , 2001;;Pellicciotti et al., 2005) and calculating the roughness length for temperature and vapour following Hock (1998).
Air temperature was distributed over the glacier surface using the approach developed by Shea and Moore (2010), which accounts for the effects of katabatic flow.In this approach, the magnitude of katabatic forcing was modelled as a function of the temperature difference ( T ) between the onglacier glacier AWS (T g ) and off-glacier ridge AWS (T a , outside the katabatic boundary layer).Temperature differences were separated into upslope (northeasterly) and downslope katabatic (southwesterly) flows, based on the wind directions of glacier AWS.Linear regression against off-glacier temperature (T a , Fig. 3) showed a positive linear increase in T , indicating the magnitude of katabatic forcing increased with increasing off-glacier air temperatures.Conversely, T did not significantly vary as a function of off-glacier temperatures during upslope flow, although temperatures above 10 • C during these episodes were rare.The elevations of both weather stations were within 100 m, and small corrections to potential temperature using a −6 • C km −1 lapse rate (as used in Stahl et al., 2008 andShea, 2010) did not produce a meaningful difference in the linear fit.
On-glacier air temperature for each grid point is modelled as a function of the katabatic temperature depression where and T * is the threshold temperature difference at which katabatic flow is observed.The magnitude of katabatic forcing for each point on the glacier, k 1 , was calculated using statistical coefficients and glacier flow path lengths (Shea, 2010;Chernos, 2014).Flow path lengths for the glacier were calculated using the Terrain Analysis -Hydrology module of SAGA GIS (Quinn et al., 1991;SAGA Development Team, 2008).During periods when wind direction was upslope, temperatures were distributed using the on-glacier temperature, T g , and a temperature lapse rate of −6 • C km −1 (Stahl et al., 2008).

93
Wind speed across the glacier was distributed as a function of katabatic forcing and ambient temperatures, following Shea (2010).When the measured wind direction at glacier AWS was downslope, wind speed at glacier AWS showed a positive linear correlation with (off-glacier) ridge AWS air temperature, while upslope wind speeds showed no discernible trend.Therefore, when the wind direction at glacier AWS is upslope, wind speed was held constant across the glacier within our melt model, using measured wind speeds from glacier AWS.
Vapour pressure is calculated from measured relative humidity and saturation vapour pressure (e sat ) which was calculated using Teten's formula (Murray, 1967).Relative humidity, measured at glacier AWS, was held spatially constant across the glacier for each time step, and saturation vapour pressure was calculated from distributed on-glacier air temperatures.

Melt contribution from rain
Energy supplied to the surface due to rain was calculated following Hock (2005): where R is the rainfall rate (m s −1 ), measured at the lake AWS (and missing values were filled with measured data from Nunatak TLC), and ρ w and c w are the density (1000 kg m −3 ) and specific heat of water (4180 J kg −1 K −1 ).The temperature of rain, T R , is assumed equal to the ambient off-glacier air temperature and was corrected for elevation using a −6 • C km −1 lapse rate.Since we observed no significant elevational or east-west precipitation gradient between Nunatak AWS and lake AWS, rainfall was held constant across the glacier.

Modelling calving flux
Calving fluxes are calculated from measured retreat rates and flow speeds, as well as estimates of ice thickness derived from bathymetry.The volume of ice discharged through calving from the glacier terminus, Q calving (m 3 a −1 ), i.e. the calving flux, is quantified as where dA T dt is the change in glacier surface area at the terminus (m 2 a −1 ), U is the terminus flow velocity (m a −1 ), and H I and W are the ice thickness (m) and glacier width (m) at the terminus.Subaqueous melt at the ice front is assumed to be negligible with respect to the magnitude of the calving flux.
The thickness of ice at the terminus was approximated by assuming that the terminus was right at the threshold for flotation.Using the height above buoyancy criterion ( Van der Veen, 1996;Benn et al., 2007b), the ice thickness (H I ) can be calculated as where H b is the height of ice above the waterline (m), D W is the water depth, and ρ w and ρ i are the densities of water and ice.During the melt season, large tabular icebergs calved and showed limited mobility, suggesting that the glacier was at or near the boundary criterion for flotation.There was a notable inflection point (Fig. 4) roughly 500 m from the endof-season terminus, where the surface slope becomes flat or slightly reclined and had remained stationary since 2012, and where we assume that the terminus transitions from grounded to floating.
The calving flux between 1984 and 2012 was computed from historical terminus positions, average retreat area, water depth taken from lake bathymetry (Fig. 5), estimated ice thickness, and measured velocity from the 2013 field season.

Historical ablation
Estimates of historical annual surface melt were derived using ELA observations and a fitted piece-wise linear mass balance gradient derived using mass balance observations from several glaciers in the region, including Bridge Glacier (Shea et al., 2013).Below the snow line, the net balance (b n ) at a point is equal to the surface melt of exposed glacier ice and was estimated using glacier hypsometry from the 2006 lidar DEM, where and is calculated for the elevation of every point, z (m a.s.l.), below the ELA.
Results from the distributed energy balance model provide a means to evaluate the mass balance gradient for Bridge Glacier that can be used to estimate surface melt below the ELA for previous years.The coefficient value (b 1 = 6.62 mm (we) m −1 ) taken from Shea et al. (2013) gives a lower estimate of surface melt for the 2013 melt season relative to that calculated with the DEBM.The mass balance gradient generated by the DEBM suggests a value of b 1 = 9.07 mm (we) m −1 (Fig. 6); this value was used for all years.The glacier area was determined from the end-ofseason calving margin.All glacial surface areas that calved prior to 2013 are estimated in Eq. ( 14) by assuming an elevation of 1400 m (a.s.l.).

Climate and retreat
The annual retreat of Bridge Glacier was composed of several stages.Retreat was slow prior to 1991, characterized by www.the-cryosphere.net/10/87/2016/The Cryosphere, 10, 87-102, 2016  small calving events along the shallow proglacial lake margin.The average rate of retreat between 1972 and 1991 was 21 ma −1 but accelerated to 144 ma −1 after 1991, punctuated by high annual retreat rates followed by years of relative terminus stability and the appearance of large tabular icebergs in the lake.The rate of retreat accelerated again after 2004 to ∼ 400 m a −1 (Fig. 7e).Since 1991, the glacier has retreated over 3.55 km, with occasional years of rapid retreat associated with calving of large, tabular icebergs, indicative of a floating terminus.The substantial retreat that Bridge Glacier has undergone since 1991 does not fully follow glacial melt predictors such as summer air temperature, winter precipitation, mean an-nual streamflow, or equilibrium line altitudes (Fig. 7).For example, air temperature anomalies became dominantly positive in the 1980s without a corresponding change in the retreat rate.Additionally, from 1988 to 1998 summer temperatures, equilibrium line altitudes, and mean annual flows from Bridge River were all above the 30-year average (Fig. 7b-d), suggesting above average melt.However, this period of elevated melt conditions did not continue into the 21st century as retreat continued to accelerate.

The 2013 surface melt
From 20 June to 12 September 2013, our model predicted surface melt ranging from 5.9 m we near the terminus to 0 at The Cryosphere, 10, 87-102, 2016 the ELA, yielding a total ablation volume below the ELA of 0.124 km 3 (Fig. 8).Melt rates were greatest along the main tongue of the glacier due to high sensible heat flux driven by persistent katabatic flow.The southernmost tributary glacier showed relatively low melt rates relative to similar elevations on the main tongue, most likely due to the fact that it remained sheltered from high winds and its north-facing aspect allowed for substantial shading throughout the melt season.
The snow line was at the terminus until 15 June, and the ablation area had become snow-covered again before 20 September, suggesting our field instrumentation captured all but 12-15 days of melt in the 2013 season.We estimate that surface melt during this period was less than 10 % of the total surface melt during the study period.
Modelled surface melt agreed within ±0.2 m we for four of the five ablation stakes (Fig. 9), representing an error of less than 5 % of the measured value.We estimate the uncertainty in our ablation stake measurements were ±0.02 m for each survey (three in total), corresponding to an estimated measurement uncertainty of ±0.06 m.Additionally, we estimate that uneven glacial melt due to heterogeneity in surface debris cover, meltwater pooling, and uneven terrain was on the order of ±0.15 m, based on observations of the glacier surface.Therefore, we estimate a total uncertainty of ±0.21 m for the melt measurements.Measured melt at ablation stake D, located roughly 400 m up-glacier (∼ 100 m increase in elevation) from glacier AWS and stake A, was up to 0.8 m less than other nearby stakes (including stake E, which is 100 m higher in elevation and further up-glacier), suggesting that there may have been localised effects shielding the stake or otherwise inhibiting melt at this site relative to the higher melt rates observed elsewhere in the ablation area, such as a locally elevated albedo.In order to quantify the potential impact of our various modelling assumptions on modelled surface melt below the ELA, we re-ran the DEBM under multiple scenarios (see Table 1).Given that albedo was only measured at glacier AWS (average for the season was 0.21), the model was also run with a seasonally and spatially constant value of 0.4 (Brock et al., 2000), considered a high envelope for the season given no snowfall was observed.In order to further test the validity of our basin-wide snow line estimates, the model was re-run by delaying snow line retreat by 1 week.In order to test the sensitivity of our sensible heat flux calculations using flow path lengths, two additional model runs were performed: one with a constant temperature lapse rate (−6 • C km −1 ) and another with a spatially constant wind speed (taken from glacier AWS).In order to test the sensitivity of our model to projected thinning of the glacier, the DEBM was run with glacier elevations artificially depressed by 50 m.Finally, as a benchmark for model sensitivity, the DEBM was also run with air temperatures increased by 1 • C.
In these model sensitivity simulations, the largest changes in modelled surface melt are due to a constant albedo of 0.4 and by delaying snow line retreat by a week.At Bridge Glacier, net radiation is responsible for 60-70 % of the melt energy supplied during the melt season, making surface melt most sensitive to changes to the radiative energy balance.Conversely, surface melt shows relatively little sensitivity to changes in glacier elevation and air temperature increases, suggesting that ignoring glacier thinning in the DEBM does not materially impact total volumes of surface melt.

The 2013 calving flux
Over the 85-day study period in 2013, a change in terminus area (dA T ) of −0.297 km 2 was measured from repeat terminus delineations, corresponding to a terminus retreat of 281 m a −1 .The average velocity at the terminus (U ) was 139 m a −1 (see Fig. 5), across a width (W ) of 1055 m, yielding a cross-sectional area of 0.0342 km 2 .The median water depth was 91 m, corresponding to a height above buoyancy of 9.9 m and an estimated ice thickness of 109 m.Combining these measurements in Eq. ( 12) yields an estimated calving flux of 0.0362 km 3 for the study period.
Adding the volume of ablation due to calving with surface melt during the same period yields a total volume of 0.160 km 3 (Fig. 10).For the 2013 melt season, calving accounts for 23 % of the total ablation at Bridge Glacier, equivalent to an additional 1.3 m of surface melt over the entire ablation area.
A 60 m uncertainty in measuring the terminus crosssection (W ) (equal to 2 Landsat pixels) is applied.The uncertainty of dA dt is estimated as 7200 m 2 a −1 (2 × 60 m × 60 m).Bathymetric error is calculated at 5.6% and was found by differencing two bathymetric models produced using a randomly selected half of the collected water depth pointmeasurements.The ice thickness uncertainty is estimated as 5.6 % plus an additional 10 m to account for changes in sedimentation and ice thickness relative to water depth.Before 1991, the terminus was not floating; therefore, an ice thickness uncertainty of 60 m is estimated to account for a range of grounded terminus geometries.Between 1991 and 2004, bathymetry has poor data coverage, and a ice thickness uncertainty of 33 m is estimated.Historical terminus velocities were assumed to be approximately equal to the average 2013 summer flow speed (140 m a −1 ), and annual calving rates are calculated with 70 m a −1 (50 %) potential variability around the 2013 mean.

Historical ablation
Between 1984 and 2013, surface melt showed a minor decrease over time, which can be attributed to the loss of surface area in the lowest reaches of the glacier due to calving and retreat (Fig. 11).Surface melt in 2013 was above the 30-year average but within 1 standard deviation of the mean (x = 0.107 km 3 a −1 , s = 0.018 km 3 a −1 ).The ELAs varied between 1926 and 2202 m during the period; in most years it was between 2050 and 2150 m.
Historical calving losses are characterized by several years of high flux and periods of relative stability.The magnitude of the calving losses increased once the glacier achieved flotation in 1991 and were minimal before then.From 1992 to 1994, the calving flux increased to 0.020-0.029km 3 a −1 (19-27 % of the total annual ablation), before a 2-year period of low flux (< 0.015 km 3 a −1 ).From 1997 to 2000, the calving flux increased again (0.023-0.052 km 3 a −1 ), while Uncertainties in calculations of surface melt below the ELA are estimated assuming a ±75 m uncertainty in ELA elevation due to timing of available Landsat images to measure the snow line or 22 % found by Shea et al. (2013), whichever is greater.The ELA uncertainty estimate is to account for errors that cannot be adequately quantified without additional historical data.For example, it is difficult to confirm the linearity or interannual consistency of the net balance gradient without several seasons of mass balance measurements (as was done by Shea et al., 2013), which changes annually depending on summer weather and winter snowpack depth and distribution.For the 2013 study period, the shape of the DEBM-derived mass balance gradient mirrors the seasonal snow line retreat rate derived from the Landsat images, where early in the season the snow line retreated quickly, then rose less than 50 m from August onwards.
Glacier hypsometry is not adjusted during the 1984-2013 study period, and it is based on a 2006 lidar survey.Although thinning invariably affects the elevation, and therefore air temperatures predicted from our lapse rate, the elevation difference between 1970 and current terminus position is estimated at less than 200 m.Moreover, sensitivity analyses for the DEBM (Table 1) show that a 50 m change in glacier surface elevation had only a 0.7 % increase in the volume of ablation below the ELA, while even a relatively large increase of 1 • C resulted in a less than 5 % increase in surface melt below the ELA.Therefore, even large changes, such as an albedo roughly double our measured value, still report surface melt differences less than the 22 % reported by Shea et al. (2013) and are used as our uncertainty bounds in our surface melt model.

Controls on calving
During the 2013 melt season, calving was a moderate contributor to ablation relative to surface melt below the ELA at Bridge Glacier.Calving losses in this system are controlled by glaciological and topographical controls that ultimately limit the magnitude of the calving flux.The glacier width at the flux gate was just over 1 km, which restricts the volume of ice that can reach the floating terminus, in turn limiting the size of calving events.In contrast, the ablation area in 2013 was 27.6 km 2 , allowing for surface melt processes to act over a much larger area and contribute a substantially larger volume of surface melt than possible from the calving front.
Relatively modest glacier flow speeds at the terminus also limit the volume of ice delivered to the terminus and calving.Flow velocity at Bridge Glacier is moderate due to gentle gradients in the lower reaches of the glacier, as well as a relatively narrow cross-sectional area.A gentle surface slope reduces the gravitational stresses, while narrow valley sidewalls constrict glacier flow by providing substantial lateral drag (Benn et al., 2007a;Koppes et al., 2011), both of which limit flow speeds.Near-terminus flow speeds at Bridge Glacier are 1 to 2 orders of magnitude smaller than those observed at larger tidewater calving glaciers in Patagonia and Alaska (Rivera et al., 2012;Koppes et al., 2011;Meier and Post, 1987;Motyka et al., 2003) and reflect a smaller mass turnover, similar to lake-terminating glaciers Mendenhall and Tasman (Boyce et al., 2007;Dykes et al., 2011).
The bathymetry of Bridge Lake also plays an important role, where interannual calving fluxes mirror average and maximum water depths at the terminus.This relationship suggests that water depths are a large-scale control on calving in lacustrine environments.In particular the onset of terminus flotation remains the largest variable responsible for initiating rapid calving losses and retreat, a finding that mirrors results elsewhere (Boyce et al., 2007;Dykes and Brook, 2010;Trüssel et al., 2013;Sakakibara et al., 2013).However, this relationship does not necessarily suggest that water depth can drive annual (or sub-annual) calving rates.While floating temperate ice tongues have been shown to be unstable (Van der Veen, 1996;Benn et al., 2007a), often leading to disintegration and dramatic retreat, several examples exist of floating termini remaining intact for multiple years.For example, at Mendenhall Glacier an unstable floating terminus remained intact for approximately 2 years (Boyce et al., 2007), while Yakutat Glacier sustained a floating ∼ 3 km terminus for over a decade (Trüssel et al., 2013).Similar results from Bridge Glacier, where the floating terminus had multiple seasons of negligible calving (2001,2002,2007), suggest that water depth offers insufficient predictive power for annual calving fluxes.

The relative importance of calving
From 1984 to 2013, the calving flux increased from an almost negligible annual value to a flux responsible for between 20 and 45 % of the glacier's annual ablation.The recent increase in calving flux closely follows water depth at the terminus, where the largest calving fluxes coincide with the terminus retreating into the deepest parts of Bridge Lake in 2003-2011.While this relationship suggests that buoyancy is a primary driver of multiannual calving at Bridge Glacier, it also implies that the high rate of calving currently observed is unsustainable over the coming decades and is instead part of a transient phase as the glacier continues to retreat up-valley and into shallower waters.
Although calving contributed less than one quarter of the total ablation from Bridge Glacier during the 2013 melt season, during 3 of the last 10 years the calving flux was on par with the volume of ablation due to surface melt below the ELA.However, large annual calving fluxes do not persist over several consecutive seasons and are instead followed by several years of minor calving losses, even though the ter-minus remained in the deepest part of the lake.The pattern of a high-magnitude calving year followed by several lowflux years is consistent with the notion that glacier dynamics respond to large calving events by alleviating terminus instability and inhibiting future calving (Venteris, 1999;Benn et al., 2007b).Following a large calving event, the glacier geometry changes, and buoyant forces can be redistributed or relieved, promoting terminus stability.
Historical reconstructions of calving and surface melt suggest that climate is the largest variable affecting long-term ablation rates at Bridge Glacier.Although calving has produced substantial ablation during the last 10 years, calving fluxes for most studied lacustrine glaciers have been shown to strongly correlate with the terminus remaining in deep water (Warren and Aniya, 1999;Van der Veen, 2002;Benn et al., 2007b).Given that Bridge Glacier is approximately 850 m from the proximal end of Bridge Lake (Fig. 4), and that the average calving rate over the last 5 years is 299 m a −1 , it is probable that calving will only remain a substantial component of ablation for another decade, suggesting that the current rate of calving is transient and unsustainable.Given that surface melt below the ELA is the primary contributor of ablation at Bridge Glacier, the glacier's future mass balance is more dependent on climatic conditions.

Bridge glacier and other lake-calving systems
Observations of the magnitude and frequency of calving at Bridge Glacier fall in the middle of a continuum of studied lake-terminating glaciers worldwide (see Table 2).The calving rate for Bridge Glacier (281 m a −1 in 2013) is larger than that for smaller glaciers in New Zealand, such as Maug, Grey, and Hooker (Warren and Kirkbride, 2003), and for Mendenhall Glacier in Alaska (Motyka et al., 2002;Boyce et al., 2007).Conversely, calving rates at the larger Patago- nian glaciers Leon, Ameghino, and Upsala are up to an order of magnitude greater than what we found at Bridge (Warren and Aniya, 1999;Sakakibara et al., 2013).Bridge Glacier's calving rate is controlled by moderate water depths and flow speeds.Higher calving rates are associated with greater water depths and significantly larger terminus velocities.Large Patagonian and Icelandic glaciers have terminus velocities of up to 1810 m a −1 (Haresign, 2004), an order of magnitude greater than what we measured at Bridge Glacier (140 m a −1 in 2013).Conversely, smaller calving glaciers in New Zealand terminate in shallow lakes (< 50 m) and many have low flow speeds (< 70 m a −1 ).Bridge Glacier's calving rate in 2013 (281 m a −1 ) also agrees well with first-order linear models relating calving to water depth (D W ) (Funk and Röthlisberger, 1989).Using the revised relationship from Warren and Kirkbride (2003), where U c = 17.4 + 2.3D W , the modelled calving rate (U C ) for Bridge Glacier is calculated as 268 m a −1 , which is within 13 m a −1 of the rate we observed in 2013.
Lake temperatures also appear to play a role in controlling the calving rate.Many Patagonian ice fields terminate in large lakes where water temperatures are up to 7.6 • C (War-ren and Aniya, 1999), significantly warmer than the wellmixed 1 • C water observed at Bridge Lake (Bird, 2014).This difference is most likely related to the surface area to depth ratio of the proglacial lakes.Bridge Lake, at 6.3 km 2 , is small relative to the much larger lakes of Southern Patagonia, while only marginally shallower.As such, many large Patagonian proglacial lakes contain vast areas that are free of the strong cooling influence of glacier runoff and trapped icebergs and can warm significantly, promoting thermal undercutting and enhancing further calving (Rohl, 2006;Rignot et al., 2010;Robertson et al., 2012).
Bridge Glacier shares similar calving characteristics with both Tasman and Mendenhall Glaciers, both of which have undergone significant retreat as they transitioned from grounded to floating termini (Boyce et al., 2007;Dykes et al., 2011).During this transition, terminus velocities increased at Tasman from 69 to 218 m a −1 (Dykes and Brook, 2010;Dykes et al., 2011), while the calving rates for both glaciers increased from 50 m a −1 to between 227 and 431 m a −1 (Boyce et al., 2007;Dykes et al., 2011); these rates are consistent with what we found at Bridge Glacier.For both Tasman and Mendenhall Glaciers, water depth and buoyancy also control the magnitude of calving (Boyce et al., 2007;Dykes et al., 2011;Dykes, 2013), suggesting that the majority of the ice discharged from the terminus is triggered by buoyant forces.
The relative contributions of calving and surface melt to ablation below the ELA at Bridge Glacier are comparable to other studies worldwide.While calving at Bridge Glacier is responsible for an average of 10-25 % of total ablation, Yakutat Glacier experienced calving fluxes between 7.9 and 16.8 % of total mass loss from 2000 to 2007 and 2007 to 2010 (Trüssel et al., 2013).These percentages are much higher than what has been observed at Mendenhall Glacier, where calving was responsible for 2.6-4.0 % of the long-term volume change (Boyce et al., 2007).
The differences in the relative contributions of calving to ablation points to different stages in a relatively uniform pattern of retreat present in lake-calving glaciers.Studies from Patagonia (Sakakibara et al., 2013), Alaska (Boyce et al., 2007;Trüssel et al., 2013Trüssel et al., , 2015)), and New Zealand (Dykes et al., 2011;Dykes, 2013) all report glacier thinning, followed by terminus flotation and a rapid step-like retreat, something that is echoed at Bridge Glacier.These findings hint at a common large-scale behaviour of retreating laketerminating glaciers and suggest a broad applicability in the region and across the globe of a pattern of transient high calving contributions to ablation as the glacier retreats across an overdeepened lake.

Conclusions
Bridge Glacier is a lake-terminating glacier in the Coast Mountains of British Columbia that has retreated over 3.55 km since 1972, with the majority of retreat occurring after 1991.This retreat was independent of regional warming trends and was enhanced by significant calving losses as the glacier terminus retreated into deeper waters.While calving has accelerated Bridge Glacier's retreat, estimates of surface melt and calving for the 2013 melt season indicate that calving was only responsible for 23 % of the glacier's total ablation.The contribution of calving to ablation was limited by modest terminus flow speeds, relatively narrow side walls in the lower glacial tongue, and lake depth at the terminus.
Estimates of calving and surface melt from 1984 to 2013 suggest that calving did not significantly contribute to surface melt before 1991.From 1991 to 2003 calving rates increased, and the calving flux was on par with ablation from surface melt below the ELA in 2005, 2008, and 2010.Although individual years had large calving fluxes, multi-year averages between 1991 and 2013 show that the calving flux only accounted for between 10 and 25 % of the glacier's annual ablation below the ELA.The rapid calving rates observed from 2009 to 2013 at Bridge Glacier are part of a transient stage in retreat as the glacier terminus passed through an overdeepened, lake-filled basin and are not expected to remain a consistently large source of ablation in the coming decades.These findings are in line with observations from other lake-calving glacier studies across the globe and suggest a common large-scale pattern in calving-induced retreat in lake-terminating alpine glaciers.Despite enhancing glacial retreat, calving remains a relatively small component of ablation and is expected to decrease in importance in the future.Hence, surface melt remains the primary driver of ablation at Bridge Glacier and, as such, projections of future retreat should be more closely tied to climate.

Figure 1 .
Figure 1.Bridge Glacier study area, instrumentation, and select terminus positions from 1973 to 2013.The DEM is from winter 2006 and contour intervals are 100 m.The 2013 end-of-season snow line was 2103 m.Insert shows the location of Bridge Glacier within southwest British Columbia.

Figure 2 .
Figure 2. Landsat imagery from 1985 to 2012, showing retreat of Bridge Glacier and opening of Bridge Lake.All images have the same orientation and scale as the upper left panel.
between 1 June and 19 September 2013.Multiple measurements of snow line altitude across the glacier surface were taken for each image and averaged to produce a basin-wide snow line elevation.Elevation data for the glacier surface were obtained using a 25 m resolution LIDAR digital elevation model (DEM) from 2006 (from C-CLEAR by M. Demuth, C. Hopkinson, and B. Menounos; see Acknowledgements).The DEM was resampled to 50 m to reduce computation time and remove unrealistic elevation changes produced at the junction of two map tiles.

Figure 3 .
Figure 3. On glacier temperature depression ( T = T a − T g ) as a function of ambient air temperatures (T a ) from ridge AWS (outside the katabatic boundary layer).The blue line is the significant fit (p < 0.01) for downslope/katabatic winds and the red line is the non-significant fit for upslope winds, while the dashed grey line demarcates no temperature depression.

Figure 4 .
Figure 4. Photograph of Bridge Glacier terminus in September 2013, showing the approximate location of the inflection point and of the proximal edge of Bridge Lake.

Figure 5 .
Figure 5. Map showing 2013 Bridge Lake bathymetry, 2013 study period flow vectors (velocity in m yr −1 at scale of the map), ablation/velocity stakes (black dots), flux gate, and historical terminus positions.

Figure 6 .
Figure 6.Modelled mass balance gradient from Shea et al. (2013) and a tuned coefficient using distributed energy balance modelling from the 2013 melt season.

Figure 8 .
Figure 8. Modelled surface melt and the location of ablation stakes (black dots) for the study period 20 June to 12 September 2013.

Figure 9 .
Figure 9. Observed (measured) melt from ablation stakes and modelled melt from the DEBM.

Figure 10 .
Figure 10.Ablation due to calving and surface melt (below the ELA) at Bridge Glacier during the 2013 melt season.
calving fluxes were small in 2001-2002.Calving fluxes were high between 2003 and 2006 (0.030-0.084 km 3 a −1 ) and from 2008 to 2011 (0.036-0.100 km 3 −1 ) with low calving rates in 2006-2007.As the calving flux increased from 2003 to 2011, surface melt below the ELA decreased slightly, resulting in the calving flux becoming a larger component of the total ablation in the 21st century.Ablation due to calving was roughly equal to surface melt below the ELA in 2005, 2008, and 2010 (44-49 % of total ablation).

Figure 11 .
Figure 11.Historical ablation from calving and surface melt (below the ELA), 1984-2013.The dark vertical line in 1991 indicates the period in which the terminus reached flotation and calving rates increased.Shaded areas correspond to calculated uncertainty.

Table 2 .
Characteristics of selected major lake-calving glaciers worldwide.D w is the mean water depth, T w is the mean water temperature (depth averaged or interannual range), U T is the terminus averaged flow speed, and U c is the calving rate.T w ( • C) U T (m a −1 ) U c (ma −1 ) Source