Surveying Drifting Icebergs and Ice Islands: Deterioration Detection and Mass Estimation with Aerial Photogrammetry and Laser Scanning

: Icebergs and ice islands (large, tabular icebergs) are challenging targets to survey due to their size, mobility, remote locations, and potentially difﬁcult environmental conditions. Here, we assess the precision and utility of aerial photography surveying with structure-from-motion multi-view stereo photogrammetry processing (SfM) and vessel-based terrestrial laser scanning (TLS) for iceberg deterioration detection and mass estimation. For both techniques, we determine the minimum amount of change required to reliably resolve iceberg deterioration, the deterioration detection threshold (DDT), using triplicate surveys of two iceberg survey targets. We also calculate their relative uncertainties for iceberg mass estimation. The quality of deployed Global Positioning System (GPS) units that were used for drift correction and scale assignment was a major determinant of point cloud precision. When dual-frequency GPS receivers were deployed, DDT values of 2.5 and 0.40 m were calculated for the TLS and SfM point clouds, respectively. In contrast, values of 6.6 and 3.4 m were calculated when tracking beacons with lower-quality GPS were used. The SfM dataset was also more precise when used for iceberg mass estimation, and we recommend further development of this technique for iceberg-related end-uses.


Introduction
Recent calving events from ice shelves and floating glacier tongues have created immense tabular icebergs in both the Arctic (e.g., Petermann Glacier:~300 km 2 ) [1] and Antarctic (e.g., Larsen-C Ice Shelf: 5800 km 2 ) [2]. Mass and dimensional estimates of these 'ice islands' are used by the offshore resource extraction and shipping industries to model the risk associated with these ice hazards, as well as that imposed by the smaller ice island fragments and icebergs produced during their deterioration [3][4][5][6]. Deterioration model output is used in operational drift models that predict the trajectory and the mass of ice hazards through time and space [6,7]. Iceberg deterioration is also modeled to improve our understanding of the biological and physical influence of meltwater input to the surrounding ocean waters [8][9][10][11][12][13]. Unfortunately, the calibration and validation of both iceberg and ice island deterioration models has been restrained by a paucity of in-situ data sets that accurately capture both the morphological change and corresponding environmental data required [14,15]. This is largely due to logistical, financial, and safety issues that accompany surveying these subjects in remote locations [16,17]. In addition, one must overcome the technical challenges associated with surveying a The two surveying methods are evaluated with consecutive, triplicate surveys over which time (<70 min) we assume that no appreciable deterioration (i.e., from melting or calving) has occurred. Our survey subjects include two freely drifting tabular icebergs, which were surveyed near Newfoundland and Labrador, Canada during a field campaign in April 2015. These icebergs were likely fragments of larger ice islands originating from the Petermann Glacier in 2012 (see Section 3). Regardless, our methods and results are transferable across the size range of small icebergs to large ice islands, as well as various iceberg morphologies.
This article uses numerous abbreviations, especially in the methods sections. Please note that a list of abbreviations is provided at the end of the article.

Field Setting
A 130 km 2 ice island calved from northwest Greenland's Petermann Glacier in 2012 and generated numerous ice islands, icebergs, and smaller fragments throughout its deterioration [1]. Figure 1 shows the location and field photos of the two freely drifting tabular icebergs that we surveyed. The icebergs ('i1' and 'i2'),were likely fragments of the 2012 Petermann Ice Island. The surface ridging of i1 was similar to previously observed ice islands that originated from the Petermann Glacier. i1 and i2 were surveyed in triplicate from the CCGS Amundsen on 23 and 27 April 2015, respectively. i1 was situated at 51 •

Marker/GPS Deployment
Markers, paired with either single-frequency GPS tracking beacons (data telemetered) or dual-frequency GPS receivers (data logged internally), were deployed on i1 and i2 via helicopter prior to aerial photography and TLS surveying. The markers, shown in Figure 2, were composed of 0.4 m 2 black burlap squares and were deployed to serve as recognizable features for the SfM image processing software during photo alignment. The GPS units were identified in multiple images to increase the quality of photo alignment. In the absence of a definitive reference survey, GPS data were used to assign each cloud point to a position in an 'iceberg reference frame' (IBRF). This is defined by the relative positions of cloud points that were collected when surveying the rotating and translating iceberg target [37,38]. The IBRF was generated with GPS data that were also used to provide scale to SfM-generated point clouds (Section 3.1.1) and to correct the TLS point clouds for iceberg drift during surveying (Section 3.1.2). The deployed markers were also used in a calculation of inter-survey error (ISE) after the individual 3D point clouds were generated and scaled (Section 3.1.1).
The combination of each survey technique (SfM or TLS) and deployed GPS units (tracking beacons on i1 or high-quality, dual-frequency GPS receivers on i2) are the four survey scenarios that we evaluated. Table 1 contains the quoted accuracies for each type of deployed GPS unit. Examples of marker-GPS pairs that were used during SfM processing are shown in Figure 2a,b. Three markers, paired with either a Canatec tracking beacon (Canatec and Associates International Ltd., Calgary, AB, Canada) (n = 1) or a Solara FieldTracker (FT) 2000 beacon (Solara Remote Data Delivery Inc., Winnipeg, MB, Canada) (n = 2), were deployed on i1. Positions were measured at 5 min and 30 s intervals, respectively. Three markers were deployed on i2 and were paired with either a Topcon Hiper V (Topcon Corp., Tokyo, Japan) (n = 2) or a Trimble 5800 (Trimble Inc., Sunnyvale, CA, USA) (n = 1) dual-frequency GPS receiver, which recorded at a 5 s interval. Data logged by the dual-frequency receivers were post-corrected in kinematic mode using precise point positioning (PPP) processing of global navigation satellite system data [45], which also yielded uncertainty information regarding positional accuracy (σ bx , σ by , and σ bz ) ( Table 1). Cantec and Associates reported their beacon accuracy as a 50% circular error probable value [46], which was converted to standard deviation units (σ bx and σ by ), while Solara beacon positional uncertainty was derived from data collected by four stationary units over a 36 h period (Table 1). Table 1. Quoted accuracies and calculated positional uncertainties used in data processing for Canatec beacons units [47], Solara FT2000 tracking beacons [48] and higher-quality, dual-frequency GPS following Precise Point Positioning (PPP) processing using Natural Resource Canada's online tool [45,49].  Figure 2c). A wireless intervalometer (CamRanger, Pinedale, WY, USA) was used to initiate a series of image acquisitions for each flight line by an operator within the helicopter. Images were acquired at a 1 s interval and georeferenced with a Nikon GP-1 GPS accessory offering a 10 m root-mean-squared horizontal accuracy and an update rate of 1 s. Flight lines were planned with a nominal flight speed (20 to 40 m s −1 ) and altitude (280 to 300 m) to ensure~80% overlap and~40% side-lap between images. This resulted in ground sampling distance between 2.4 and 3.4 cm. Three consecutive surveys were completed for i1 over 26 min; however, only two aerial photography surveys were conducted over 14 min for the smaller i2, as inclement weather forced the helicopter to land. All surveys pertaining to either i1 or i2, conducted with either aerial photography or TLS, were completed in a short timeframe over which no appreciable deterioration was expected to occur. Details on the instruments and field execution for both aerial photography and TLS surveying are contained in Table 2.

Terrestrial Laser Scan Surveying
For TLS surveying, an Optech Ilris-HD ('Intelligent Laser Ranging Imaging System') laser scanner (Optech Inc., Toronto, ON, Canada), which used a 1535 nm wavelength, had a 10,000 Hz pulse repetition rate and a 40 • field of view. It was integrated with an INS consisting of an iIMU-FSAS (iMAR GmbH, St. Ingbert, Germany) motion sensor and a synchronous position, attitude, and navigation (SPAN) system (NovAtel Inc., Calgary, Canada). The latter included two NovAtel GPS antennas for heading measurement ( Table 2). The TLS and INS were situated on a custom-built telescoping steel frame, mounted on the ship's barge, with a 3.3 m baseline between the two global navigation satellite system (GNSS) antennas (Figure 2d). High sea ice concentrations impeded the deployment of the barge on which the survey instruments were mounted. Data were therefore collected with the instruments as originally positioned but while the barge was stowed on the CCGS Amundsen's deck. The laser scanner was situated at an elevation of~15 m above the sea surface at this location. INS and TLS data were acquired and integrated by Hypack ® (Version 2015; Xylem Inc., Middletown, NY, USA) software and the position of the TLS survey points was recorded in Universal Transverse Mercator (UTM) projection. A PPP solution, found with the POSPac MMS (Version 7.1; Applanix, Richmond Hill, ON, Canada) software, improved the accuracy of the points at a later date.
Both icebergs were surveyed by TLS over three consecutive circumnavigations by the CCGS Amundsen. Surveying was conducted over 68 and 45 min for i1 and i2, respectively. The distance between the ship and the icebergs ranged between 20 and 110 m for i1, and between 60 and 140 m for i2. A geolocation error was found in the i2 INS dataset, which likely stemmed from a poor sky view and multi-pathing of the satellite signal associated with GPS configuration on the CCGS Amundsen. This was corrected by using a time-based offset between the ship's C-Nav3050 (Oceaneering International, Inc., Houston, USA) differential GNSS receiver and the survey INS position. All survey data and processing scripts can be found under Polar Data Catalogue entry 12818. See the supplementary materials section for further information.

Structure-from-Motion Photogrammetry
Photoscan Pro (V. CC 2014) (Agisoft LLC, St. Petersburg, Russia) was used to derive 3D models from aerial photography surveys with SfM processing methods. The output is referred to here as 'SfM point clouds'. This software contains algorithms for feature matching and point cloud generation, and provides post-processing tools for scale transformation and coordinate assignment in the IBRF [29,30]. Photos that captured a portion of the survey target were imported into the software and aligned with the highest accuracy alignment setting. Photo resolution was upscaled by a factor of four with bilinear interpolation to aid the recognition of surface features (Agisoft, personal communication). During the alignment step, camera positions and orientations were also calculated and a sparse point cloud was generated through iterative bundle adjustment. This adjustment iteratively optimizes scene coordinates and camera locations to minimize reprojection errors [21,50,51]. Since photo alignment does not rely on georeferenced photos, there was no need to correct for drift [32]. Point cloud distortion was decreased with an optimization procedure, which improves the estimation of point coordinates and camera orientations based on reference coordinates assigned to the deployed GPS units [50]. To do this, the GPS units were identified in each image, and a dense point cloud was then constructed with the calculated camera orientations at the highest quality setting. These basic processing steps, as well as those involved with point cloud comparison, are illustrated in Figure 3.
Point clouds generated from SfM processing were scaled by assigning a distance value between the identified GPS units. The mean Euclidean distance between the three GPS units (the 'inter-GPS distances') was calculated from the positions logged over the total duration of a survey type (i.e., all aerial photography surveys of i1). The GPS units deployed on i1 did not record at a constant interval, so a cubic spline interpolation at 1 s intervals between transmissions was conducted to calculate 2D inter-GPS distances between corresponding timestamps. 3D inter-GPS distances were calculated from the positions logged at 5 s intervals by the GPS units deployed on i2 since elevation data were provided by the dual-frequency GPS receivers. Points that represented the surrounding ocean and sea ice were removed manually following a visual assessment. The first of the three surveys was designated to be a reference point cloud for the purpose of survey comparisons. This point cloud was cropped to a wider extent than subsequent survey point clouds to improve the chance that all points in the compared clouds would have corresponding points to be compared to on the reference cloud (Section 4.3).
A point cloud registration error, representing the uncertainty of the positions assigned to a 3D model during photo alignment and point cloud georeferencing [52], is usually calculated with respect to a definitive reference survey (c.f., Westoby et al. [53]). This error could be introduced into the workflow by incorrectly identifying the center of GPS units in imagery or photo-alignment issues in AgiSoft Photoscan (Section 4.2.1). A definitive reference survey to which we could compare our 3D iceberg models was not available. Therefore, we used the marker centers as common points in each survey [54] and derived a proxy registration error which we refer to as the inter-survey error (ISE). We examined orthophotographs (2.8 cm pixel spacing) and DEMs (5.7 cm pixel spacing) after image alignment and point cloud scaling to extract the x, y (from orthophotographs), and z (from DEM) position of the three marker centers for each survey. We then calculated the difference in the orthophoto/DEM-derived marker position between the three surveys and took the mean of these differences (n = 9; three marker positions compared across three surveys) to represent ISE. We assume these errors to be isotropic and uniform across the survey target's surface (c.f., Westoby et al. [53]).

Terrestrial Laser Scanning
TLS surveys were corrected for iceberg drift with custom R (V. 3.0.2) and Matlab (V. R2015A) scripts. This procedure involved determining distances and offset angles between cloud points and deployed GPS units by matching timestamps and then re-projecting the cloud points in an IBRF. Canatec and Solara beacon data were first splined with a natural cubic curve to constant time intervals due to slight logging interval variations. Positions from two GPS units per iceberg were Kalman filtered (forward-stepping) and Kalman smoothed (backward-stepping) to improve estimates of these GPS locations. An isometric Kalman gain was applied and visually tuned for each beacon type. The Kalman filtered data was passed through a natural cubic spline with a 0.1 s interval output to improve timestamp matching between the GPS and cloud points. Figure 4 provides a simple schematic of the offset angles and distances, calculated below, that were used in the drift correction procedure. Points were matched with the GPS position recorded with the closest timestep. The distance (d) between the positions of one GPS unit and each cloud point was calculated with Equation (1), while Equations (2) and (3) were used to calculate the heading between the cloud point and a second GPS unit (h pc ). It is irrelevant which GPS units are used for the offset calculations so long as they are used consistently in their respective roles. The differences in x and y Cartesian coordinates are denoted with dx and dy, respectively.
(2) c is equal to 1 or −1 based on the sign of the corresponding dx of the survey point, with a value of 1 assigned if dx = 0. Finally, if h pc < 0: A heading in degrees was found by multiplying h pc by 180/π. The offset (h OFF ) between h pc and the heading between the two GPS units (h GPS ), and calculated d, were used in Equations (4) and (5) to correct cloud points to x', y' positions in an IBRF defined with an origin at x o , y o and a heading of 0 • . The assigned values of x o , y o , and the initial heading are arbitrary as the distance and angle offsets to each cloud point will remain consistent.

Point Cloud Comparison and DDT Calculation
The corrected point clouds, one corresponding to each SfM point cloud or TLS survey, were opened in CloudCompare (V2) where background sea ice returns were manually removed from the TLS point clouds.
CloudCompare is an open-source point cloud comparison software program that has been successfully used to detect landform changes [43,53,55,56], and was the software program used in this study to compare point clouds representing the individual surveys of each tabular iceberg.

The Multiscale Model-to-Model Cloud Comparison Algorithm
The precision of each surveying technique was determined by calculating the linear distance between the drift-corrected point clouds in CloudCompare with the Multiscale Model-to-Model Cloud Comparison (M3C2) algorithm [43,53]. By eliminating the need to grid or mesh by directly comparing point clouds, no spurious distance values are returned due to varying point density (ρ p ), data gaps and noise [43,57].
The M3C2 algorithm creates cylindrical volumes around a subset of core points that are oriented normal to the surface of the reference point cloud. The reference point cloud was subset by a factor of 5 to 10 to generate a group of core points to decrease processing time, as suggested by Lague et al. [43]. The distances between the core points in the reference cloud and analogous points in the compared point cloud (PC dist ) were then found ( Figure 5). The direction of the normal vectors were determined based on a plane that was fit to all points within a user-defined diameter, D, from each core point ( Figure 5). D was set to be sufficiently large (>20-25 times the roughness of the point cloud surface) to ensure that small-scale surface roughness variations did not bias the calculation of normal vectors [43,53]. The normal scale assessment criterion, ξ, was used to ensure that D was adequately assigned: where σ 1i (D) is the surface roughness of the reference cloud, calculated as the standard deviation of each point to the plane found at scale D centered around each core point (i) [43,53]. We aimed to assign a value of D that satisfied the condition ξ i > 20 (equivalent to a normal orientation error of < 2%) for 98% of core points, as did Lague et al. [43] and Westoby et al. [53]. The projection scale (p) and maximum search distance (δ) defined the diameter and length (respectively) of the cylinder that was projected in the direction of the surface normal at each core point ( Figure 5) [43,53]. δ was set to decrease computation time and can be thought of as an outlier removal as points beyond this distance will not be included in PC dist or u calculations. The value of p was set to ensure that an average of five or more points were within the cylindrical projection space based on the ρ p of the reference cloud [43]. The locations of all cloud points in the normal direction and contained within the cylinder were averaged for both the reference and compared point clouds by the M3C2 algorithm and were used to calculate the PC dist between the two point clouds at each core point [43,58]. In step 1, the surface normal (N) is determined by fitting a plane to the points (open circles) within a user-defined diameter, D, from each core point (grey circle). In step 2, a cylinder with diameter p is centered on the core point and projected along the direction represented by N. The algorithm projects the cylinder both above and below the core point, but the figure shows the projection in a single direction for simplification. The cylinder length is the maximum search distance (δ) and this defines a space where matching points in both point clouds are averaged and differenced to calculate the point cloud distance (PC dist ). The open circles represent cloud points. The grey circles represent core points. Adapted from Lague et al. [43].
In addition to computing PC dist , M3C2 output includes a 95% confidence interval (or uncertainty; u) for each PC dist value. This considers the fact that two surveys will always have slightly different surface roughness at any given location, which could cause spurious distances to be calculated. Equation (7) considers this in the calculation of u with the inclusion of surface roughness (σ 1i and σ 2i ) of the reference and compared point cloud as well as the ρ p (ρ 1p and ρ 2p ) within the projected cylinder [43]: where n = number of points within the projected cylinder. The ISE (traditionally, the registration error) value is an optional parameter which is added to the confidence interval in the calculation of u i . The inclusion of point cloud alignment error along with surface roughness in this uncertainty calculation is a strength of the M3C2 algorithm [43,56], which resulted in a more conservative uncertainty estimation and final DDT calculation. Our study uses u and PC dist to detect an iceberg-specific DDT. This was done through a separate bootstrapping analysis which is described in Section 3.2.4.

Structure-from-Motion Photogrammetry Derived Point Clouds
Optimal M3C2 parameters were determined using the most complete survey for both icebergs (i2, survey 1) and these settings (D = 2.3 m, p = 0.29 m, δ = 8.2 m) were used for all SfM point cloud comparisons. The value of D satisfied the normal assessment criterion conditions set above and is in line with the 2 m value used by Westoby et al. [53] for a blue ice moraine complex in Antarctica. The reference point cloud was sub-sampled at 0.15 m spacing to select core points. As described in Section 4.2.1, each SfM point cloud was scaled using common inter-GPS distances.
Any potential scaling error which may result from GPS position error is therefore identically assigned to all point clouds.

Terrestrial Laser Scanning Point Clouds
Optimal M3C2 parameters (D = 14 m, p = 1.8 m, δ = 24 m) for all TLS point cloud comparisons were determined from survey 1 of i1. The reference cloud was sampled at 0.89 m spacing for core point assignment. D = 14 m resulted in 95% of points showing ξ > 20. The value of D would need to be increased to 25 m to meet the recommended 98% threshold [43]. We decided not to increase D to 25 m because the quality of the TLS surveys was obviously poor upon initial inspection due to a combination of low ρ p and noise imposed due to physical characteristics of the ice, instrument noise, varying viewpoints while surveying, and drift correction (Table 3). It was impossible to assign an ISE value in the M3C2 comparisons since there were no known common points between TLS surveys. However, the positional accuracy of TLS point cloud data hinges on the accuracy of the GPS data used for drift correction and, in this manner, potential positional error was incorporated into the TLS comparison.

Deterioration Detection Threshold Computation
The minimum amount of iceberg deterioration between successive surveys that can be reliably detected by a given survey scenario is its DDT [43]. Due to the low number of surveys which we have to conduct these comparisons, DDT is expressed as the right-tailed 95% confidence interval of the bootstrapped statistic (ε) in Equation (8) [59,60]: where ε represents the precision corresponding with each individual core point (denoted by a subscript i). The ε i values associated with each survey scenario, found with the output of all point cloud comparisons associated with a given survey technique and GPS unit type, were amalgamated for the 1000 bootstrap sample iterations. This was conducted to determine associated DDT values.

Mass Estimation
Iceberg volume was estimated with CloudCompare's 2.5D Volume algorithm [61] for both survey techniques. Reference point clouds, which had previously been cropped to a wider extent for the M3C2 analysis, were cropped to match the other point clouds. A 0.1 × 0.1 m grid spacing was assigned, and the maximum height for points within each grid cell was used to compute the volume. Empty grid cells were filled with a linear interpolated height value. The total mass (M) of each iceberg was then calculated from the sail (i.e., above water) volume (V) for each individual survey using the buoyancy calculation from McKenna and Ralph [62]: where ρ w is water density (1025 kg m −3 ) and ρ i is ice density. The latter was assigned a value of 873 kg m −3 . This was the average density of a 1 m ice core retrieved from a grounded ice island during the field campaign. Like the two surveyed icebergs, this ice island is believed to originate from the Petermann Glacier. The value of ρ i is usually assumed in analyses of Greenland ice tongues, ice islands and icebergs and values between 900 and 920 kg m −3 have been assigned in the past [12,63]. While our assigned ρ i value is derived from a single ice island, we feel that it is best to use available in-situ data when possible. We note that our mass estimations do not account for density gradients, which may exist due to saturation below the waterline. This would result in underestimation of the iceberg mass.
The precision relative to the iceberg mass (i.e., relative uncertainty) is found with the difference in mass estimates between individual surveys. The relative uncertainty percentages are the same for both mass and volume, as the two measures are directly proportional. In addition, the relative uncertainty will not change based on the assigned ρ i . We state both mass (tonnes) and volume (m 3 ) when reporting direct magnitudes.

Iceberg Characteristics and Mass Estimation
The length (L), width (W), and 2D surface extent of each survey target was determined in ArcGIS from SfM generated geotiffs (Table 3). i1 had a surface area of 5.30 × 10 4 m 2 and a relatively level surface with a mean freeboard (height above waterline) of~12 m. The smaller iceberg, i2, had a surface area of 7.00 × 10 3 m 2 and a sloped surface from a maximum freeboard of~9 m to a minimum of 2 m. i1 had a mass of 4.2 × 10 6 tonnes and volume of 4.8 × 10 6 m 3 . i2 had a mass of 7.5 × 10 5 tonnes and volume of 8.6 × 10 5 m 3 . An average mass difference of 22% was found for the estimates generated from the i1 SfM point clouds. Mass estimates from the i2 SfM point clouds, which had the best coverage from any survey scenario, differed by only 4%. However, there were only two surveys from which comparisons could be made for this survey scenario. The TLS surveys of i1 and i2 had average mass differences of 19% and 11%, respectively.

Survey Quality and DDT Values
Illustrative examples of CloudCompare's M3C2 algorithm output for the SfM and TLS point cloud comparisons are provided in Figures 6 and 7, respectively. Table 4 contains summary statistics of both SfM and TLS generated point cloud comparisons and the distribution of ε for each survey scenario is shown in Figure 8 along with the calculated DDT values. Table 4. Precision assessment of the aerial photography with structure-from-motion photogrammetry processing (SfM) and terrestrial laser scanning (TLS) surveys for i1 and i2 by survey scenario. Each scenario is determined by the combination of the surveying technique and GPS type used. TB = tracking beacon, GPS = dual-frequency GPS receiver, PC dist = distance between point clouds, u = uncertainty, DDT = deterioration detection threshold.   . (a, b) show the distance between triplicate point clouds (PC dist-i ) for i1 (a) and i2 (b), respectively. Uncertainty (u i ) in these distances at the 95% confidence level is illustrated in panels (c, d) for i1 and i2, respectively. Histograms to the right of the colour scale bar depict the frequency of the corresponding values. The insets in panels (a, b) contain example point clouds of each iceberg before M3C2 point cloud differencing. Inset colours represent height values where red represents the higher elevations. Observed ice rubble was located in the region pictured at the top of the inset in (b). The ram (i.e., an underwater 'foot' or 'skirt' that often forms underneath the waterline) likely became exposed at the opposing side of i2. The low elevations at the bottom of the inset represent the ram. Note the variation in scale between panels and in comparison to Figure 7.  . (a, b) show the distance between triplicate point clouds (PC dist-i ) for i1 (a) and i2 (b), respectively. Gray points represent 'NA' values where corresponding points between the two clouds did not exist. The uncertainty (u i ) in these distances at the 95% confidence level is illustrated in (c, d) for i1 and i2, respectively. Histograms to the right of the colour scale bar depict the frequency of the corresponding values. Point size was augmented in comparison to Figure 6 to aid visualization. The insets in (a, b) contain example point clouds of the respective ice island before M3C2 point cloud differencing. Inset colours represent height values where red represents the higher elevations. Survey line misalignment can be seen in (a) and is an artifact of the point cloud correction procedure. Note the variation in scale between panels and in comparison to Figure 6.

Aerial Photography with SfM Processing
The i1 point cloud comparison yielded a mean absolute PC dist of 1.1 m. There was a large dispersion in the individual PC dist-i values (σ = 1.7 m) (Figures 6a and 8a), and both the mean PC dist and u were greater than those of i2 (Figures 6b and 8b; Table 4). This was mostly due to the poorer quality of the GPS units deployed on i1, which increased the value assigned as the registration error (our ISE) in the M3C2 algorithm. The influence of GPS accuracy can also be seen in panels c and d of Figure 6, which illustrate the difference in u i values between the i1 and i2 surveys. It is also apparent that there were gaps in survey data for the upper horizontal surface of i1 in spite of generally high ρ p values for both survey targets (Figure 6a; Table 3). Gaps in the coverage of i1's upper horizontal surface coincided with 'featureless' regions in the original photographs. These gaps may have been caused by shadowing, smooth surfaces resulting from meltwater freezing, or illumination angles that created 'wash-out' or glare.
The DDT values for i1 and i2 were 2.5 and 0.40 m (Table 4), respectively, which again underscores the benefit of using the higher-quality dual-frequency GPS receivers over the tracking beacons. In addition, the PC dist-i values were more evenly distributed over the surface of i2 in comparison to i1 (Figure 6a,b). However, the distribution in u i (Figure 6d) showed that the distances calculated along the sidewalls and lower elevations of i2 were relatively more error-prone than those associated with the horizontal upper surface. These locations were likely less well-covered during aerial photography surveying due to the nadir-viewing angle of the camera or were associated with ice rubble that was observed at lower elevations. With such low PC dist-i , the DDT value for i2 was predominately a result of the associated u i values.

Terrestrial Laser Scan Surveying
Over the three consecutive TLS surveys, i1 translated 1280 m and rotated 27 • , while i2 translated 820 m and rotated 9 • . This equates to a drift of~0.3 m s −1 for both icebergs. In comparison to the SfM point clouds, the TLS data had greater PC dist values. There were also relatively more regions that could not be compared across surveys due to incomplete survey coverage and lower ρ p . The average ρ p of all TLS surveys was an order of magnitude lower than that of the SfM point clouds (Table 3). This can be visually assessed by comparing the insets of Figure 7a,b to those of the SfM point clouds in Figure 6a,b. The poor coverage of the upper surface in both i1 and i2 (Figure 7) was likely due to specular reflection or weak returns given the high incidence angle, reflectivity of the ice surface and increased beam footprint (c.f., Landy et al. [64], Soudarissanane et al. [65]). These factors resulted in incomplete comparisons between surveys and larger DDT values relative to SfM point clouds. DDT values for i1 and i2 were 6.6 and 3.4 m, respectively (Figure 8c,d; Table 4).

Beacon Data
GPS accuracy determined the variability of inter-GPS distances. The mean inter-GPS distances were calculated with data collected over the duration of all surveys and these values were used for SfM scale assignment. The average relative standard deviations of the inter-GPS distances were 2.7% for the tracking beacons and 0.01% for the high-quality GPS units. Because the same values were consistently used to scale all point clouds corresponding to a certain survey target, it was not necessary to consider the error imposed by variability in the inter-GPS distances. However, this would need to be considered if the GPS units were repositioned between surveys. Lower SfM point cloud ISE values, assigned for the M3C2 point cloud comparisons, were associated with the more precise dual-frequency GPS receivers ( Table 3).
The relatively high potential positional error in the Solara (σ bx = 2.9 m, σ by = 2.9 m) and Canatec (σ bx = 7.9 m, σ by =11 m) tracking beacons (Table 1) negatively impacted the quality of the TLS drift correction. This is seen in the larger absolute mean PC dist value of 2.2 m for TLS point clouds which were corrected with the tracking beacon data, in spite of the Kalman filtering and smoothing. The higher recording frequency and lower potential positional error of the dual-frequency GPS receivers (σ bx = 0.03 m, σ by = 0.04 m; Table 1) resulted in improved Kalman filter solutions and the subsequent drift correction (absolute mean PC dist = 1.4 m).

Detection Thresholds and Monitoring Iceberg Deterioration Processes
To measure the magnitude and spatial distribution of iceberg deterioration with a given survey technique, the deterioration that transpires between surveys must exceed DDT . This, in turn, hinges on the surveying technique (SfM vs. TLS) and geolocation accuracy (the four 'survey scenarios'), the rate of the deterioration process, the time elapsed between surveys, and the extent of the affected area. We discuss the implications of DDT and survey precision further using values from our study's survey scenarios and grouping deterioration processes into those resulting in 'chronic' or 'acute' morphological changes.
Chronic deterioration represents general and relatively slow changes over a large surface area (Figure 9). Ballicater Consulting [66] estimated the rate of horizontal melt at the waterline from wave erosion, one chronic deterioration mechanism, to be~2.7 m day −1 in conditions likely to be encountered offshore of Newfoundland and Labrador, Canada. Deterioration at this rate could be confidently detected after 1.3 days with repeat TLS surveying drift-corrected with GPS data logged by dual-frequency receivers given the calculated DDT values (Table 4). However, this same deterioration could be detected within 4 h if the aerial photography with SfM processing technique was employed. Chronic sidewall recession caused by radiation and forced convection is much harder to detect. Assuming a rate of 0.3 m day −1 [16,66] is representative, deterioration from this process could be reliably measured if employing dual-frequency GPS receivers within 8 days with TLS surveying or 1.3 days with SfM processing of aerial photography. Figure 9. Example of idealized ice loss due to the two general ice island deterioration types: (a) 'chronic' deterioration occurring between Time 1 (T 1 ) and Time 2 (T 2 ), assuming here that all sidewalls recede at a slower rate than (b) 'acute' deterioration that occurs at a specific location between T 1 and T 2 . The example of acute deterioration in (b) is the calving of an overhanging slab that is created through chronic wave erosion at the waterline.
Equation (10) can be used to estimate the amount of ice that must be lost before detecting chronic deterioration with repeat surveys given the L, W and freeboard (f ) (in m) of an idealized rectangular ice island/iceberg (Figure 9a): The equation assumes that surface recession occurs uniformly over the sidewalls of an iceberg sail. These examples relate to deterioration occurring to iceberg sidewalls that result in the change in horizontal (L and W) extents. The chronic lowering of the surface from radiation and turbulent heat fluxes [16], which affects the vertical dimension, is also possible to detect with these surveying techniques. However, the influence of buoyancy on freeboard height would need to be considered when reporting deterioration magnitudes.
Unlike chronic deterioration, acute deterioration will be immediately detectable with repeat surveying, as long as the dimensions of the calved piece are greater than the survey technique DDT (Figure 9b). Equation (11) can be used to estimate the mass loss associated with the minimum detectable acute deterioration event given the surface area (A, in m 2 ) of the calved piece, the DDT and ρ i .
The value of DDT greatly influences the mass loss required for detection. For example, assuming a 4 × 5 m growler were to calve from an iceberg sail, SfM processing of aerial photography surveys with dual-frequency GPS correction would detect this event if it weighed more than 7.0 tonnes (8.0 m 3 ), but it would need to be more than 59 tonnes (68 m 3 ) using the TLS technique with dual-frequency GPS.

M3C2 Algorithm
The ISE is used in Equation (7) to calculate the corresponding u i for each PC dist-i calculation associated with the SfM comparisons, as are ρ p and the σ i of the reference and compared point cloud surfaces. The σ i will be influenced by instrument noise (e.g., incidence angle and range noise), genuine surface roughness, GPS and INS errors, and the correct orientation of the normal vector with respect to the surface of both point clouds. The size of D is crucial for the determination of surface normal orientation. It needs to be large enough that it is not influenced by small scale surface roughness, which can result in the incorrect orientation of normal vectors and spurious distance calculations. However, large scale surface orientation changes (i.e., naturally occurring corners), such as the transition from the vertical sidewall to the horizontal upper surface of a tabular iceberg, can be occluded if D is too large.
Since no deterioration likely occurred between consecutive surveys in our study, each unique PC dist-i calculation should have the same normal orientation even if calculated with a large D patch size. If the compared survey did have a different normal orientation in comparison to the reference point cloud, this would be reflected in a large value of u i through the inclusion of σ 2i in Equation (7) ( Figure 10). Likewise, if D occludes a face orientation change on the reference point cloud this will be reflected by an increase in σ 1i . A lower ρ p for either the reference or compared point clouds will also be reflected in increased u i values (Equation 7). The possibility of occluding surface orientation changes should be considered by future surveyors who aim to detect deterioration with any surveying technique (e.g., SfM, TLS, or sonar). A method of obtaining reasonable ρ p at edges should be favoured so that u and DDT can be minimized in these critical locations. Misalignment between the surface orientations of the reference (top) and compared (bottom) point clouds is represented by thin arrows. This misalignment is captured by the large surface roughness of the compared cloud (σ 2i ) and will increase the uncertainty (u i ) of the calculated distance between the two point clouds (PC dist-i ). The reference cloud's surface roughness is represented by σ 1i . The open circles represent cloud points. The grey circle represents a core point. The black circles denote the end of the thicker line representing PC dist-i .
One of the large benefits of the M3C2 algorithm is that all points are retained for the normal orientation computation while the use of core points used for the calculation of PC dist and u decreases processing time [56]. While the use of core points was adequate for our evaluation, which did not expect deterioration to occur between surveys, future users of the M3C2 algorithm should be attentive with this subsampling step in the event that they are interested in discerning very small-scale deterioration magnitudes.

Deterioration Detection and Mass Estimation
Our analysis based the DDT calculation on the precision of 3D models generated by aerial photography surveying with SfM processing and TLS surveying. We are not aware of another study that has conducted such an analysis or determined this threshold for drifting icebergs or ice islands. The DDT values derived from the two aerial photography with SfM processing survey scenarios in our study (0.4 and 2.5 m) envelope those reported in past studies that used this technique to assess changes in the vertical dimension of their cryospheric research subjects [28,30,31]. Robe and Farmer [33] estimate a 10.5% error in their mass calculation method by propagating the errors associated with each step of their stereo-imagery processing workflow. However, these error estimates were based on static dimensions for L and W in a simplified volume estimate equation and cannot be expected to adequately represent the complex surface of an iceberg. Smith and Donaldson [67] used a sextant and sidescan sonar to survey icebergs and report a 10% error for their mass estimates. Our mass estimates generated from surveys associated with low quality tracking beacon deployment had similar error values to those reported previously (~10% to 20%). However, we saw in our study that surveying an iceberg target with SfM, paired with dual-frequency GPS unit deployment, provides an impressively high precision for mass estimation (4%).
The increasing spatial and temporal resolution provided by new satellite programs will improve the ability to monitor the deterioration of numerous icebergs and ice islands over large areas or timespans. For example, Enderlin et al. [12] use 0.55 m horizontal resolution WorldView stereo-image pairs to estimate daily linear area-averaged melt rates of icebergs with uncertainties ranging from 8% to 100%. Regardless, these estimations and the associated freshwater flux values are useful to constrain the role of iceberg meltwater on ocean properties over relatively large regions. In contrast, our technique assessment is necessary before further field based studies can be conducted to study the physical processes linked to observed deterioration.

Surveying Workflows
Future surveyors using TLS for iceberg and ice island study should consider deploying marker-GPS pairs where they can be seen and identified in the TLS point cloud for the determination of a registration error (ISE). Along with Ryan et al. [30], who found that a limited number of marker-GPS pairs inhibited the relative accuracy of SfM-generated glacier DEMs, we recommend that a number of strategically placed markers with high-quality GPS reference (i.e., GCPs) be deployed before future aerial photography surveys for georeferencing validation during SfM processing. In this study, markers were included in 60% to 80% of the images for i2's two surveys, but only between 20% and 60% for the larger i1.
Eisen [37], Fugro [68], Hamilton et al. [14], McGuire et al. [39], Younan et al. [32], and Zhou et al. [69] have all previously generated 3D iceberg maps from iceberg keel surveys conducted with multi-beam echo-sounders. Our drift correction approach can also be utilized for correcting these keel surveys. Kimball and Rock [9] developed a vessel-based terrain relative navigation method for conducting autonomous underwater vehicle (AUV)-mounted sonar surveys to map the keels of large, tabular Antarctic icebergs. Estimates of the iceberg's translation and rotation, derived relative to that of the AUV, were used to correct Doppler sonar range data and create the keel map. Alternatively, McGuire et al. [39] corrected TLS iceberg surveys with average drift velocity and rotation magnitudes derived by aligning consecutive surveys with an iterative closest point matching algorithm [39]. While GPS deployment is not a requirement when using the approaches of Kimball and Rock [9] or McGuire et al. [39], these techniques must estimate iceberg translation and rotation. The 3D model generation of a moving target is typically improved when a motion model is incorporated in the survey processing workflow to account for variations in the rate of translation and rotation [9,70]. Corrections that employ a constant drift velocity may be adequate for ice management purposes where it is necessary to rapidly process survey data to generate information regarding the size and shape of hazards [39]. However, for deterioration studies, more accurate positioning will likely be required to guarantee high-quality results. Table 5 includes a number of important considerations associated with each technique during iceberg or ice island surveying field campaigns. A valuable characteristic of SfM that was apparent in this study is the potential to log an order-of-magnitude greater ρ p than TLS [71]. This contributed to the superior deterioration detection and mass estimation capabilities of the technique. The imagery also provides a user with additional information to aid in the interpretation of the 3D models and point cloud comparisons [70]. To decrease the chances of gaps in survey coverage on this surface (Figure 6a), we recommended conducting aerial photography surveying during overcast conditions to improve contrast. This will optimize image quality and provide an adequate number of tie-points for SfM processing [19]. Another alternative is to use bathymetric LiDAR, which uses a green wavelength and would not be affected by melt ponding [31] or shadowing, and could be used to map both the ice island and shallow surrounding bathymetry (<50 m) depth during an airborne survey. Table 5. Field campaign considerations for aerial photography surveying (with subsequent structure-from-motion (SfM) processing) and terrestrial laser scanning (TLS) surveying. This list includes some of the factors that surveyors may contemplate when deciding between the use of these two survey techniques for future research programs. Costs are approximate. Further discussion on the benefits and limitations of these techniques for glaciological studies can be found in [30]. We recommend acquiring both oblique imagery with a vessel-mounted camera or with an additional aerial pass with an off-nadir camera angle along with vertical photographs to adequately cover steep surfaces [21]. Improving coverage should help to reduce the u i values associated with PC dist-i values calculated along sidewalls and edges, as seen in the SfM point comparison output for i2 ( Figure 6d). If an oblique survey is not conducted, surveyors should be cautious when interpreting small deterioration magnitudes observed around edges and sidewalls due to the increased uncertainty at these locations, as the DDT values that we provide are an average determined from the PC dist-i and u i values determined across the entire ice surfaces.

Survey Conditions & Requirements Aerial Photography/SfM Processing TLS
If future surveyors utilize vessel-based TLS surveying, it will be best suited for the observation of sidewalls and developing wave-notches due to instrument mounting location and viewing angle. A benefit of vessel-based TLS surveying is that the resulting point cloud processing takes less time compared to SfM processing once a drift correction workflow is developed [31]. Additionally, vessel-based surveying is not hindered by poor sky conditions to the same extent as helicopter-based aerial surveying. However, we believe that the benefits of SfM outweigh these limitations. Exploring the use of uninhabited aerial vehicles (UAVs) is recommended, as their use could alleviate the dependency on costly, weather-sensitive helicopters, and data collection via fixed-wing UAVs is already increasingly utilized in glacier [27,72] and iceberg or ice island [17] investigations. To truly be independent of helicopter use, one would additionally need to use an alternative method for marker/GPS deployment such as that developed by McGill et al. [17] in the Southern Ocean.
A different method for assigning scale or collecting data, such as that documented by Wang et al. [41] or Younan et al. [32], would need to be utilized if one wanted to conduct aerial photography surveying without marker/GPS deployment. It is important that the error of these alternative methods be reported for the continued development and comparison of iceberg and ice island surveying methods. For high precision surveying, we recommend that marker-dual-frequency GPS pairs are permanently deployed or re-deployed in the same position for each survey. This will ensure that the point clouds associated with each survey are consistently scaled.

Conclusions
This research provides the first quantification of the DDT for drifting icebergs or ice islands with observations made through SfM processing of repeat aerial-photo surveys or TLS surveys. In addition, it details a method of iceberg drift correction using GPS data as well as protocols for point cloud generation and comparison. The quality of GPS units deployed on the iceberg surface for drift correction greatly influenced the precision of the survey data, which yielded a large range (0.40 to 6.6 m) in DDTs (Table 4). We also calculate the relative uncertainty of our iceberg mass estimates based on our repeat survey data and found that we can repeatedly derive an iceberg's mass within 4% when using aerial photography and SfM processing techniques. These survey techniques can now be utilized in the offshore environment to provide upper and lower bound estimates of an iceberg's mass based on the uncertainties associated with the particular surveying scenario.
The two surveying techniques differed in terms of their capability to detect specific deterioration processes, along with logistical, instrumental, and platform considerations for field campaigns (Table 5). Based on our experience, we recommend that ice island or iceberg field campaigns utilize SfM processing of aerial photography surveys for deterioration observations after considering survey target characteristics, environmental conditions, and campaign objectives. This is due to the decreased cost of equipment, additional visual information provided by imagery, high point density, relatively lower noise of point clouds, the lower DDT values, and better mass estimation capabilities of this technique that will allow for confident reporting of morphological changes in greater detail. To reduce reliance on helicopters, we also recommend the integration of one or more of the following into an iceberg aerial photography survey system: (1) oblique vessel mounted cameras, (2) UAV platforms for aerial photography survey, and (3) UAV tracking beacon delivery or an alternative method for scaling SfM point clouds, such as reported by Wang et al. [41]. This will dramatically decrease cost, allow for the use of survey vessels without helidecks, and alleviate some weather limitations associated with helicopter flying. An off-nadir circumnavigation of the survey target's sidewalls should also be added to improve sidewall coverage for waterline deterioration detection. Surveying campaigns may also include the use of multi-beam sonar [9,39] to both validate the volume and mass estimations made with the sail surveys as well as to capture the wave and thermal erosion occurring at the ice-water-air interface.
Future deterioration measurements via repeat surveys from a dedicated vessel, along with concurrent environmental observations (e.g., current, wind, temperature), will provide much needed in-situ data for iceberg drift and deterioration analysis and model development. These efforts can then be applied to mitigate the risks associated with offshore industrial operations, as well as for fundamental research into the physical and ecological implications of the influx of iceberg meltwater into the ocean. Offshore stakeholders who are concerned with ice hazards in a more immediate sense can utilize the surveying techniques to provide mass estimates while being aware of the error bounds of these estimates. With this information and the auxiliary considerations discussed in this paper, ice island and iceberg surveyors can now plan field campaigns (e.g., technique, platform, GPS quality) based on the subject of interest (i.e., deterioration modelling or mass estimation) and the required quality of the end product.
Supplementary Materials: The data and computer scripts used in this study are publically available on the Polar Data Catalogue and can be found at https://www.polardata.ca/pdcsearch/PDCSearch.jsp?doi_id=12818. laser scanner and assisted during TLS surveys. The advice and assistance during survey preparation from Josh Sampey and Mike Mutschler of Seahorse Geomatics was essential for our field success. The external helicopter mount was provided by the Center for Earth Observation Science at the University of Manitoba, Winnipeg, Manitoba, Canada. We recognize the continued support for ice island research provided by the Canadian Ice Service (CIS), Environment and Climate Change Canada, Ottawa, Ontario, Canada. CIS employees Matt Arkett, Angela Cheng, Hugo Drouin, and Leah Braithwaite provided updates on potential survey target locations throughout the field campaign. The expertise of Luc Desjardins and Ron Saper from the Water and Ice Research Lab at Carleton University, Ottawa, Ontario, Canada also contributed to the successful campaign. We thank Greg Crocker, Doug King, and Koreen Millard as well as several anonymous reviewers for insightful comments on earlier versions of this paper. Greg Crocker also provided computer code and background information regarding the Kalman filter.