Automatic detection of calving events from time-lapse imagery at Tunabreen, Svalbard

Calving is an important process in glacier systems terminating in the ocean and more observations are needed to improve our understanding of the undergoing processes and be able to parameterise calving in larger scale models. Time-lapse cameras are good tools for monitoring calving fronts of glaciers and they have been used widely where conditions are favourable. However, automatic image analysis to detect and calculate the size of calving events has not been developed so far. Here, we 5 present a method that fills this gap using image analysis tools. First, the calving front is segmented. Second, changes between two images are detected and a mask is produced to delimit the calving event. Third, we calculate the area given the front and camera positions as well as camera characteristics. To illustrate our method, we analyse two image time series from two cameras placed at different locations in 2014 and 2015 and compare the automatic detection results to a manual detection. We find a good match when the weather is favorable but the method fails with dense fog or high illumination conditions. 10 Furthermore, results show that calving events are more likely to occur (i) close to where subglacial melt water plumes have been observed to rise at the front and (ii) close to one another.


Introduction
Tidewater glaciers are one of the main contributors to sea-level rise (Church et al., 2013;Gardner et al., 2013)  2013). Chapuis and Tetzlaff (2014) studied the calving events sizes of Kronebreen, Svalbard, and concluded that calving events variability is inherent to calving dynamics even under stable external conditions. Event-size distribution follows a scale-invariant power-law that has been further discussed by Åström et al. (2014), who classify termini of calving glaciers as self-organised critical systems such as earthquakes. They recommended to focus on quantifying the effects of external forcing on the critical state of calving margins, which makes event size and frequency analysis from observation necessary. 5 Time-lapse cameras are convenient monitoring tools that have been used for different purposes in glaciology such as daily digital elevation models (James et al., 2014), propagation of flexion zones up-glacier (Murray et al., 2015), glacier velocities (Ahn and Box, 2010;Messerli and Grinsted, 2015;James et al., 2016), monitoring of supraglacial lake levels (Danielson and Sharp, 2013) and meltwater plume surface area (How et al., 2017) or glacier surges (Kristensen and Benn, 2012). Time-lapse imagery has also been used to quantify calving events manually, with different resolution and different scaling methods, for 10 2 Geosci. Instrum. Method. Data Syst. Discuss., https://doi.org /10.5194/gi-2018-5 Manuscript under review for journal Geosci. Instrum. Method. Data Syst. Discussion started: 2 July 2018 c Author(s) 2018. CC BY 4.0 License.
Columbia glacier, Alaska (Walter et al., 2010), Paierlbreen, Svalbard (Åström et al., 2014), Tunabreen, Svalbard (Åström et al., 2014;Westrin, 2015) and Rink Isbrae, West Greenland (Medrzycka et al., 2016). Alternative approaches have also been used to estimate the size and frequency of calving events at tidewater termini such as direct visual detection (e.g. O'Neel et al., 2007;Bartholomaus et al., 2012;Chapuis and Tetzlaff, 2014), icequake detection (e.g. Bartholomaus et al., 2012;Köhler et al., 2016) or satellite imagery (e.g. Schild and Hamilton, 2013). 5 To manually detect calving events from time-lapse cameras is a very laborious task that requires days of work and becomes too difficult when the amount of images is too large. For example, if the manual processing of an image takes two minute, the operator would need 17 days (8 hours of work per day) to process a month of images with 10 minutes interval. As a consequence, few studies use the potential of time-lapse cameras and do not analyse more than a few days. Automatic methods have the advantage to enable longer time series but, to date, no satisfying automatic calving detection from time-lapse cameras 10 have been developed.
Here, we utilize a method, first presented by Adinugroho (2015a) and Adinugroho et al. (2015b), to automatically detect calving events from sequential time-lapse imagery. The algorithm follows five main steps in order to achieve this: image registration, segmentation of the calving front, change detection, mask reconstruction and size calculation. The results from this automatic method are compared to a manual detection of calving events to verify and evaluate its accuracy. We use observations 15 from two time-lapse cameras placed at two different locations in front of a tidewater glacier in Svalbard, Tunabreen, during two consecutive summers, 2014 and 2015. The calving event locations at the front are compared to the position of two rising plumes observed at the front and the distribution of event-sizes is compared to a power-law curve.

Study area
Tunabreen is a 27 km long surging tidewater glacier, draining from the Lomonosovfonna ice cap terminating at the head of 20 Templefjorden in central Svalbard (see Fig. 1). Its drainage area is approximately 174 km 2 (Nuth et al., 2013). This glacier is known to have surged in 1930, 1970 and lately in 2002-2005, experiencing multiple retreats and advances and leaving submarine footprints (Forwick et al., 2010;Flink et al., 2015). It has retreated from its maximum extent, in 2004, until 2016 when it started to surge again. The 3 km wide terminus is roughly 70 m thick and grounded in 40 m deep water (Flink et al., 2015). At the front of Tunabreen, there is one main subglacial drainage portal (see Fig. 1) that can also be seen on the pictures 25 (see Fig. 10 and Fig. 1b). The main subglacial plume, plume 1, is described in detail in Schild et al. (2017). A second subglacial plume, plume 2, less visible, is intermittently present on the pictures (see Fig. 1b). There are also two terrestrial melt-streams, one on each side of the glacier.

Time-lapse cameras
In 2014, the time-lapse camera was installed in front of Tunabreen calving front, in the moraine field, at the coordinates (N78 • 27.084' E17 • 16.195' 73 m) as shown in Fig. 1. The camera (Canon EOS 450D) was placed on a tripod in a waterproof plastic box (Pelican Storm IM2075) with a drilled hole in front of the camera (see Fig. 1). The intervalometer is a Digisnap 5 2700 from Harbortronics with a low temperature modification. The system was powered by a 12 V alkaline battery pack, placed in a plastic box covered by stones. The time interval between pictures was 14 min.
In 2015, an additional time-lapse camera was installed on Ultunafjella, a rock outcrop to the west of the glacier terminus.
The camera (Canon EOS 700D) was enclosed in a custom-made peli-case box along with a Harbortronics Digisnap 2700 intervalometer. This system was powered by an external 12 V battery and a 10W solar panel. The camera system was installed 10 on a tripod, and the tripod was buried into the ground and anchored by rocks.
Both years, the cameras were set in aperture priority (set aperture value and automatic shutter speed). Camera properties are presented in Table 1.

Image registration
Detecting changes in a sequence of images requires perfect alignment of the images so that non-moving objects are at the same location for every image in the sequence. However, that condition may not hold in our case. Due to weather conditions, the camera can slightly move and rotate. This may cause false change detection if the images are not geometrically aligned. Thus, image registration is an important step in calving event detection.  The first step in automatic detection of calving events is to isolate the calving front from the surroundings to avoid detection of changes at the surface of the glacier or the ocean. First, the image is cropped around the glacier front geometry so that most noise is removed. The cropping region is determined from the first image. Second, we use a region-based active contour 25 method, the Chan-Vese model (Chan and Vese, 2001), to isolate the front and the ocean.
The Chan-Vese model uses a region-based active contour, which works based on region information instead of edge information (boundary detection) as used by an edge-based active contour. We adopt the Chan-Vese model because it does not rely heavily on edge detection, which is often hard to find between the surface and the front of the glacier. It uses an evolving curve to detect objects in a given image u 0 . To do that, the Chan-Vese segmentation minimises an energy function defined as where µ, ν, λ 1 and λ 2 are non-negative parameters, and C is a curve that separates the image into two regions. c 1 and c 2 denote mean intensity values of the two regions, inside C and outside C, respectively. Eq. 1 reaches its minimum if and only if the 5 curve C lies at the boundary of two homogeneous areas.
At the first iteration, the user is required to delineate two initial masks, around the ocean, M ocean process, if the visibility is bad we use the mask from the previous iteration.

Change detection
The goal is to detect changes between two grayscale images masked by the intersection of The rotational dependent Local Binary Pattern (LBP) is a simple visual descriptor method to extract the image texture T i of image i (Ojala et al., 1996(Ojala et al., , 2002Pietikäinen et al., 2011). The center pixel, (x c , y c ), in grayscale value g c , is compared to its P neighbours values situated at a radius R of coordinates {(x c + R cos(2πp/P ), y c + R sin(2πp/P )), p ∈ [0, P − 1]}. The grayscale value, g p , is linearly interpolated if not falling at the center of a pixel. The texture value of the center pixel (x c , y c ) is 20 then defined as with δ(x) the Heaviside step function, which gives 0 if x < 0 and 1 if x ≥ 0. Here we use P = 20 and R = 5 for 2014 and R = 6 for 2015. The different values depend on the characteristics of the camera (focal length, distance to the front, resolution, etc.) and the choice is empirical. Parameters P and R have been sampled in [8,10,15,20,30,40] and [5,10,15] respectively. 25 The quality of the settings has been assessed by comparing errors between the automatic and the manual methods using the comparison metric presented in subsection 3.5. For both year, best results were achieved with P = 20 disregarding the value of R. The value of the radius R is sensitive to the size of the image pixel.
To determine the amount of change in the texture of two consecutive images, we use the Combined Difference Image and K-means clustering (CDI-K) as proposed by Zheng et al. (2014). Two operators, the absolute difference and the logarithmic difference, are applied to the texture images T i−1 and T i , and produce two change-images, D s and D t respectively The change-images are normalised to the range [0, 255]. We then apply a mean filter (11×11) to D s in order to remove isolated pixels and a median filter (3×3) to D t in order to remove isolated pixels but preserve edges. The resulting images, D s and D t are combined, the former achieving sleekness and the latter maintaining edge information, using a weight parameter w > 0 We set the weight parameter w = 0.1 as suggested by Zheng et al. (2014).
To decide if there is change or no-change, we choose to use a local neighbour threshold because of its deterministic nature in contrast with the randomness of K-means. A median function is applied to D on a 25×25 pixel window.
When the weather is too foggy in one image, the whole calving front is detected as calved and it takes a long time for the 15 algorithm to perform. We therefore remove images that are falling in this category by calculating the coefficient of variation (CV) of the calving front, which gives an idea of the intensity distribution of the image. If the CV is above or below a certain threshold, determined beforehand from a set of images, the image is removed. The threshold has been determined by comparing errors between manual and automatic methods using the comparison metric presented in subsection 3.5. Above the threshold, the image is highly illuminated and below the threshold, the image is covered by dense fog. The front is almost not recognisable 20 so calculation is impossible and errors are systematically very large and the size of the automatically detected calving events is close to the size of the front.

Change mask reconstruction
The end product of the change detection method detailed above is generally a set of clustered points of change (see Fig. 2a) and we use a mask reconstruction method, called the α-shape method (Edelsbrunner et al., 1983), to transform the result into 25 polygons (see Fig. 2b). We can determine the set of points situated on the boundary of any empty open disk of radius α (the points situated on the blue curve in Fig. 3). The polygon is then constructed from these points (in red in Fig. 3). α is a very sensitive parameter and has to be chosen with care. If α tends towards zero, only the points themselves are kept. The smaller is α, the smaller are the empty open disks around the points and the smaller is the α-shape (the reconstructed area). If α tends towards infinity, the α-shape is the convex-hull of the points. We further test the sensitivity of the mask reconstruction 30 parameter α. A polygon is kept if the number of pixels it contains is higher than a threshold otherwise it is considered as noise. The size threshold depends on the real size of the pixels (described in the next section). To simplify, we use the real size of a pixel averaged over the whole front and a threshold of 50 m 2 .

Size calculation
To estimate the pixel size we use the photogrammetric concept that the real pixel size is dependent on the focal length of 5 the camera, the distance from camera to the front and the physical size of the pixel at the sensor. The distance to the front is calculated from georeferenced satellite images Landsat 8 OLI/TIRS C1 Level-1 downloaded from USGS EROS (ten images from July-September 2014 and one image from August 2015).The real pixel area of the image pixel ( where d i X and d i Z are the real lengths of the pixel in the horizontal and vertical directions of the image, respectively (see Fig. 4).  The front is outlined with a user control tool in the image so we can determine the distance of the front in pixels, d. We outline the front on the satellite image and project it into the camera horizontal direction (black line in Fig. 4a) to get a number of d evenly distributed points. When projected back to the front (red line in Fig. 4a), the distance between each point (X i , Y i ) is the real pixel size in the horizontal direction. To perform the projection we use reference points on the front, camera position (X cam , Y cam ) and the focal angle θ = 2 arctan(x s /(2l f )) where x s is the sensor size in the horizontal direction and l f is the 5 focal length (see Table 1). In the vertical direction, we assume a vertical front and calculate the distance, L i , from the camera to the front. We assume no distortion of the projected pixels in the vertical direction so that each pixel on a vertical column has the same size (same distance from the camera) and we correct for the elevation, Z cam , of the camera location. The real size of the pixel in the vertical direction, d i Z , is then 10 where L i = ((X cam − X i ) 2 + (Y cam − Y i ) 2 + Z 2 cam ) 1/2 is the distance from the camera to the front at the sea-level and y p x is the real size of a pixel in the y-direction. Lens distortions are neglected here. Given the radial distortion coefficients of the lens, the focal length and the pixel size, the distortion error for the front pixel situated the farthest from the camera was smaller than 0.2 % for both lenses. 15 In order to assess the performance of the automatic detection method, we compare the results to a manual method on the same set of images using human visual detection. This does not give the absolute accuracy of the method in relation to true calving since this would require an independent dataset. In 2014, 1100 images were used for comparison of the automatic detection for the period 26 August to 6 September. In 2015, 469 images were used for the period 1 September to 4 September. which represents a calving event, is outlined manually and the coordinates are saved. The manual detection for 2014 has been performed by Westrin (2015).

Comparison metrics
The location and size of each calving event are the main results we want to retrieve. For the location, we construct a confusion matrix, which gives the agreements and disagreements between the automatic and the visual detection for each pair of masked 5 images (e.g. Stehman, 1997) and which contains:  The number of pixels automatically detected as calved and non-calved are P auto = T P + F P and N auto = T N + F N respectively. The number of pixels visually detected as calved and non-calved are P visual = T P + F N and N visual = T N + F P .
The total number of pixels is m = T P +F P +T N +F N . It is important to note that, for the case of glacier such as Tunabreen with limited calving event size, the number of pixels labeled as calved is generally small compared to non-calved pixels. In this sense, the accuracy measure Acc = (T P + T N )/m is not appropriate since T P << T N as stated by Kubat et al. (1998) 15 for imbalance classes. Here we need a method to detect changed regions, instead of unchanged ones. They recommend the use of F-measure to measure a classification problem with imbalanced classes. However, F-measure has a property that is invariant under the change of TN (Kubat et al., 1998), thus it is less sensitive to changes in unchanged pixels. Alternatively, the Matthews correlation coefficient (Matthews, 1975), 20 is a better measure for imbalanced classes and it takes into account both true and false positive. Its interpretation is similar to the Pearson correlation coefficient between the observed and predicted binary classifications. A value of +1 represents perfect match while -1 total disagreement. Perfect match is only happening when no change is detected by both methods. When there is no overlapping or a calving event is only detected by one of the methods, the M cc will give a low value by definition whatever the size of the detected event. To be able to compare the magnitude of the error, we can look at the positive difference as the 25 difference between detected calved zones from both methods divided by the total number of pixels, to assess the importance of the mismatch.
The weather conditions and illumination change the pixel intensities and can inhibit detection in some cases. We determined five categories: Normal conditions (N), light fog (LF), dense fog (DF), high illumination (HI) and low illumination (LI). Some examples are shown in Fig. 5. When comparing manual and automatic detection, we visually determine the categories qualitatively. Each image has also been cut in six different sections placed in a weather condition category.

Comparison assessment
An example of calving detection from the automatic and manual methods is shown in Fig. 6.  5 (2015) and can be considered as good matches. When M cc = 0, the difference between pixels, P dif f is close to zero. On the contrary, dense fog (DF), high illumination (HI) conditions and any combination with one or the other 5 show a relatively low M cc and high P dif f particularly. If such a configuration is combined with any of the others, the result is also poor.
We decided to remove all pictures with high illumination and dense fog from the analysis. To determine whether a picture is falling in this category, we look at the standard deviation and mean of pixel intensity for each section of the image after calibrating these values for normal conditions on a set of images. Pictures with high mean intensity compared to normal  range for each section fall into high illumination category. Pictures with low CV (ratio between standard deviation and mean intensity) fall into dense fog category. Moreover, given the calving characteristics of Tunabreen (generally no glacier-wide calving events), we assume that if the detected calving size is greater than a certain value depending on the front size, the detection is not satisfactory and we thus remove it from the results. Fig. 8 shows the results for 2014 and 2015. This method has limitations since even if it performs well at removing bad detection, it also removes good ones.

5
To test the sensitivity of the mask reconstruction parameter α, we look at five different values and compare the results after removing unfavorable weather. This is shown in Fig. 9. If α is too small, the event sizes from the automatic detection tend to be smaller than from the manual detection and vice-versa for too big α. In Fig. 9, the best fit between manual and automatic detected calving events is when α = 10 given the M cc and P dif f results.  10 The initial number of pictures to be analysed in 2014 was 6292 but because of weather conditions and camera settings, only 3497 usable images remain, of which 2084 showed no calving. In total, 2575 calving events have been detected ranging from 20 to 3500 m 2 . An example of mask delineation and calving detection is given in Fig. 10a. To facilitate the analysis, the front is divided in six sections as shown in Fig. 10a In 2015 3242 images were analysed, but 571 of these were removed as the glacier terminus was obscured by water droplets on the porthole cover of the time-lapse enclosure. Because of this, there is no coverage between 23-27 August. After weather filtering, 1495 images are left including 898 with no calving detection. In total, 1647 events were detected (222 in section 1, 327 in section 2, 228 in section 3, 234 in section 4, 398 in section 5 and 238 in section 6). An example of mask delineation and calving detection is given in Fig. 10b.

Spatio-temporal size distribution
In Fig. 11, we present the total area detected per image (in blue) and the daily moving average (in red). It is difficult to establish seasonal pattern because of the data gaps during unfavorable weather conditions.
14 Geosci. Instrum. Method. Data Syst. Discuss., https://doi.org /10.5194/gi-2018-5 Manuscript under review for journal Geosci. Instrum. Method. Data Syst. Instead, we explore the relative size distribution for the different sections to estimate which section undergoes the largest calving events per image (see Fig. 12a-c). When calving events are detected on an image, we look at the size proportion of each section. If, for a pair of images, calving is detected in only one section, this section gets 100% of the total. If other events occur within the same time frame in other sections, these sections get a proportional fraction corresponding to their size compared with the total calved size of the whole calving front. The primary event section receives the largest share. In Fig. 12b-d, for 5 each primary event section, we show the share of the secondary event.
If a section gets 100% of the total size during a time frame, it is more often section 2, followed by section 5, in 2014, whereas, in 2015, it is more often section 5 followed by section 2. In general, section 2 for 2014 and section 5 for 2015 have the most occurrences of primary calving event size. Also, when primary calving event size occurs in section 5, secondary calving occurs particularly in section 2. For 2015, section 5 is also secondary when section 2 is primary. It is also where the most margin 10 retreat occurs in the season. Furthermore, there seems to be a link between sections. When a primary calving event occurs in a section, secondary calving event sizes are usually found in adjacent sections.  The spatio-temporal information such a method provides has the potential to enhance understanding of certain processes or validate calving theories. Here we used two time series of two different years from two different locations and camera settings.

Event-size distribution
Some problems still need to be addressed and are discussed in the following section. The long time-lapse and the weather-5 related issues do not allow direct conclusion on the spatio-temporal calving behaviour of Tunabreen. Still, it appears that most of the large calving events happened near plume locations (sections 2 and 5) with larger events at the most active plume location, plume 1 (section 2). The plume in section 5 has not been detected elsewhere in the literature so this conclusion has to be taken with care and this should be investigated further. It also appears that when a large event is triggered at the less active plume location, plume 2 (section 5), a smaller one is simultaneously triggered at of plume 1. Moreover, events in one section 10 seem to trigger smaller events in adjacent sections, which confirms the destabilisation of the local neighborhood of the calving region observed by Chapuis and Tetzlaff (2014).
A self-organised critical (SOC) system in nature, is a system that exhibits a slow and steady accumulation of an instability followed by a rapid relaxation triggered from a single point and, independently to external conditions, leading to the collapse of any possible size (Jensen, 1998). Åström et al. (2014) showed that the calving front behaviour of tidewater glaciers had the 15 characteristics of an Abelian sandpile model, a simple SOC system (Dhar, 1999). The probability distribution at the critical point approaches a power-law with an exponent close to −1.2, which is also the case for the results showed in Fig. 13. Cut-offs for large sizes have also been reported in Åström et al. (2014) for smaller glaciers and they are due to the glacier front size. For smaller sizes, the relative abundance is below the power-law estimate. This is due either to the long time-lapse between images that does not allow small events to be detected, or to the spatial resolution of the camera. The system always tends towards criticality without regard to external factors. However, external conditions have the potential to change the system from sub-to supercritical. Investigation of relationship between calving events and environmental controls will be the subject of a separate paper.

Errors and limitations 5
The main limitation is clearly the poor performances due to the weather: in presence of fog or when the calving terminus is strongly illuminated by the sun. The illumination problem is worst during daytime when the sun is reflected on the ice. At those times, mostly in the middle of the day, the temperature is high leading to surface melt, subglacial discharge and active plume.
The marine melt feedback is therefore difficult to analyse with this method if the weather condition problem is not solved.
Nevertheless, in normal conditions or low illumination, the method performs well compared with the manual method and is a 10 good tool to assess the location of the calving events. Errors in detection, however, may arise because of different factors. First, the method requires different parameters at each step and most of them are defined empirically (trial and error), particularly the change detection and mask recognition parameters but also the threshold for unfavorable weather condition and minimum/maximum calving size. Second, the method detects frontal changes and is not able to recognise external objects coming into view of the camera, such as birds or icebergs, nor make the difference between frontal ice changing position, because of 15 gravity for example, and actual calving. Third, while the glacier front calves, it also advances due to ice flow and both processes have an effect on the front position. Because of that, the front outline needs to be manually refined periodically, depending on the flow speed of the glacier and the camera view angle. For Tunabreen, refinement was needed every 1000 pictures or every 8-10 days. Fourth, no subaqueous calving can be detected with this method although, for certain glaciers, this could account for a large proportion of the frontal mass loss. Finally, the manual detection is subject to error because of a failure to correctly 20 identify or digitise events.

Recommendations and future improvements
Despite inherent limitations, the method can still be improved. First, the camera settings and camera installation have an impact on the result and can be optimised. For the setup in 2014, the large aperture setting (f/2.8) and the lack of polarised filter rendered more pictures compared to 2015 to fall into the HI category. It is important to consider all situations the automatic 25 settings of the camera will have to handle. Also, instead of fixing the aperture width, one could fix the aperture speed to let the lens adapt to the amount of light available or adapt the ISO, as in Kwasnitschka et al. (2016).The position of the camera also influences the results and our two different camera locations show both benefits and drawbacks. A camera positioned at low elevation and in front of the glacier (2014) is favourable for the segmentation (refinement has to be done less often) whereas a camera situated more to the side and at higher elevation (2015) is favourable to recognise larger events but parts of the front 30 may be occluded. In some images, calving events are recorded at different locations and having a shorter time-lapse interval (to be set according to the characteristics of the glacier) could help to distinguish single events. Storing the images in lossy file can be changed or if there is a service of automatic downloading of images, it would help to store the pictures in raw format so that more postprocessing is available.
Second, the method can be improved in detecting or enhancing images with unfavorable weather. Fog detection is implemented for different applications (e.g. automatic car navigation) and could be useful (e.g. Hautière et al., 2014). If weather observations are available, it can also be good to combine these data to detect total whiteout of an area for example (Jiskoot 5 et al., 2015).
Third, the size calculation only gives the 2-dimensional area of the calving event and a third dimension could be calibrated from observation (Chapuis and Tetzlaff, 2014) or from other cameras placed at different locations. Moreover, the distance between the camera and the front is estimated from satellite data with rather low spatio-temporal resolution and data with higher resolution could improve the size calculation (e.g. satellite, other time-lapse cameras or unmanned aerial vehicle photogram-10 metry).
This automatic method has been developed for a particular glacier, Tunabreen, but can be calibrated for other glaciers, even if the spatio-temporal scale is different. For instance, glaciers in Greenland, where calved icebergs are larger and detached more often, are a good candidate. Time-lapse cameras are easy to place in front of a glacier to monitor the calving state of it and this method can help to automatise an otherwise difficult task.