2.1. Study site

Beacons were deployed from the Canadian Coast Guard Ship (CCGS) Amundsen, a research icebreaker, during the Offshore Newfoundland Research Expedition (ONRE) from 17 April to 4 May 2015. The deployment of the tracking beacons took place on 29 April 2015 offshore northern Newfoundland, where a large ice island (most likely from the Petermann Glacier in NW Greenland) was detected in Landsat imagery at the NE entrance to the Strait of Belle Isle. At this time, the ice island was grounded near the southeastern coast of Labrador and locked in 9+/10th sea-ice concentration. The dimensions of the ice island were estimated to be ~1.6 km in maximum waterline length (with 6% uncertainty) and ~10 m in average freeboard; the latter was estimated using the hydrostatic balance for the parent ice island (Fig. 2).

Fig. 2. The ice island in this study, shown surrounded by sea ice at the NE end of the Strait of Belle Isle on 29 April 2015.

2.2. Tracking beacon deployments and characteristics of ice island fragments

Four beacons (Canatec and Associates, Ltd) enclosed in waterproof, plastic cases were deployed by hand via the ship's helicopter on the ice island in a quadrilateral array 500–700 m apart from one another toward the center of ice island (Table 1). The beacons were secured to the ice with boards that had beds of nails in them, to prevent the beacons from being blown off the ice island by the wind. The ice-tracking beacons transmitted Universal Time Coordinated (UTC) time, internal temperature and latitude/longitude GPS readings via Iridium satellite to a data management website (with 1.8 m Circular Error Probable). Initially, the four beacons were remotely set to transmit data at 6-hour intervals during the time that the ice island was grounded. The ice island broke up into several fragments (Fig. 1) on 6 May 2015, which began to drift SW further into the Strait of Belle Isle. The four tracking beacons were then set to transmit hourly data and remained on the ice island fragments until the fragments melted (or experienced calving or rolling) (Table 1). The tracking beacons continued to transmit data as surface current drifter buoys, but this study only presents the analysis during the time that the beacons remained on the fragments. To determine whether the beacon was on the ice island or freely drifting in the ocean, diurnal variability in beacon temperatures was analyzed. While on ice, the beacon temperature records consistently displayed large diurnal variations. However, when the beacons fell into water, their temperature records reflected surface water temperature with minor fluctuations.

Table 1. Summary of ice island tracking beacon deployments (times and locations where the beacons were first deployed) on a large ice island grounded at the NE end of the Strait of Belle Isle

^a International mobile equipment identity.

^b Coordinated universal time.

2.3. Reduction in surface area

The surface area, length and width of each ice island fragment in this study were estimated from the available (29 April 2015–21 July 2015) SAR images during the fragment drift, using Focus PCI Geomatics software (version 2013). These dimensions were interpolated between satellite observations based on the assumption that their reduction in size was linear. However, it is possible that the fragments may have undergone episodic calving between observations.

2.4. Surface ablation

The surface ablation of the ice island fragments was estimated using a TIM model, in which the surface melt rate is calculated by (Hock, Reference Hock2003):

(1)$$V_{{\rm sa}} = DDF\mathop \sum \nolimits PDD$$

where DDF (Table 2) represents the degree-day factor, which is a region-specific parameter that the TIM model is sensitive to (Crawford and others, Reference Crawford, Mueller, Humphreys, Carrieres and Tran2015). Reported values of degree-day factors range from 2 to 11 mm d⁻¹ °C⁻¹ for glaciers in various sites (Hock, Reference Hock2003). This study uses an average of 4 mm d⁻¹ °C⁻¹ for the degree-day factor (Table 2), which was estimated based on the assumption that the ice island fragments became completely melted (through surface and basal ablation) at the same time that the beacons fell into water. While the beacons could have fallen into the water due to the break-up or capsize of the fragments, this is unlikely given that the available SAR images did not show any remaining fragments beyond the time we assumed the beacons fell into water. The variable PDD in Eqn (1) represents the positive degree days, which is the sum of the daily mean positive air temperatures (>0°C) during the analyzed period (Hock, Reference Hock2003).

Table 2. Description and values of the parameters and variables

^a A 40% higher value (0.56) was used for c _a to implicitly account for the effect of ocean surface waves (Smith, Reference Smith1993; Keghouche and others, Reference Keghouche, Bertino and Lisæter2009).

2.5. Basal ablation

The thinning rate through the base of the ice island (V _ba) was determined using the forced convection basal ablation model that was calibrated with ice island melt rate observations and recorded atmospheric and oceanic data, given by (Crawford, Reference Crawford2018):

(2)$$V_{{\rm ba}} = k\vert {{\vec{v}}_{\rm w}(b)-\vec{u}} \vert ^{0.8}\left( {\displaystyle{{T_{\rm w}(b)-T_{{\rm mi}}} \over {L^{0.2}}}} \right)$$

where k is an empirical coefficient (Table 2) given by Crawford (Reference Crawford2018), $\vec{v}_{\rm w}(b)$ is the ocean current velocity at the keel depth, u is the ice island fragment velocity, T _w(b) is the water temperature at the base of the fragment, T _mi is the melting temperature of the ice and L is the length of the ice island fragment. The fragment velocity was calculated between successive beacon position/time records. The calculated values represent the velocity at the midpoint between these positions, so they were linearly interpolated in space and time back to the beacon record positions. The melting temperature of the ice was estimated using an empirical correlation to account for the ice melting point depression due to the interfacial salinity at the base of the ice island fragments (Løset, Reference Løset1993; Kubat and others, Reference Kubat, Sayed, Savage, Carrieres and Crocker2007), given by Josberger (Reference Josberger1978):

(3)$$T_{{\rm mi}} = T_{\rm f}{\rm exp}\lpar {-0.19\lsqb {T_\infty -T_{\rm f}} \rsqb } \rpar $$

where T _∞ indicates the ambient water temperature at a distance beyond the influence of the ice island fragment. It should be noted that T _∞ and T _w(b) are effectively the same in our model given that the extracted ocean water temperature from CMEMS has a resolution of 1/12° (~9 km). So the ambient water temperature at some distance away from the ice island fragments lies within the same gridpoint as the ice island where T _w(b) was extracted. T _f is the freezing temperature of the water, which was estimated as a function of the water salinity (S) given by (Fujino and others, Reference Fujino, Lewis and Perkin1974):

(4)$$T_{\rm f} = -0.036\,-\,0.0499S\,-\,0.000112S^2$$

The daily water salinity and potential temperature at different depth layers were extracted using Global Ocean Physics Reanalysis Glorys12v1 model from Copernicus Marine Environment Monitoring Service (CMEMS). Potential temperature was then adjusted based on the water pressure and salinity at each depth to calculate the actual water temperature. The velocity of ocean currents over the depth of ice island fragments was extracted from CMEMS Global Ocean 1/12° Physics Analysis and Forecast model. The extracted data were all interpolated to the fragments' observed locations and then integrated into Eqn (2) to allow for the calculation of basal ablation rates.

To present the total change in the thicknesses of ice island fragments over time, the cumulative vertical melt was estimated according to the above surface and basal ablation models and subtracted from the initial thickness. The initial thickness of the original ice island was estimated using the ice island hydrostatic balance at the time it started to drift. It was assumed that once the beacons started drifting, the ice island would have a draft slightly less (by 5 m) than the shallowest water depth under the ice island.

2.6. Mass calculation

The mass of the fragments was calculated by:

(5)$$M = \rho _{\rm i}V_{\rm i}$$

where ρ _i is the density of the ice island fragments. While the density of pure (bubble-free) ice is ~917 kg m⁻³ (Enderlin and Hamilton, Reference Enderlin and Hamilton2014), ice islands have lower density. This study uses the ice density of 873 kg m⁻³ for the fragments, as measured by Crawford and others (Reference Crawford, Mueller and Joyal2018b). The variable V _i in Eqn (5) represents ice island volume, calculated by:

(6)$$V_{\rm i} = A_{{\rm ha}}h_{\rm i}$$

where A _ha and h _i represent the ice surface area and thickness of the ice island, respectively. The mass of each fragment was recalculated over time based on the interpolated surface area and modeled thickness.

2.7. Drift equation

The equation describing the drift of an ice island fragment is stated as (Lichey and Hellmer, Reference Lichey and Hellmer2001):

(7)$$M\displaystyle{{{\rm d}\vec{u}} \over {{\rm d}t}} = \vec{F}_{\rm a} + \vec{F}_{\rm w} + \vec{F}_{\rm c} + \vec{F}_{{\rm ss}} + \vec{F}_{{\rm si}}$$

where M is the mass of the ice island fragment. The five terms on the right side of Eqn (7) represent forces due to air drag, water drag, Coriolis effect, sea surface tilt, and sea ice, respectively. To evaluate these forces when the ice island fragments were not grounded, the atmospheric and oceanic data were extracted from available databases; the 3-hourly values for 10-m wind velocity and 2-m air temperature were extracted from NARR dataset, however, sea surface height and sea-ice velocity were extracted from CMEMS Global Ocean 1/12° Physics Analysis and Forecast model. There were time periods in the extracted ocean currents data for the four fragments during which the ocean current velocities were missing. Therefore, the extracted data were linearly interpolated to fill out the gaps in the velocity of ocean currents. While the mass of the ice island was assumed to be constant in some drift models (e.g. Lichey and Hellmer, Reference Lichey and Hellmer2001), it varied in our model (Eqn (7)) to improve the model accuracy.

Air drag and water drag force (Eqn (7)) both on vertical walls and horizontal surfaces of the ice island fragment were calculated by (Smith and Banke, Reference Smith and Banke1983):

(8)$$\vec{F}_{\rm a} = \left[ {\displaystyle{1 \over 2}\lpar {\rho_{\rm a}c_{\rm a}A_{{\rm va}}} \rpar + \rho_{\rm a}c_{{\rm da}}A_{{\rm ha}}} \right]\vert {{\vec{v}}_{\rm a}-\vec{u}} \vert \lpar {{\vec{v}}_{\rm a}-\vec{u}} \rpar $$

(9)$$\vec{F}_{\rm w} = \left[ {\displaystyle{1 \over 2}\lpar {\rho_{\rm w}c_{\rm w}A_{{\rm vw}}} \rpar + \rho_{\rm w}c_{{\rm dw}}A_{{\rm hw}}} \right]\vert {{\vec{v}}_{\rm w}-\vec{u}} \vert \lpar {{\vec{v}}_{\rm w}-\vec{u}} \rpar $$

where ρ _a is the air density, c _a is the coefficient of resistance (form drag) for air and c _da is the air drag coefficient (skin drag). In Eqn (8), the vertical and horizontal areas of the ice island fragment exposed to air flow $(\vec{v}_{\rm a})$ are denoted by A _va and A _ha, respectively. The parameters and variables in Eqn (9) are analogous to those in Eqn (8), but are indicated by the subscript w for water. Given the tabular shape of ice islands, the average current velocity $(\vec{v}_{\rm w})$ in Eqn (9) was defined as the average of current velocity profile over the fragment draft. A _va and A _vw were estimated from the fragment length (L), width (W), sail height (h _s), and draft (h _k) by:

(10)$$A_{{\rm va}} = h_{\rm s}\displaystyle{{\lpar {L + W} \rpar } \over 2}$$

(11)$$A_{{\rm vw}} = h_{\rm k}\displaystyle{{\lpar {L + W} \rpar } \over 2}$$

The resulting force from the Coriolis effect and sea surface tilt are given by (Lichey and Hellmer, Reference Lichey and Hellmer2001):

(12)$$\vec{F}_{\rm c} = 2M\Omega \;\sin \phi \vec{k} \times \vec{u}$$

(13)$$\vec{F}_{{\rm ss}} = M\vec{g}\sin \alpha $$

where Ω (7.27 × 10⁻⁵ rad s⁻¹), ϕ and $\vec{k}$ are the Earth's angular speed, the latitude of the ice island fragment and the unit vector perpendicular to the surface of Earth, respectively. In Eqn (13), g is the acceleration due to gravity 9.81 m s⁻² and α represents the sea surface slope. The sea surface tilt (sinα) was estimated directly from CMEMS satellite altimetry data.

The force due to the surrounding sea ice was calculated using two different approaches; the residual method and Lichey and Hellmer's model (Reference Lichey and Hellmer2001). From Eqn (7), the residual force $(\vec{F}_{{\rm res}})$ in the presence of sea ice was estimated by:

(14)$$\vec{F}_{{\rm res}} = \vec{F}_{{\rm si}} + \vec{e}_{\rm i} = M\displaystyle{{{\rm d}\vec{u}} \over {{\rm d}t}}-\lpar {{\vec{F}}_{\rm a} + {\vec{F}}_{\rm w} + {\vec{F}}_{\rm c} + {\vec{F}}_{{\rm ss}}} \rpar $$

where $\vec{F}_{{\rm si}}$ is the sea-ice force and $\vec{e}_{\rm i}$ is the error term that accounts for the error in the input data. In the presented analysis, it was assumed that the error term, $\vec{e}_{\rm i}$, was sufficiently small, then $\vec{F}_{{\rm res}}\~\vec{F}_{{\rm si}}$. This allowed us to estimate the contributions from sea ice based on the residual force. In the case that there was no sea ice present $(\vec{F}_{{\rm si}} = 0)$, Eqn (14) still holds, but the residual force represents only the error term between the net force and the estimated force components. The residual force in open water is therefore given by:

(15)$$\vec{F}_{{\rm res}} = \vec{e}_{\rm i} = M\displaystyle{{{\rm d}\vec{u}} \over {{\rm d}t}}-\lpar {{\vec{F}}_{\rm a} + {\vec{F}}_{\rm w} + {\vec{F}}_{\rm c} + {\vec{F}}_{{\rm ss}}} \rpar $$

The force due to the surrounding sea ice was calculated by (Lichey and Hellmer, Reference Lichey and Hellmer2001):

(16)$$\vec{F}_{{\rm si}} = \left\{ {\matrix{ 0 \hfill & {C \le 15\%} \hfill \cr {\displaystyle{1 \over 2}\rho_{{\rm si}}c_{{\rm si}}A_{{\rm si}}\vert {{\vec{v}}_{{\rm si}}-\vec{u}} \vert \lpar {{\vec{v}}_{{\rm si}}-\vec{u}} \rpar \;} \hfill & { 15\% \lt C \le 90\%} \hfill \cr {-\lpar {{\vec{F}}_{\rm a} + {\vec{F}}_{\rm w} + {\vec{F}}_{\rm c} + {\vec{F}}_{{\rm ss}}} \rpar } \hfill & {C \ge 90\% \; \;{\rm and}\;P \ge P_{\rm s}} \hfill \cr}} \right.$$

where ρ _si is the density of sea ice, c _si is the coefficient of resistance for sea ice, assumed to be 1.0 (Bigg and others, Reference Bigg, Wadley, Stevens and Johnson1997; Lichey and Hellmer, Reference Lichey and Hellmer2001), A _si is the estimated contact area between sea ice and the ice island fragment (defined as the product of sea-ice thickness and ice island waterline length; Lichey and Hellmer, Reference Lichey and Hellmer2001), $\vec{v}_{{\rm si}}$ is the sea-ice velocity and C is the concentration of sea-ice cover in percent. According to Eqn (16), sea ice with concentrations <15% are considered as open water and exert no force on the ice island fragment. However, for sea-ice concentrations between 15 and 90% (open drift), sea ice impose additional drag on the ice island. For concentrations over 90% (close pack) where the sea-ice compressive strength (P) is greater than its threshold value (P _s), sea-ice envelopes the ice island fragment causing it to drift at the sea-ice velocity. The sea-ice compressive strength can be estimated from its thickness (h _si) and concentration by (Hibler, Reference Hibler1979):

(17)$$P = P^{^\ast}h_{{\rm si}}{\rm e}^{-a(1-C)}$$

where P* and a are empirical coefficients with values of 20 000 N m⁻² and 20, respectively (Lichey and Hellmer, Reference Lichey and Hellmer2001). While Eqn (16) presents the sea ice forcing for different sea-ice concentrations, only the relationships associated with open water and open drift were used in this study as the concentration of sea ice in the vicinity of ice island fragments during the analyzed period were never beyond 90%. Sea-ice concentration and thickness were extracted manually from the CIS daily ice charts and cross-checked against the SAR images at the location of the ice island fragments over their drift period. A summary of all the parameters and variables used in this study, along with their values, is given in Table 2.

2.8. Sensitivity analysis

While commonly used values of skin and form drag coefficients for air and water (c _da = 0.00025, c _dw = 0.0005, c _a = 0.4, and c _w = 0.85) were used in this study (Table 2), it is important to examine how sensitive the residual sea-ice force results are to a change in these coefficients. Due to the roughness variation in the surfaces, bases and walls of the ice island fragments, the skin drag coefficients were varied between 0.00025 and 0.008, while the form drag coefficient values ranging from 0.05 to 2 were tested. These values were selected to account for the variation among the reported air/water drag coefficients on icebergs in other studies (e.g. Lichey and Hellmer, Reference Lichey and Hellmer2001; Kubat and others, Reference Kubat, Sayed, Savage and Carrieres2005; Eik, Reference Eik2009; Keghouche and others, Reference Keghouche, Bertino and Lisæter2009; Rackow and others, Reference Rackow2017), where values in the ranges of 0.00025–0.0025, 0.0005–0.0055, 0.1–2 and 0.05–1.5 were used for the air skin drag, water skin drag, air form drag and water form drag coefficients, respectively. To examine the individual influence of a varying drag coefficient, the other three coefficients were fixed at the nominal values mentioned earlier. These values were integrated into the residual sea-ice force model (Eqn (14)) to examine which of these four drag coefficients the model is most sensitive to.

Due to the uncertainty in the estimated masses and ocean currents, the sensitivity analysis was also performed to study the residual sea-ice force results to the variation of ice island fragments' masses, as well as the mean ocean currents speed and direction. Given that the dimensions of the fragments were estimated from satellite imagery data with resolution of 100 m, the propagated error for mass was estimated and integrated into the sensitivity analysis. Also, the mean error in the magnitudes (speed) and directions of surface currents from CMEMS Global Ocean 1/12° Physics Analysis and Forecast model was reported to be 0.08  m s⁻¹ and 7.2°, respectively (Lellouche and others, Reference Lellouche2018). Assuming the same error for the currents in different layers beneath the water surface, the residual forces (Eqn (7)) were calculated using the estimated average error in speeds (± 0.08 m s⁻¹; high and low) and directions (± 7.2°; high and low) of the currents to evaluate the model sensitivity to the error in the CMEMS ocean currents.