Initiation of a major calving event on the Bowdoin Glacier captured by UAV photogrammetry

In this paper, we analyse the calving activity of the Bowdoin Glacier, north-western Greenland, in 2015 by combining satellite images, UAV (unmanned aerial vehicle) photogrammetry and ice flow modelling. In particular, a highresolution displacement field is inferred from UAV orthoimages taken immediately before and after the initiation of a large fracture, which induced a major calving event. A detailed analysis of the strain rate field allows us to accurately map the path taken by the opening crack. Modelling results reveal (i) that the crack was more than half-thickness deep, filled with water and getting irreversibly deeper when it was captured by the UAV and (ii) that the crack initiated in an area of high horizontal shear caused by a local basal bump immediately behind the current calving front. The asymmetry of the bed at the front explains the systematic calving pattern observed in May and July–August 2015. As a corollary, we infer that the calving front of the Bowdoin Glacier is currently stabilized by this bedrock bump and might enter into an unstable mode and retreat rapidly if the glacier keeps thinning in the coming years. Beyond this outcome, our study demonstrates that the combination of UAV photogrammetry and ice flow modelling is a promising tool to horizontally and vertically track the propagation of fractures responsible for large calving events.


Introduction
In the last decades, many ocean-terminating outlet glaciers of the Greenland ice sheet (GIS) experienced thinning and rapid retreat (Joughin et al., 2010;Pritchard et al., 2009), which contributed to the global loss of ice and affected sea level rise.While the general retreat is triggered by warming temperatures, local conditions determine the timing and pace of retreat (Kjeldsen et al., 2015).For instance, the tidewater glaciers of north-western Greenland started to rapidly retreat only in recent years, i.e. about 20 years after glaciers of the south.Approximately half of the ice ablation of the GIS is due to calving, i.e. the release of icebergs at the edge of glaciers (Enderlin et al., 2014).The calving mechanism is still not entirely understood, mostly because of the complex interconnection between the involved processes (Benn et al., 2007b).Among them, the sharp acceleration of the ice flow at the calving front in response to high buoyant forces generates deep crevasses, which precondition the future break off of icebergs.Since calving events are triggered by fracturing processes, calving is often considered a random process, which must be analysed over large data sets of events.Therefore, relatively few studies have been dedicated to a description of calving at the level of individual events (e.g.O'Neel et al., 2003;Chapuis and Tetzlaff, 2014, and references therein).However, observations in Greenland (e.g.Medrzycka et al., 2016) show that large-scale, occasional calving events often contribute more to total calving loss than small-scale frequent events.
Calving primarily results from the propagation of fractures upstream of the calving front in response to high stresses ( Van der Veen, 1998).Opening and sustained growth of cracks occur when the normal and the shear stress components exceed a certain threshold (Colgan et al., 2016).Yet, the normal component is often assumed to be a primary control over the shear one (Benn et al., 2007b).For a given normal strain rate, Nye (1957) derived a function for calculating the depth of a crevasse by assuming that it will extend until the normal strain rates balance the creep closure rate resulting from the ice overburden pressure.Based on this function, Benn et al. (2007a) derived a calving law by assuming that calving occurs when a crevasse penetrates the water line.Earlier, Vieli et al. (2001) proposed another criterion based on the height-above-buoyancy variable.However, these semiempirical approaches assume closely spaced crevasses and do not include the effect of the stress concentration at the tip of cracks.To overcome this problem, one has to implement more complex continuous approaches based on linear elastic fracture mechanics ( Van der Veen, 1998), or continuum damage mechanic (Pralong and Funk, 2005), which merge in the same framework the micro-scale formation of cracks and the macro-scale behaviour of the damaged ice flow.Despite a large number of contributions on calving modelling in recent years, relatively few studies have attempted to model calving events individually (Åström et al., 2013); this is presumably because microscopic approaches can hardly be coupled to traditional ice flow, unlike macroscopic approaches.
The main obstacle to implementing continuous models is the lack of data to constrain the parameters related to the fracturing of ice (Pralong and Funk, 2005).In particular, no direct measurements of failure stress values within the ice are available except for laboratory ones which, however, cannot reproduce the observed orders of magnitude.By contrast, critical strain rates can be measured.Although 0.01 yr −1 is often given as reference threshold (Meier, 1958), observed threshold normal strain rates for the initiation of crevasses span over 2 orders of magnitude (Colgan et al., 2016).The estimation of critical strain rates has been rendered easier through the use of automatic camera and aerial photogrammetry combined with feature-tracking techniques (Colgan et al., 2016;Messerli and Grinsted, 2015).In particular, recent developments of unmanned aerial vehicles (UAV) and photogrammetry by structure-from-motion allows high-resolution ice flow velocity fields and associated strain rates to be inferred (Ryan et al., 2015).
In this paper, we use UAV photogrammetry, feature tracking and ice flow modelling to analyse, in detail, the propagation of a major crack on the Bowdoin Glacier, north-western Greenland, before it collapses and generates a major calving event representing about 14 % of the annual calving in terms of area.From this analysis, we investigate the calving pattern at the terminus of the Bowdoin Glacier in May and July-August 2015 observed on satellite images.This paper is organized as follows.First, we describe our data and the methods we have employed to lead this study.Then, we present our results related to satellite images, UAV-inferred ice flow field, strain analysis and ice flow modelling.Finally, we use these results to revisit and explain the calving pattern of the Bowdoin Glacier observed in 2015.
2 Data and methods

Study site
The Bowdoin Glacier (77 • 41 N, 68 • 35 W) is an oceanterminating glacier, which belongs to a network of outlet glaciers located in the north-west of the Greenland ice sheet.It is approximately 10 km long and 3 km wide with an average surface slope of less than 1 • .The glacier discharges into the Bowdoin Fjord through a 3 km wide, most-likely grounded calving front (Sugiyama et al., 2015).At the centre of the calving front, the Bowdoin Glacier was about 250 m thick, and its maximal flow speed was about 1.5 m day −1 in 2013.A medial moraine, that we often use as a reference line in what follows, runs parallel to and about 1 km away from the left glacier margin (see Fig. 1).
The front of the Bowdoin Glacier has been fairly stable since the end of the 19th century.Its present-day position is no more than 2 km upstream of its 1897 position (Chamberlin, 1897).However, from July 2008 to September 2013, a rapid retreat of the calving front at a rate of 220 m yr −1 was recorded, while no significant changes occurred during the previous 20 years (Sugiyama et al., 2015).From 2013, the calving front has stabilized to about its current position.The rapid retreat recorded during the period 2008-2013 was attributed to a local depression in the bedrock, which is known to make the grounding line unstable (Schoof, 2007).
In July 2013, 2014 and 2015, we carried out field measurements on the Bowdoin Glacier.This included ice radar and GPS measurements (Sugiyama et al., 2015;Tsutaki et al., 2016), seismic records (Podolskiy et al., 2016), automatic camera installations and borehole drilling (to record in-ternal ice deformation, englacial temperature and water pressure).

Borehole temperature and deformation data
A vertical temperature profile was obtained from borehole thermistors installed 2 km upstream of the calving front, see Fig. 3. Temperature measurements (unpublished data) were averaged over the 9-month period from 14 October 2014 to 14 July 2015 after complete diffusion of the local temperature perturbation caused by the hot water drilling on 20 July 2014.The resulting profile shows ice at −10 • C at the surface, −3 • C to a depth of 25 m, −6 • C to the depths of 90 and 150 m, −3.5 • C to a depth of 210 m and 0 • C at the 260 m deep bedrock.As another relevant outcome for this study, we found that the ice deformation represents less than 10 % of the total motion in the borehole (unpublished data).

Satellite images
To analyse the large-scale calving pattern and ice flow field of the Bowdoin Glacier in 2015, we used the Landsat 8 Operational Land Imager panchromatic images.In particular, these data allowed us to infer the positions of the calving front (see Fig. 1) and estimate the calved surface area during the largest events.

UAV photogrammetry
During the 2015 field campaign, we flew an unmanned aerial vehicle over the terminus of the Bowdoin Glacier to collect aerial images and reconstruct, in high resolution, the 3-D geometry of the calving front using photogrammetry.Although we closely followed the methods described by Ryan et al. (2015), we shortly review the techniques we have employed, from UAV implementation to image post-processing.
For the UAV, we used the 2 m wide fixed-wing "Skywalker X8" equipped with the "Pixhawk" open-source autopilot (https://pixhawk.org/)running with the APM "ArduPlane" firmware (http://ardupilot.org/ardupilot/).The UAV also carried a camera (Sony Alpha 6000), which had a 24 megapixel sensor and a 16 mm lens.The UAV flew autonomously twice on 11 and 16 July, following a pre-programmed sequence of waypoints located at the two extremities of the glacier front so that the calving front was covered by four parallel flight lines.For each flight, the UAV flew 25 km at an altitude of about 300 m above the ground and collected about 1000 overlapping pictures of the calving front, with a resolution of 7 cm per pixel.An overlap of 95 % in flight direction and 70 % in cross-flight direction was chosen to obtain a stable aerotriangulation.For reliable georeferencing, Ground Control Points (GCPs) were installed on the left border of the glacier and on the moving glacier surface.The position of the moving GCPs was recorded repeatedly with a static differential global positioning system (DGPS) so that their absolute position at the time of each flight was linearly approximated.
The pictures of each flight were post-processed through the software Agisoft PhotoScan (http://www.agisoft.com/).Using the image correlation techniques (structure-from-motion) of Photoscan, a high-resolution 3-D model and orthoimage (covering about 2 km 2 ) of the calving front was reconstructed for each flight (see Fig. 2a and b).

Feature tracking
The orthoimages derived from satellite and UAV data are used to infer the surface displacement fields of the ice by using a feature-tracking method.There exist a variety of image matching methods to derive glacier surface speed (Heid and Kääb, 2012).Here, we have used an algorithm of crosscorrelation in frequency domain on orientation images (Abe et al., 2016;Heid and Kääb, 2012;Sugiyama et al., 2015).On the one hand, a 30 m resolution velocity field, averaged between 14 June and 1 August (see Fig. 3), was inferred from satellite images.On the other hand, from a native resolution of 10 cm for the UAV orthophoto, a template searching chip of 200 × 200 pixels and a step number of 10 × 10 pixels, we obtained a 1 m resolution velocity field averaged between 11 and 16 July (see Fig. 2c).

Strain analysis
The first-order control for triggering calving events is related to the strain rates (Benn et al., 2007b).For this reason, we computed the horizontal part of the strain rate tensor: where the derivatives are approximated by finite difference from the horizontal velocity field (u x , u y ) inferred by feature tracking.However, D depends on the system of coordinates (x, y) and thus cannot be interpreted directly.Instead, we computed the eigenvalues of D since they do no depend on the system of coordinates.The highest eigenvalue (called the maximal principal strain) indicates the maximal normal strain rate among all possible directions, while its associated eigenvector (called the maximal principal direction) is the direction in which the normal strain is maximized (see Fig. 2d).

Modelling
The UAV can only capture the ice dynamics at the glacier surface.To gain insight into the dynamical processes within the glacier, the terminus of the Bowdoin Glacier was modelled using the Elmer/Ice code (Gagliardini et al., 2013).
In fact, we modelled an enlarged domain covering the last 4 km of the glacier and a reduced domain which only encompasses the last few hundred metres, see Fig. 3.For each domain, we built a 3-D mesh of the glacier front from basal and surface ice topographies using the GMSH mesher (Geuzaine and Remacle, 2009).The basal topography was www.the-cryosphere.net/11/911/2017/The Cryosphere, 11, 911-921, 2017 July, with the resulting velocity field inferred using feature tracking (c) and the maximum principal directions of the strain rate (d), respectively.For the sake of visualization, only the maximum directions with a magnitude above 5 yr −1 are drawn on (d).The vector lengths scale with the magnitude of the maximal principal strain rate.A scale corresponding to 100 yr −1 is shown as reference.Blue and red arrows stand for principal directions on the upstream and downstream edge of the crack, respectively.
extrapolated with a parabolic profile from ice radar measurements (Sugiyama et al., 2015), while the surface topography was derived from the UAV-inferred 3-D model and completed by a larger available Digital Elevation Model.
Elmer/Ice implements a full-Stokes model (Greve and Blatter, 2009), and ice is considered as a non-Newtonian fluid governed by Glen's flow law: where and τ denote the strain and deviatoric stress tensor, τ II denotes the second invariant of τ , A is the temperaturedependent Arrhenius factor, A(T ) (Paterson, 1994), and E is an enhancement factor which controls the stiffness of ice.A relationship between englacial temperature and depth was interpolated from borehole temperature measurements (Sect.2.2) and generalized over the whole modelled do-main.For both (enlarged and reduced) modelled domains, the model was supplemented by the following boundary conditions.No force applies to the top ice surface, while the following condition applies at the calving front.
where σ ij is the Cauchy stress tensor, ρ w denotes the density of water, g is the gravity and (n x , n y , n z ) is the outer normal vector.At the back and lateral boundaries of the modelled domain, we impose the Dirichlet boundary condition, where u i,meas is the UAV-inferred velocity field at the glacier surface, b is the bedrock elevation, s is the ice surface elevation, f = (ρ w /ρ i )×(−b/ h) is the flotation ratio (dimensionless) -which is 0 when no buoyant force applies and 1 when ice is floating (this is never the case in the present study), ρ i denotes the density of ice and h = s − b is the ice thickness.Equation (4) says that the horizontal velocity decreases linearly with depth to reach the fraction f of the surface motion at the bedrock.Note that we have f ∼ 90 % at the borehole, in agreement with the measurements.At the glacier bed, we apply Budd's friction law (Budd et al., 1984), which links basal shear stress (τ b ) to basal sliding (u b ) through the relation where C > 0 is a constant sliding coefficient and α > 0 is a given exponent.Since only the terminus of the glacier was modelled, the buoyancy under the glacier bed could be reasonably calculated as a simple function of depth (z) assuming a fully distributed hydrological system.In that case, the effective pressure (N), which is defined by the bed normal stress in the ice minus the buoyancy, is related to the flotation ratio (f ) through the relation N = ρ i gh(1 − f ), approximating the bed normal stress as hydrostatic.Thus, Eq. ( 5) can be understood as follows: a high buoyancy renders the effective pressure (N) or (1 − f ) small and favours sliding.Parameter α in Eq. ( 5) therefore controls the degree of influence of the flotation ratio on sliding from no-influence when α = 0 (Weertman's law) to a quadratic influence when α = 2.In this study, we restrict ourselves to the range α ∈ [0, 2], as tuned values reported in the literature (e.g.Budd et al., 1984;Nick, 2006) most often remain in this interval.
The model was run in two steps.First, the ice flow parameters E, α and C in Eqs. ( 1) and ( 5) were tuned so that (i) the modelled surface velocities match the measured ones as well as possible and (ii) the ratio between vertical ice deformation and total ice motion at the borehole is about 10 %, as measured in 2014-2016 by borehole inclinometers (see Fig. 4, top panel).Since the borehole is located about 2 km upstream of the calving front (see Fig. 3), we first modelled an enlarged domain and used the satellite-inferred velocity field (see Fig. 3) instead of the UAV-inferred ones for boundary conditions since it covers a larger area.Once the three ice flow parameters were tuned, we restricted our domain to the glacier front and optimized it (to improve the consistency between modelled and measured velocity fields) according to the two geometrical parameters (h lift and d frac ) of the glacier front (see Fig. 4, middle and bottom panels).First, the bedrock topography was lifted laterally, from the moraine position to the south-eastern glacier side, to mimic a bump presumed responsible for the sharp transition between fast and slow ice flow.The height of this bump is denoted h lift .Second, a cut of constant depth (d frac ) was applied over the surface at the location of the main fracture to simulate the real depth, which cannot be inferred from aerial images.It must be stressed that this cut directly affects the mesh.Thus, opposed hydrostatic forces (boundary condition, Eq. 3) apply to both sides of the fracture (assuming it is filled with water up to sea level, see Sect. 4) with the result of intensifying its opening.
Lastly, we again performed the optimization of parameters C, h lift and d frac as described above but for a fixed value α = 1, different enhancement factors E ∈ {1, 2, 4, 8}, and ignoring the 10 % measured vertical ice deformation data, see Fig. 5.This last experiment serves to assess the impact of the uncertainties related to the Arrhenius factor, A(T ), on our results.These uncertainties might be caused by anisotropic ice, damage ice or errors in the temperature profile.

Calving pattern
The analysis of the 2015 satellite images reveals that the substantial ablation by calving was due to a few monthly spaced calving events.In particular, four May and July-August events (which occurred on 3 May, 19 May, 27 July and 9 August) contributed to 0.13, 0.13, 0.17 and 0.07 km 2 of surface loss, respectively (see Fig. 1).These numbers must be compared to the annual loss by calving that we estimate to be 1.2 ± 0.2 km 2 .This estimate was obtained by integrating the ice flux across the calving front from the satellite-inferred velocity field (see Fig. 3), with the calving front being fairly stable since 2013.In total, the four events contributed to about   40 % of the annual calving.Intriguingly, the events of May and July-August show similar patterns: a first slice of ice splits on the south-east side, however, without touching the glacier border and a second slice separates further north-west tens of days later (see Fig. 1).

Calving event of 27 July
During the summer 2015 field campaign, field observations and UAV images allowed us to reconstruct, step-by-step, the early stages preceding the calving event, which occurred on 27 July.On 11 July, the UAV orthoimage shows no major fracture in the vicinity of the calving front and gives no signs that a major event is about to occur (see Fig. 2a).One day later, team members of the field campaign reported that a few-metres-wide crack appeared about 100 m upstream of the calving front across the medial moraine.On 13 July, the crack was visible from an automatic camera, which monitors the calving front from about a 2 km distance.On 16 July, the UAV orthoimage clearly indicates that the crack was about 750 m long (see Fig. 2b).On 27 July, about a 1 km long slice of the front (i.e. about one third) collapsed (see Fig. 1).

Velocity and strain rate fields
Figures 3 and 2c show the horizontal ice flow velocity fields inferred by feature tracking from satellite and UAV orthoimages, respectively.Although the satellite-inferred data covers a larger domain and time period than the UAV-inferred data, we obtain a good agreement where the two results overlap.Interestingly, the UAV result shows a clear discontinuity along the crack, which was identified on the orthoimage of 16 July (see Fig. 2b).Additionally, Figs. 3 and 2c reveal a sharp transition in the ice flow speed: homogeneously fast flow in the central section and very slow motion in the southeastern part, within about 700 m from the glacier's margin.Finally, Fig. 2d displays the maximum strain rate's principal directions computed from the velocity field.As a matter of fact, the main fracture and the highly damaged zone delimiting the slow and fast bands concentrate the highest maximal principal strain rates.

Modelling results
First, we calibrated the ice flow and sliding parameters, E, C and α (related to Eqs. 1 and 5), to simultaneously minimize the misfit between modelled and measured surface velocities and to match the deformation rate of ice recorded at the borehole.In fact, we found that the enhancement factor (E) could be tuned independently of the other parameters, and the value E = 4 proved to give a consistent deformation rate at the borehole (i.e. about 10 % of vertical ice deformation), irrespective of the other parameters.For each α ∈ {0, 0.5, 1, 1.5, 2}, we found a unique parameter C minimizing the misfit between modelled and measured surface velocities (see Fig. 4, top panel).However, it was not possible to clearly identify the best couple (C, α), since all combinations give comparable misfits (about 0.3-0.4m per day).Finally, we found a good linear fit between α and the logarithm of the sliding coefficient (C) through the following empirical relationship.
which is in the same order of magnitude as the value C = 6.06 × 10 4 × h α Pa m −1/3 s 1/3 obtained for Columbia Glacier in Nick (2006), and α = 3.5/3 when the ice thickness (h) spans between 100 and 300 m.
Second, we optimized the bump height (h lift ) on the southeastern side of the glacier front (see Fig. 6) to best reproduce the high-shear zone shown by the UAV-inferred velocity field (see Figs. 2c and 7). Figure 4 (middle panel) displays the misfit between modelled and measured surface velocities as a function of α (and the optimal C given by Eq. 6) and the bump height h lift .Note that the measure of the misfit excludes the detaching slice of ice to better focus on the high-shear zone.Interestingly, as in the previous optimization, one cannot clearly constrain α from this optimization since several values give comparable misfits (see Fig. 4, mid- dle panel).However, we observe that bumps are slightly thinner, with high values of α.Indeed, basal motion contributes to most (90 %) of the total motion so that a bump on the bed further slows down the ice flow when the influence of the reduced effective pressure (N) on basal sliding is the highest (i.e. for high α).As a result, we found that an h lift = 125 m high bump was necessary to generate a high-shear zone comparable to the one derived by UAV when using the mean α = 1, with an uncertainty of ±25 m when accounting for the limiting cases α = 0 and α = 2.After including such a bump on the bedrock, the misfit dropped from 0.2-0.3 to 0.05 m per day (as compared to the first optimization) for any value of α.
Third, we optimized the depth of the main fracture (d frac ) to best reproduce the discontinuity shown by the UAVinferred velocity field at the crack location, assuming an h lift = 125 m high bump over the bedrock.Figure 4 (bottom panel) displays the misfit between modelled and measured surface velocities as a function of α (and C given by Eq. 6) and the fracture depth (d frac ).The measure of the misfit is restricted to the area near, and downstream from, the crack.As a result, only a 175 ± 35 m deep crevasse shows a result consistent with UAV-inferred observations (see Figs. 4,bottom panel,and 7).Interestingly, this result is not too affected by the considered sliding law parameterization (α and C).As a matter of fact, including the crevasse depth in the model reduced the misfit by about 65 % (see Fig. 4).Additional modelling setups based on non-constant crack depths (d frac ) were run, see Fig. 8 bottom panel.However, constant depth runs yielded to the best match between modelled and measured velocities (results not shown).
Lastly, Fig. 5 shows that the uncertainties related to the Arrhenius factor, A(T ) (caused, for instance, by errors in the temperature profile), do no impact our results much: less than 10 m for the bump height h lift and 20 m for the crack depth d frac , not accounting for the case where E = 8, which is beyond any realistic enhancement factor.

Discussion
The orthoimage of 16 July reveals a main fracture (Fig. 2b); however, it is not clear whether the fracture reaches the calving front or not, due to the presence of seracs at the glacier surface.The discontinuity in the velocity field (Fig. 2c) along the fracture settles this issue confirms that the fracture extends to the front and accurately describes its path.This discontinuity suggests that a slice is detaching from the glacier on the south-eastern side and remains only attached to the other side.Corroboratively, Fig. 2d shows very high maximal principal strain rate components near the crack (between 50 and 150 yr −1 ).As a matter of fact, the directions of maximal strain along the crack are not always perpendicular to the fracture, but tilt up to 45 • toward the centreline of the glacier.This indicates that without control from the southeastern side, the entire slice tends to rotate around its extremity, which is still anchored.This is confirmed by modelling results: only about a 175 m deep and water-filled crevasse (i.e. about two thirds of the ice thickness) can reproduce the observed jump in velocities (see Fig. 4, bottom panel, and Figs.7 and 8).Additional modelling results showed that the crack depth has to be rather constant along the crack path.Because the crack was connected to the sea (see Fig. 2a), the water level within the crevasse most likely coincided with the sea level.In that configuration, we can estimate the pressure exerted by ice and water at the bottom of the crack to about 1.56 and 1.42 MPa, respectively.On the other side, the tensile stress must have been in the order of magnitude of tensile strengths, i.e. around 0.8 ± 0.4 MPa for clean temperate ice (Schulson and Duval, 2009).Yet, the crack is expected to deepen further if the pressure of water plus the tensile stress exceed the pressure of ice (Van der Veen, 1998).As a result of this delicate balance, we can reasonably assume that the crack propagated downward after 16 July in response to small perturbations, like, for instance, a slight increase in the water level within the crevasse.Such an increase could have been triggered by the intense meltwater production during the first half of July, consistently to warm temperatures, and the cloudless sky observed on the field.The deepening of the crack must have further increased the stress concentration at the crack tip and caused the crack to extend further laterally before the collapse (see Fig. 1).Unfortunately, we cannot corroborate our crack depth value with the formula given in Benn et al. (2007b), which provides an analytic expression for the crevasse depth as a function of the tensile strain and the depth of water in the crevasse, since it does not account for stress concentration effects, which are significant in the present situation.
The velocity field displayed in Fig. 2c clearly shows that there is a highly sheared zone located upstream of the junction between the crack and the calving front, with the velocities experiencing a strong variation between slow/resistive lateral and fast/floating central flow.Such a pattern is typical for glacier margins, where the lateral drag induced by the fjord's buttressing is the highest.However, this zone is not directly in contact with the glacier margin but is instead located between 500 and 1000 m away from it.In addition, this area is characterized by many large crevasses, which are 45 • tilted with respect to the flow direction, a typical characteristic of highly sheared margin zones (Paterson, 1999;Colgan et al., 2016).Such a situation can only be the result of a strong irregularity of the bedrock or of the basal conditions.Ice radar measurements of the basal topography across a profile about 800 m upstream of the calving front indicate a rather flat bed in the central part and a constant slope profile between the medial moraine and the margin (Sugiyama et al., 2015).Such configuration might cause transversal gradients of the velocity (the central part of the glacier being nearly floating while the margin experiences more basal resistance), but not to the extent shown by Fig. 2c.Unfortunately, no information about the bedrock is available near the calving front since it is far too crevassed to make ice radar measurements.Yet, our modelling results (see Figs. 4 and 7) show that a bump on the bed of almost half of the glacier thickness (about 125 m) from the moraine position to the south-eastern glacier side can cause the observed shear on its own (see Figs. 8 and 6).Indeed, a basal bump would substantially reduce the buoyancy, the basal sliding and the surface motion.Other field evidence corroborates this theory.First, the presence of an up-to 20 m high bump in the centre of this area of the glacier surface (see Fig. 6) probably reflects a feature in the basal topography.Interestingly, this surface hump has already been reported in Chamberlin (1897), presumably because this intriguing feature was further pronounced at the end of the 19th century due to thicker ice: "The slope just beyond the crevassed area [. . . ] and between that area and the medial moraine inclines backward [. . .] and a brook flows in that direction.The notable feature of the surface is the crevassing."Finally, the border side of the glacier is characterized by a ridge, which is oriented toward the sheared area (see Fig. 6).If the ridge extends under the ice, the bedrock is expected to be shallower there than elsewhere.An alternative explanation for the high-shear zone would be a locally frozen bed preventing basal sliding near the glacier margin.However, such an asymmetry in basal temperatures would likely be the result of an asymmetric bedrock.As a consequence, an intermediate case of a shallower but frozen bump over the bedrock cannot be excluded, so the 125 m bump height found in modelling experiments must be understood as an upper bound.
The presence of a freshwater plume at the level of the moraine (see Fig. 2a) raises the question of its role as a potential trigger of the 27 July calving event.Indeed, the plume indicates the presence of a freshwater discharge at the foot of the calving front.Yet, plumes are characterized by turbulent flow and enhanced submarine melting along the cliff.In such conditions, the calving face below sea level might have been overhanging, causing the front to lean forward and favour the opening of the crevasses (Kimura et al., 2014).However, additional simulations (not shown) accounting for local overhanging of the calving front did not allow us to establish any clear link between the presence of a plume and the opening of the main crack.Indeed, the observed jumps in the velocity field are remarkably uniform along the crack (see Fig. 2d), while melting the foot of the calving front locally would render the modelled jumps of velocity strong near the plume but weaker elsewhere.
Based on the analysis of the three previous paragraphs, we can now revisit chronologically the proceeding of the mechanisms, which led to the calving event of 27 July 2015.First, a shallow and possibly frozen bedrock near the calving front between the south-eastern margin and the medial moraine (see Fig. 8) slows down the ice flow in this area.In response to high resulting strain rates (see Fig. 2d), many crevasses form at the transition between the slow and fast zones.In turn, this destabilizes the vicinity of the calving front and a major crack appears there on July 12.When in contact with the sea, tidewater enters the crack and triggers an irreversible deepening, with the cryostatic pressure no longer being able to balance the water pressure at the bottom of the crevasse.The crack deepening gradually detaches the slice of ice from the glacier trunk and causes some stress concentration at the crack tip.In turn, this causes the crack to propagate laterally over 1 km to the north-western part about 100 m upstream of the front before the collapse on 27 July (see Fig. 1).A remaining chunk collapses on 9 August during a second major calving event.This proceeding of calving events is very similar to the set of events which occurred in May (see Fig. 1), and we can reasonably assume that events of this type will occur again in the future as long as the calving front remains over the presumable bump.However, the Bowdoin Glacier has experienced a rapid retreat of ca. 2 km between 2007 and 2013 before stabilizing at the current position despite a thinning at a rate of ca. 4 m yr −1 (Tsutaki et al., 2016).Thus, we assume that the shallow and possibly frozen bedrock between the medial moraine and the south-eastern margin has been playing a crucial role in stabilizing the calving front.Consequently, if the glaciers keeps thinning in the future, one must expect a rapid, unstable retreat of the calving front shortly after this zone is overpassed because of a likely bedrock overdeepening.

Conclusions
Using UAV photogrammetry, we captured high-resolution orthoimages of the calving front of the Bowdoin Glacier before and after the initiation of a large fracture which induced a major calving event.A detailed analysis of the diswww.the-cryosphere.net/11/911/2017/The Cryosphere, 11, 911-921, 2017 placement field and the resulting strain rates allowed us to accurately reconstruct the path taken by the crack.Combined with modelling results, we could determine that the crack was more than half-thickness deep, filled with water up to sea level and getting irreversibly deeper when it was captured by the UAV.Later on, the crack deepening caused stress concentration around the tip, with the crack propagating laterally, and finally a collapse.Modelling results also indicated that the crack was likely triggered in a highly crevassed area near the front caused by a shallow and possibly frozen bedrock.The asymmetry of basal conditions at the front explains the calving pattern observed in May and July-August 2015, while symmetric conditions would have rendered calving events less predictable.More importantly, our results indicate that the calving front of the Bowdoin Glacier is likely stabilized by a shallow bedrock under the south-eastern glacier margin.As a consequence, the glacier might pass over a basal depression if the glacier keeps thinning, so that a rapid unstable retreat of the front must be expected in the next years.
In this paper, we have used UAV photogrammetry and feature tracking by closely following the techniques described by Ryan et al. (2015).By combining this approach to ice flow modelling, we additionally analyse, in detail, the horizontal and vertical propagation of large water-filled crevasses near the calving front before they collapse and generate a calving event.This approach is especially relevant for large calving events given that few of them can contribute more to the global ice ablation than small and frequent ones (Medrzycka et al., 2016).However, only two snapshots of the entire fracturing process were available in the present study.In particular, this is not enough to determine whether the crack kept propagating vertically after 16 July or if the stress concentration at the tip was sufficient to allow for a lateral extension as well.By contrast, frequent UAV-inferred orthoimages and the resulting flow fields between the first crack initiation and the final collapse of this particular event would have certainly allowed us to accurately infer the whole proceeding.As a consequence, future studies focusing on monitoring individual crevasses prior to iceberg calving should include repetitive UAV flights to capture all phases associated with the event.
Data availability.Data are available from the authors upon request.
Author contributions.Guillaume Jouvet designed the study, ran the Elmer/Ice model and wrote the paper with support from Julien Seguinot, Martin Funk and Shin Sugiyama.Guillaume Jouvet and Yvo Weidmann assembled, tuned and flew the UAV during the field campaign in 2015.Yvo Weidmann processed the aerial images by photogrammetry.Takahiro Abe and Daiki Sakakibara inferred the velocity fields and the calving front positions from UAV and satellite orthoimages.Hakime Seddik derived the Digital Ele-vation Model of the bedrock from ice radar measurements.Martin Funk and Shin Sugiyama led and organized the field campaign on the Bowdoin Glacier in July 2015.All authors contributed to the paper.

Figure 1 .
Figure 1.Positions of the calving front of the Bowdoin Glacier inferred from satellite images before and after the calving events of 3 May, 19 May, 27 July and 9 August.

Figure 2 .
Figure 2. Calving front of the Bowdoin Glacier: orthoimages were obtained by UAV on (a) 11 and (b) 16July, with the resulting velocity field inferred using feature tracking (c) and the maximum principal directions of the strain rate (d), respectively.For the sake of visualization, only the maximum directions with a magnitude above 5 yr −1 are drawn on (d).The vector lengths scale with the magnitude of the maximal principal strain rate.A scale corresponding to 100 yr −1 is shown as reference.Blue and red arrows stand for principal directions on the upstream and downstream edge of the crack, respectively.

Figure 3 .
Figure 3. Average velocity field between 14 June and 1 August inferred by feature tracking from LANDSAT satellite images.The back dot indicates the position of the borehole.The dashed and continuous lines delimit the enlarged and reduced modelled domains, respectively.

Figure 4 .
Figure4.Level sets of the absolute misfit (in m per day) between modelled and measured surface velocity magnitudes as a function of the parameters α (y axis), C (x axis, upper panel), h lift (x axis, middle panel) and d frac (x axis, bottom panel).Each black dot corresponds to one model realization from which the level sets were extrapolated.For each value of α, the symbol + displays the minimum with respect to C, h lift or d frac , while the dashed line connects the + symbols.

Figure 5 .
Figure 5. Level sets of the absolute misfit (in m per day) between modelled and measured surface velocity magnitudes as a function of the enhancement factors E (y axis) and h lift and d frac (x axis), respectively.Here, we used a constant parameter α = 1, and C was tuned to measured surface velocities in the extended model domain.The symbols have the same meaning as in Fig. 4.

Figure 6 .
Figure 6.Contour lines of the initial bedrock topography (a), the optimized bedrock topography (b) and the digital elevation model (c) in the neighbourhood of the calving front.The thick line in panel (b) indicates the reduced modelled domain in which the bedrock was optimized.Panel (c) reveals the presence of a 20 m high bump (shown by the symbol) on the glacier surface and of a ridge on the south-eastern side of the glacier.

Figure 7 .
Figure 7. Magnitude of the modelled velocity field after the first (a), the second (b) and the third (c) optimization, i.e.(a) corresponds to the original setup with optimal parameters (E, C, α), (b) accountsfor the optimal bump height on the bed while (c) accounts for the optimal crack depth.For the sake of comparison, (d) displays the measured UAV-inferred velocity field on the surface.

Figure 8 .
Figure 8. Top: inferred position of the crack on the orthoimage of 16 July (dashed line).Bottom: bedrock and ice surface along the main crack starting from the tip before and after optimizing the bump height and the crack depth to the UAV-inferred ice flow velocities and tested varying crack depths.The dotted rectangle indicates the position where the crack crosses the moraine.