1. Introduction
Climate change has had a dramatic impact on the Arctic, which is warming three to four times faster than that of the global average (Bush and Lemmen, Reference Bush and Lemmen2019; IPCC, 2022; Rantanen and others, Reference Rantanen2022). This rapid increase in atmospheric and oceanic temperature has resulted in many changes to cryospheric features such as ice shelves and floating ice tongues, which now calve more frequently compared to in previous decades (Rignot and others, Reference Rignot, Velicogna, van den Broeke, Monaghan and Lenaerts2011; Bevis and others, Reference Bevis2019). At the northern coast of Ellesmere Island, Canada’s ice shelves underwent nine major calving events between the years 2000 and 2020 (Mueller and others, Reference Mueller, Copland, Jeffries, Copland and Mueller D2017; Vincent and Mueller, Reference Vincent and Mueller2020), but used to calve approximately once per decade (Jeffries and Sackinger, Reference Jeffries, Sackinger, Murthy, Paren, Sackinger and Wadhams P1990). The reduction of multi-year land fast sea ice, a protective barrier to the ice shelves fringing northern Ellesmere Island (Copland and others, Reference Copland, Mueller and Weir2007; Mueller and others, Reference Mueller, Copland, Jeffries, Copland and Mueller D2017), and open water at the ice interface are thought to play a role in accelerating the calving activities in this region (Reeh and others, Reference Reeh, Thomsen, Higgins and Weidick2001; Copland and others, Reference Copland, Mueller and Weir2007). Increased calving frequency of western Greenland’s floating glacier tongues has also been noted, with ice shelves and ice tongues calving icebergs into the Arctic Ocean more frequently (Bigg and others, Reference Bigg2014; Benn and others, Reference Benn, Cowton, Todd and Luckman2017; Shepherd and others, Reference Shepherd2020).
Ice islands, a type of iceberg that is large and tabular, originate from ice shelf and floating ice tongue calving events, such as the those from Petermann Glacier in 2008, 2010 and 2012. These events occur erratically compared to icebergs, but they can fracture into many smaller ice islands over a wide area (Crawford, Mueller, Desjardins, and others Reference Crawford, Mueller, Desjardins and Myers2018; MacKie and others, Reference MacKie, Millstein and Serafin2024), which poses a risk to offshore infrastructure and shipping vessels (Eik and Gudmestad, Reference Eik and Gudmestad2010; Bailey and Phillips, Reference Bailey and Phillips2018). Though ice islands may remain adrift or grounded in the Arctic Ocean, a considerable number of them drift as far south as the Grand Banks region of Newfoundland where oil-drilling and production platforms are located (Fuglem and Jordaan, Reference Fuglem, Jordaan, Copland and Mueller2017). As the Arctic warms, extensive reductions in sea ice during summers (Stroeve and others, Reference Stroeve2008; Overland and Wang, Reference Overland and Wang2013) may spark interest in new offshore oil exploration projects at higher latitudes where collisions with ice islands are a heightened risk (Fuglem and Jordaan, Reference Fuglem, Jordaan, Copland and Mueller2017). Further, the shipping industry is poised to take advantage of shorter trans-oceanic routes through the Northwest Passage or directly across the Arctic Basin where these glacial hazards may also exist (Smith and Stephenson, Reference Smith and Stephenson2013).
Collisions between ice hazards and vessels needs to be considered in the context of these environmental changes (Pizzolato and others, Reference Pizzolato, Howell, Derksen, Dawson and Copland2014; Dawson and others, Reference Dawson, Pizzolato, Howell, Copland and Johnston2018) as historically, accidents are not unprecedented. The Ship Collision with Iceberg Database contains over 670 records of incidents over the last two centuries (Hill, Reference Hill2001). In March of 2017 an iceberg drifted within 180 m of the SeaRose floating production, storage and offloading (FPSO) vessel as it was extracting oil 350 km east of St. John’s, NF (O’Neill-Yates, Reference O’Neill-Yates2018). Operators of the SeaRose had knowledge of the incoming hazard but protocols to avoid a collision were not followed (O’Neill-Yates, Reference O’Neill-Yates2018). Though no accident occurred, the lack of action led to suspension of the SeaRose for a year (O’Neill-Yates, Reference O’Neill-Yates2018), bringing attention to the seriousness of ice hazard risk management and regulation.
Ice island deterioration processes play an important role in risk assessment and management. Ice islands sporadically calve smaller fragments that can be difficult to detect, but are still considered a hazard to some offshore activities (Van Wychen and Copland, Reference Van Wychen, Copland, Copland and Mueller D2017; Akbari and Brekke, Reference Akbari and Brekke2018). Sizes of calved fragments range in length from very large (>200 m long), large (121–200 m), medium (61–120 m) and small icebergs (15–60 m) as well as bergy bits (5–15 m) to growlers (<5 m) (Canadian Ice Service, 2005). Continuous (e.g. thinning from ablation) and frequent small-scale deterioration processes (such as edge-wasting) are relatively well documented (Crawford, Reference Crawford2013), however, large-scale calving is rare and episodic. Consequently, there are only a few studies of these events (Wagner and others, Reference Wagner2014; Zeinali-Torbati and others, Reference Zeinali-Torbati, Turnbull, Taylor and Mueller2021), with many focussing on Antarctic icebergs (which are typically tabular and analogous to ice islands (Stern and others, Reference Stern, Adcroft, Sergienko and Marques2017; Bouhier and others, Reference Bouhier, Tournadre, Rémy and Gourves-Cousin2018; England and others, Reference England, Wagner and Eisenman2020). Filling this knowledge gap would improve deterioration models and, coupled with drift models, improve our understanding of how these ice hazards are created and distributed.
2. The ‘footloose’ calving mechanism
Exposure of ice island sidewalls to waves and water temperatures above the freezing point at the waterline leads to the formation of a notch in the ice, which grows until the overhanging ice can no longer be supported and breaks off (Wagner and others, Reference Wagner2014). This creates an underwater terrace called a ‘ram’ (Fig. 1), which increases in size as the surface of the ice island sail shrinks from progressive edge-wasting above the notch (Wagner and others, Reference Wagner2014). The development of a ram creates a local hydrostatic disequilibrium as the buoyancy of the completely submerged ram pushes upward, causing a bending stress at the bottom of the ice island. If the ram is large enough, localized tensile stresses (membrane stress) may exceed the yield strength of the ice and result in flexural failure whereby the resultant fracture propagates upwards to the top of the sail, creating separate ice fragments. This calving process was identified by Diemand and others (Reference Diemand, Nixon and Lever1987), and dubbed the ‘Footloose Mechanism’ by Wagner and others (Reference Wagner2014). Wagner and others (Reference Wagner2014) simulated this calving process in 1-D and Mosbeux et al. (Reference Mosbeux, Wagner, Becker and Fricker2020) did so in 2-D.
Diagram of an ice island. The sail is the mass that remains above the waterline and the keel is the mass submerged below the water. The freeboard is the height of the ice island sidewall above the waterline while the draft is the vertical depth below it. The ram length is the distance it protrudes away from the sidewall of the ice island while the extent is the length of the perimeter of the ice island that the ram occupies. Thickness of the ram is roughly equal to the keel depth or draft. Note that the horizontal extent of an actual ice island is much greater than depicted here in this diagram.
The footloose mechanism is also thought to be a driver of iceberg calving from glaciers and ice shelves (Mosbeux and others, Reference Mosbeux, Wagner, Becker and Fricker2020). For example, Becker and others (Reference Becker, Howard, Fricker, Padman, Mosbeux and Siegfried2021) and Scambos and others (Reference Scambos, Sergienko, Sargent, MacAyeal and Fastook2005) observed depressions near the edge of the Ross and Ronne ice shelves in Antarctica using ICESat altimeter measurements. This, they postulated, was related to upward bending of ice from underwater rams at their termini. As buoyancy drives the rams upward, a depression forms in the upper surface of the ice some distance from the glacier edge parallel to the calving front. Sartore and others (Reference Sartore, Wagner, Siegfried, Pujara and Zoet2025) also speculate this mechanism is a significant driver of smaller-scale calving events of these ice shelf fronts.
While there is a strong interest in the footloose process (sometimes referred to as ‘super buoyancy calving’) from Antarctic ice shelves, there is less research on this calving mechanism in their thinner Arctic counterparts. Trevers and others (Reference Trevers, Payne, Cornford and Moon2019) postulated a similar mechanism could explain the rapid calving of icebergs from Ilulissat Glacier (Jakobshavn Isbrae) in Greenland. Benn and others (Reference Benn, Cowton, Todd and Luckman2017) modelled calving events from undercuts occurring at the base of tidewater glacier termini, which creates stress via a hydrostatic imbalance at the bottom of the glacier near the calving front. Our approach to investigating the footloose mechanism is broadly applicable to glaciers and ice shelves in the polar regions.
A 3-D finite element analysis (FEA) model was developed by Sazidy and others (Reference Sazidy, Crocker and Mueller2019) to study the effect of rams on ice island calving. This new model can be used to evaluate the influence of different configurations of ram length, width and extent on stress build-up and fracture within ice islands to better understand footloose calving events. Ice island calving events are rarely observed in situ, with only a single study reporting an observed footloose calving event (Wagner and others, Reference Wagner2014). While some work has been done to ascertain the shape of keels and rams in icebergs (Wang and others, Reference Wang, deYoung and Bachmayer2015; McGuire and others, Reference McGuire2016; Zhou and others, Reference Zhou, Bachmayer and deYoung2019), those of ice islands remain poorly understood. This is, partly because they are less common than icebergs, but they are also more remote, making them time-consuming to reach while their large size makes comprehensive surveying of their sides and underside technically challenging (Forrest and others, Reference Forrest2012). However, the novel Canadian Ice Island Drift, Deterioration, and Detection (CI2D3) Database catalogues calving events and observations of drifting ice islands in the Eastern Arctic Ocean between 2008 and 2013 (Crawford, Crocker, and others, Reference Crawford, Crocker, Mueller, Desjardins, Saper and Carrieres2018), and presents an excellent opportunity to evaluate the influence of ram size and shape on many observed ice island calving events in spite of a lack of measurements of these features.
We leverage the observations of ice islands before and after calving in the CI2D3 Database along with an FEA model that simulates the physics of footloose calving. To overcome the current lack of ice island ram observations in the literature, we develop a custom script to simulate ice islands with various ram configurations and compare the results to observed calving, analogous to back analysis, used in geotechnical engineering to determine the likely value of parameters following failure events in natural materials (Sakurai, Reference Sakurai2017).
Our study examined the location and magnitude of stress in simulated three-dimensional ice islands using FEA. To summarize the relationships between ice island morphology and stress magnitude we derived an empirical model from our FEA results to explore which ram characteristics are most important and to work towards improved predictions of footloose-type calving.
3. Methods
To simulate ice island stress and fracture, ice islands were selected from the CI2D3 Database (Section 3.1) and their 3-D shape was rendered into a mesh using an assumed thickness and ram shape (Sections 3.2 and 3.3). The hydrostatic pressure on these ice island models was computed (Section 3.4) and they were assigned the various parameters (Section 3.5) required to compute the evolution of stress in 3 dimensions within each mesh element over simulation time (Section 3.7). Here, we focus on the maximum principal stress (MPS), which is a scalar representing the highest amount of stress within a given element at a given time, where positive MPS is a tensile stress and negative MPS is compressive stress. Ice failure is assumed to occur when the MPS exceeds a prescribed ice strength parameter. When not modelling failure, we can calculate the maximum MPS (mMPS) across all elements and times in a given simulation. Since the mMPS is a single number representing the stress due to a particular simulated ram configuration, we created an empirical relationship between stress and the shape characteristics of the ice island as a way of summarizing the results, and inferring plausible ram shapes and sizes (Section 3.9). Stress and deformation in a test beam and a selection of smaller, simplified rectangular ice island meshes were also created for comparison to analytical solutions.
3.1. Ice island selection
The CI2D3 Database (V1.1; https://www.polardata.ca/pdcsearch/PDCSearch.jsp?doi_id=12678) catalogues the drift and deterioration of ice islands from three major calving events from Petermann Glacier and other floating ice tongues of northwestern Greenland (Ryder, Steensby and CH Ostenfield glaciers) (Crawford Crocker, and others, Reference Crawford, Crocker, Mueller, Desjardins, Saper and Carrieres2018). The database contains over 25 000 ice island geospatial polygons, traced from RADARSAT-1, RADARSAT-2 and ENVISAT images, as they drifted from Nares Strait to as far as Newfoundland (Crawford Crocker, and others, Reference Crawford, Crocker, Mueller, Desjardins, Saper and Carrieres2018).
In the CI2D3 Database, ice island observations are linked to their ‘parent’ fragment through the lineage field in the CI2D3, which was queried for all ice island observations that underwent calving prior to their subsequent database entry (n = 397). The 2-D polygon of each parent ice island was plotted along with the two or more ‘child’ fragments that were produced via calving to visualize the results of each fracture event in a ‘calving plot’ (Fig. 2a). These were analyzed qualitatively to eliminate all fracture events that were not relevant to this study. For example, ice islands that fractured near the centre were discarded since this mode of calving is thought to result from large internal flaws unrelated to the footloose mechanism (Diemand and others, Reference Diemand, Nixon and Lever1987). Only parent ice islands that calved near their edges and produced a child fragment that was approximately 20% or less of the surface area of the original parent were retained. Ice islands trapped in sea ice or that were grounded were not analyzed, since calving may have been influenced by substantial non-hydrostatic loads which are not accounted for in the fracture model we used. To lower mesh processing to a reasonable computational time, the ice island subset was reduced to ice islands with a length of ≤7.5 km, leaving 26 ice islands in total.
Example of a calving plot and an ice island mesh for ice island BSS. (a) Plan view of the 1.37 km2 ice island BSS on 5 July 2011 prior to calving, with child fragments NUA (1.24 km2) and NCJ (0.12 km2) on 8 July 2011 after calving. The dashed red line delineates the ends of the BSS calving edge based on the two fragments, while the blue dotted line indicates the edge of NCJ. (b) Isometric view of a completed mesh of BSS with a ram length of 60 m (green arrow) along the full extent (blue arrow) of its calving edge.
3.2. Ram extent
We refer to the distance the ice island ram measures along the ice island edge as the ram extent (Fig. 1), although we recognize there is no universal terminology for this feature, with (Crawford and others, Reference Crawford, Crocker, Smith, Mueller and Wagner2024) referring to it as the failure width. The length of the ram refers to the distance it protrudes away from the ice island perpendicular to the edge of the sail. Two broad types of ram extents were simulated in this study: (1) uniform rams—which fringe the entire perimeter of an ice island and (2) isolated rams—which extend along a particular portion of the ice island perimeter. To recreate the footloose calving, isolated rams were simulated where the ice island ultimately calved, along the calving edge. Here, we define the calving edge as the edge of the ice island that is lost to the smaller child fragment(s) following a calving event. Calving edges were determined by examining the calving plot of each ice island under study (Fig. 2a and b).
Uniform rams were simulated initially, but our study focused on the effect of isolated rams, which were simulated in two different ways: a full extent where the ram was created along the entire calving edge and a half extent version spanning only the middle portion of the calving edge.
3.3. Meshing and ram creation
A custom workflow in R (Version 3.3.2), a programming language for statistical computing and graphics generation (R Core Team, 2022) was used in conjunction with LS-PrePost (a software program that creates meshes and visualizes FEA model results for LS-DYNA Solver (https://www.ansys.com/) to pre-process the polygons of ice island extents into 3-D meshes. Each mesh is composed of quadrilateral elements that are delineated by vertices, referred to as ‘nodes’. The R scripts used to process and solve ice island meshes may be accessed at https://github.com/jsmith2-wirl/icemesher.
The 2-D polygons of 26 selected ice islands were extracted from the CI2D3 Database using R and converted to a format that LS-DYNA could compile. For each experiment, three ram lengths were prescribed (rl = 20, 40 and 60 m). The 2-D shape was meshed at a 20 × 20 m resolution (horizontal mesh size = 20 m) and extruded to prescribed thicknesses (ice island thickness = 80 or 100 m) at a vertical mesh size of 10 m via an R script. The 80 and 100 m thickness values were based on representative thicknesses for Petermann ice islands in the literature (Halliday and others, Reference Halliday, King, Bobby, Copland and Mueller2012; Smith, Reference Smith2020).
To create uniform ram ice island meshes, a concave hull function (from the R package concaveman V1.0; using the default concavity setting) was used to identify nodes and their associated elements along the perimeter of the uppermost surface of the mesh. Mesh elements in this layer were then deleted along the periphery of the ice island iteratively inward to correspond with the specified ram length (Fig. 3). All nodes that were not referenced by the remaining elements were deleted from the mesh to guarantee a stable solution (LS DYNA, 2007). This resulted in an ice island mesh with the top of the synthesized ram occurring at approximately its natural waterline (with either 87.5% or 90% of the thickness submerged below the waterline).
Creation of a localized ram on ice island ‘IKY’, which calved in September 2011. The buffered surface extent of the ice island is delineated with a black outline. The synthesized uniform ram is represented by the element centroids (red points) of the polygon, while the ice island is represented by teal points. The calving edge was delineated by selecting two points along the ice island polygon corresponding to the region that subsequently broke off. An isolated ram can be created by computing a line (black dashed) between these points and deleting all ram elements on the larger portion of the ice island. This leaves a ram that extends only along the ice island calving edge.
To create an ice island with a full isolated ram, the coordinates at both ends of the calving edge were passed to an R function that deleted all mesh elements associated with a uniform ram along the perimeter of the ice island meshes except for the calving edge (Fig. 3), resulting in a ram fringing only the calving edge. The half extent variation of the isolated ram was created as above, but only considering the second and third quarter of the calving edge.
3.4. Segment set and height of the waterline
Gravity and hydrostatic pressure are predicted by the FEA model. The hydrostatic pressure equation requires specification of a segment set (basal plane) to indicate the surface of loading and the position of the waterline above the segment set (Sazidy, Reference Sazidy2020). Hydrostatic pressure (Phs) was calculated from the following equation:
\begin{equation}{P_{hs}} = {\rho _w}g{h_z}\end{equation}where ρw is the density of water, g is the acceleration due to gravity on Earth (9.8 m s−2) and hz is the ice island draft.
The segment set was assigned to each mesh using a custom function in R that returned the coordinates of nodes at the lowest mesh elevation. The position of the waterline above the bottom surface of the ice island was determined based on the ice volume, water density and ice density, and was listed in the reference plane field (ref-z) in the model’s hydrostatic curve section (Sazidy and others, Reference Sazidy, Crocker and Mueller2019). Since ice islands are very extensive and their thickness does not vary in the same way as it does for most iceberg shapes, we assume that their thickness is constant, and therefore (ignoring any rams) that their upper and lower surfaces are parallel to the water surface. During simulation, these surfaces will likely not remain parallel to the water due to ice bending or tilting caused by the presence of rams.
3.5. Material properties, forces and simulation settings
All of the required FEA parameters are defined and organized within the same LS-DYNA file that also contains the geometry of the meshes (LS DYNA, 2007; Table 1). The segment set (the plane at the bottom surface), position of the ice island waterline and damping coefficient (see Section 3.6) are parameters unique to each mesh and were derived independently for each ice island. In addition to the dimensional parameters mentioned above, parameters describing the material properties of the ice and sea water (density of ice (ρi), density of water (ρw), Young’s Modulus (E), Poisson’s Ratio (v) and the simulation parameters (termination time, timestep, tssfac) were transferred into the files. Termination time is the total length of the simulation run and timestep is the interval at which the calculations are made. tssfac is a time step scaling factor that further adjusts the time step according to the minimum element size. This is set to 0.9 by default, but our simulations ran with tssfac of 0.6 for improved accuracy (LS DYNA, 2007). The model assumes ice is an elastic material, and is compatible with the element erosion method, which deletes elements that exceed a stress threshold. The plastic nature of ice under high pressures was ignored since the ice thicknesses used in this study (80 and 100 m) were unlikely to cause significant creep deformation in ice islands. As well, Wagner and others (Reference Wagner2014) and Diemand and others (Reference Diemand, Nixon and Lever1987) reported realistic solutions with ice treated as an elastic material. The model calculates deformation in the ice at every time step and the ice island is allowed to move in response to imbalances between the forces of gravity and buoyancy. For more details on the static FEA model parameters and hydrostatic assumptions, refer to Sazidy and others (Reference Sazidy, Crocker and Mueller2019). Crack propagation in our model is simulated through element erosion, which was meant to visualize where stress exceeded a given threshold and to release stress in the surrounding ice. This method does not realistically represent the propagation of cracks, which would occur at size scales well below the mesh resolution we used, propagate over time-scales outside that of our short simulations, and involve complex interactions including the effect of hydrostatic pressure within the cracks (Jezek, Reference Jezek1984; Duddu and others, Reference Duddu, Bassis and Waisman2013; Ranganathan and others, Reference Ranganathan, Robel, Huth and Duddu2025). In our simulations, when elements eroded at the basal surface, progressive changes in stress caused propagation of eroded elements which suggest that rapid and complete failure would occur, and so was used as a visualization tool to identify the likely location of a fracture.
FEA model parameters
3.6. Damping coefficients
In oscillating physical systems, energy from resistive or frictional forces must be removed for the model to achieve realistic, smoothed results (Sazidy, Reference Sazidy2020). This is referred to as ‘damping’ and is accomplished with a coefficient term (Ds; unitless) whose value is related to the mass or volume of the object. Ds can be found by taking the time of two successive stress peaks in a time series plot of the element which reaches the mMPS, which is assumed to be in tension (Sazidy and others, Reference Sazidy, Crocker and Mueller2019):
\begin{equation}{D_s} = \frac{{4\pi }}{T}\end{equation}where T is the difference in time between successive oscillation peaks of mMPS in an undamped solution. Mosbeux and others (Reference Mosbeux, Wagner, Becker and Fricker2020) discuss in detail why simulated forces cause oscillations in deflection in this type of model, which is the cause of the oscillations we observed in the stress time series.
An empirical model was derived to predict damping coefficients of the meshes to avoid the need to fit this coefficient to all model permutations in this study (n = 312) by hand. Eight simplified rectangular meshes with volume ranges corresponding to those from the CI2D3 Database subset were created and solved. Damping coefficients were plotted as a function of mesh volume, and a logarithmic-linear model was fit to these points (R2 = 0.96, RMSE = 6.36). The model was used to predict an appropriate damping coefficient of each ice island model based on its mesh volume, and an R script wrote this value into the appropriate field of the damping card in the K file (Sazidy and others, Reference Sazidy, Crocker and Mueller2019; Smith, Reference Smith2020).
3.7. Sensitivity analysis
FEA resolution is typically optimized by selecting the coarsest resolution possible while maintaining adequate model skill (Reddy, Reference Reddy2005; Patil and Jeyakarthikeyan, Reference Patil and Jeyakarthikeyan2018). Since element resolution and aspect ratio are known to influence the accuracy of FEA results (Patil and Jeyakarthikeyan, Reference Patil and Jeyakarthikeyan2018), a sensitivity analysis was conducted to evaluate if the model was robust to a variety of horizontal and vertical mesh element sizes.
A simple idealized rectangular ice island (1040 × 480 × 80 m) with no ram was meshed at four different resolutions and each was run through the FEA model. A custom R function identified the specific element that reached the mMPS in each solution and calculated the mean MPS for all elements located within a 60 m radius of that element’s location. The element stress histories of each of these ‘peak stress zones’ were aggregated into a table and summary statistics were recorded to determine the influence of mesh size on modelled stress. This was used to evaluate the coarsest mesh size that could be used with minimal compromise of model skill (see Section 4.1 for details).
3.8. Compiling solutions
LS-DYNA was used to solve for MPS within ice islands under the prescribed ram configurations. This was done to examine the general stress pattern (spatial variation in MPS as it evolved over the several seconds of simulation time), to compute the mMPS and to model fracture. To simulate ice island fracture, the element erosion feature was activated, which deleted elements when their MPS exceeded the flexural failure threshold of ice island. This was set at 500 kPa in this study, the same threshold used by Wagner and others (Reference Wagner2014), allowing our results to be comparable to this earlier study. Failure in compression near the upper surface of the ice island is not considered here because of the relatively high failure threshold in comparison to tensile failure at its basal surface. The simulated fracture patterns were then compared with the CI2D3 Database polygons of the observed fracture locations (calving plots, e.g. Fig. 2a).
Other simulations were run with the element erosion feature de-activated to analyze the evolution and pattern of MPS regardless of whether the failure threshold was exceeded or not. These results were then used to fit an empirical model relating ice island geometry to mMPS.
Stress and fracture simulations within ice islands were generated for uniform rams at various ram lengths to determine if this ram configuration was a plausible cause of footloose fracture. For each of the 26 selected ice islands, MPS (and mMPS) was simulated for the 12 unique combinations of 2 ice island thicknesses (rt = 80 or 100 m), 3 ram lengths (rl = 20, 40 or 60 m) and 2 ram extents (full calving edge and half calving edge). Along with these model parameters the total ice volume along with the percentage of total ice volume contained within the ram were also computed from the ice island mesh.
Each K file was run through LS-Solver to determine the evolution of stresses within each mesh over time and examined in the LS-PrePost Graphical User Interface (GUI). Animations of the MPS for each mesh were viewed from different angles to observe the stress history and general behaviour of the ice island over time. The location and value of the highest stresses were recorded in a spreadsheet along with general notes about the simulation results.
3.9. Empirical modelling
It is assumed that an mMPS greater than the yield strength of ice will cause fracture and that a propagation of failing elements will lead to a calving event. A model that could predict mMPS may be useful as a tool to determine the stability of a given ram configuration. Such a model was postulated to help distil the important size and shape characteristics of ice islands that contribute to calving events and may help inform operational deterioration model development.
The results from the ice island meshes and stress simulations were used to statistically model the mMPS. Meshes from simulations that could not reach a peak stress within a computationally achievable simulation run time were discarded, as it was not possible to determine a representative mMPS of these ice islands. The dependent variable (mMPS) was extracted from the simulations, while all independent variables were extracted from the meshes (ice island volume, volume of the ram, ram extent and the percentage of total ice volume contained in the ram), and were then examined for normality using histograms, QQ-plots, symmetry (skewness) and the Shapiro–Wilk test. To satisfy parametric model assumptions of normality and equal variance, variables that were not normally distributed were transformed and then re-examined as above. The relationship between mMPS and each of the variables was assessed with a series of plots, linear models and correlation matrices. Five FEA model results with mMPS >1.7 MPa were removed from the dataset (an upper limit measured in empirical studies by Gagnon and Gammon (Reference Gagnon and Gammon1995)). Some explanatory variables were highly correlated with each other (i.e. collinear variables) and provided much of the same information. The collinear variables identified as having a Variable Inflation Factor (VIF) of 2 or more (Zuur and others, Reference Zuur, Ieno and Elphick2010) were removed from further analysis.
A Multiple Linear Regression (MLR) model was created with all available non-collinear explanatory variables and their first- and second-order interaction terms. Superfluous explanatory variables were removed from this maximal model to produce a model that maximizes explanatory power while retaining the fewest possible parameters (i.e. most parsimonious). Variables that were not statistically significant or caused the Akaike Information Criterion (AIC) of the model to increase (Crawley, Reference Crawley2007) were removed one at a time until the optimized model was determined. To estimate prediction ability in a larger population, bootstrapping was used. Model parameters were re-calculated using 10 randomly selected sub-samples from the original dataset. The Root Mean Square Error (RMSE), Mean Absolute Error (MAE) and the coefficient of determination (R2) for each of these model permutations were averaged. These bootstrapped skill metrics are expected to be lower than those in the optimized model and more representative of how the model would perform on new datasets.
3.10. FEA model validation
The FEA model used in this study was validated by Sazidy and others (Reference Sazidy, Crocker and Mueller2019) against an analytical solution of a semi-infinite beam on elastic foundation theory. The maximum stress and its location from the edge of the FEA ram and beam were found to be in very good agreement with each other (530 kPa vs 511 kPa, and 331 and 355 m, respectively). Further, Sazidy and others (Reference Sazidy, Crocker and Mueller2019) compared the analytical solution of a 2-D iceberg ‘Frances’ provided by Diemand and others (Reference Diemand, Nixon and Lever1987), to a 3-D version of Frances in the FEA model. Stress values were comparable, reaching 800 and 870 kPa, respectively, near the bottom centre of the iceberg.
To validate vertical deflection in the model, we re-examined the FEA model validation beam simulation of Sazidy and others (Reference Sazidy, Crocker and Mueller2019) to investigate vertical deflection along the horizontal axis at a steady-state time in the simulation (t = 150 s). The results were combined with the stress and location values in Table 5 along with the equivalent analytical solutions for deflection and solved using equation 19A from Hetényi (Reference Hetényi1946).
To further explore if smaller ice islands tilt to overcome ram-induced imbalances in hydrostatic forces while larger ice islands build up stress due to bending in the model, we created 3 simplified rectangular meshes with 40 m long rams. These were 80 m thick, 280 m wide and were 1960, 960 and 100 m-long, for the long, medium and short simplified meshes, respectively. We calculated deflection from horizontal, maximum stress, and the location of the maximum stress, to examine how much bending and/or tilt might occur in meshes at steady-state.
4. Results
The 26 ice islands retained for use in this study had a waterline length of 1.2 to 7 km and were between 0.37 and 14.7 km2 in surface area at the time of calving. Their locations ranged in latitude from 49°48′N to 82°09′N and in longitude from 47°28′W to 93°15′W, spanning thousands of kilometres across the Arctic and North Atlantic oceans. Figure 4 shows a map of their distribution throughout the eastern Canadian Arctic and North Atlantic regions, and Table 2 shows descriptive statistics of the calved fragment(s).
Map of the distribution of ice islands modelled in this study.
Ice islands from the CI2D3 Database examined in this study
The date is the last observation prior to calving. Unique three letter names listed here are extracted from the full ice island name in the database.
4.1. Sensitivity analysis results
Deviation between the resolution that was selected for use in this study (20 × 20 × 10 m) was not appreciably different from the finest modelled resolution of 10 × 10 × 10 m and was therefore justifiable for our FEA (the mean, maximum and minimum MPS were within 0.5% of each other). The RMSE of the element which reached the mMPS at a 20 × 20 × 10 m resolution when compared to the finest resolution (10 × 10 × 10 m) was 284 Pa, with a range of 0 to 147 kPa.
4.2. Uniform ram results
In ice islands with uniform rams, the stress and fracture patterns that propagated throughout the meshes during the simulations did not match those seen in the calving plots. In almost every case, the peak simulation stresses accumulated at the bottom centre of the meshes (Fig. 5a), and when the element erosion method was activated, cracks propagated outwards from the centre in different directions, causing the ice island to break into multiple pieces of similar sizes (Fig. 5b). This type of fracture is inconsistent with the footloose-type fracture patterns seen in the selected ice islands calving plots, which calved relatively smaller fragments at short distances (typically hundreds of metres) inwards from the ice island calving edge (Fig. 5d; Section 5.1).
Simulations of maximum principal stress (MPS; kPa) at the underside of ice island BSS with an assumed thickness of 80 m. (a) BSS with a 40 m uniform ram at simulation step 131 just before fracture. Note the contiguous ram around the entire perimeter of the ice island concentrates stress in the middle of the ice island. (b) The same simulation at step 207 following fracture. Tensile stress has been relieved by the erosion of fractures which would result in several large fragments of similar size. (c) BSS with a 40 m long ram isolated along the full extent of the calving edge (blue arrow here and in Figure 2) at simulation step 186. Erosion (fracture) was turned off for this simulation but the stress (which peaked later in the simulation at over 700 kPa) concentrates along where the fracture was observed in the CI2D3 Database (see Figure 2). (d) The same simulation at step 193 but with element erosion turned on. The fracture can be seen proceeding in the same direction as the simulation in (c) and joining the two ends of the calving edge. Compare this fracture to the observed calving event (Figure 2a). (e) Vertical displacement along a 1250 m cross section of ice island BSS simulated just prior to fracture. Displacement in cm is shown in colour with a maximum upward displacement of 17 cm at the ram. Fractures first appeared near the asterisk, approximately 268 m from the outer edge of the ram. The position of the cross section is shown with a dashed line in (d). Note that the displacement scaled by a factor of 50 relative to the scale of the cross section for visualization.
4.3. Full and half extent isolated ram results
When the element erosion feature was activated, the ice island meshes with isolated rams induced fractures that generally matched the calving patterns (similar fracture location) observed in the calving plots (Figs. 2a and 5d). These fracture patterns were consistent whether the rams extended along the full calving edge or half of it. However, greater ram lengths were needed before fractures occurred when the ram extent was reduced.
When the element erosion feature was deactivated, stress patterns similar to the previous cases were observed. Magnitudes of MPS greater than the 500 kPa threshold developed in regions corresponding to where the cracks developed in the element erosion simulations. Again, the magnitude of mMPS along these fracture planes were considerably higher with full-extent rams compared to those with half extent. Summary statistics for the results of both the full and half extent ram models are shown in Table 3.
Summary statistics for the percentage of ice islands that exceeded the 500 kPa stress threshold (upper table) and maximum of maximum principal stress (mMPS) range for those that did not exceed the threshold
n/a signifies not applicable.
4.4. Stress simulation results
When the element erosion method was de-activated, the mean mMPS were 473, 850 and 1104 kPa for 20, 40 and 60-m long rams with a thickness of 70 m, respectively. For the 90 m thick ice island rams, the mMPS values were 417, 706 and 945 kPa, for the respective ram lengths. This indicates that the stress increases considerably for every 20 m elongation of the ram regardless of thickness, and the effect is particularly large when the ram increases in length from 20 to 40 m. The thicker ice islands (100 m) had an average mMPS 12–17% lower than the thinner ice islands (80 m). Figure 6a shows that as ram length increases, the mMPS does as well. A linear regression between ram length and mMPS is significant (p < 0.001) with a slope of 14 kPa m−1. mMPS is inversely related to ice island thickness, with stress decreasing on average by 5 kPa m−1 of thickness (p < 0.05; Fig. 6b). These results suggest that in general, as the size of the ram increases relative to the ice island, internal stress will increase. However, this rule-of-thumb does not hold true for the smallest ice islands, since, in cases where the ram is already very large compared to the size of a small ice island, stress may be reduced due to tilting (Section 5.4).
Violin plots of the relationship between ram morphology and the maximum of maximum principal stress (mMPa). Mean and standard deviation are denoted by the circle and vertical lines. (a) Ram length versus mMPS. (b) Ram thickness versus mMPS. Stress increased with ram length and decreased slightly as thickness increased. The failure threshold used in some simulations is denoted by the red dashed line.
4.5. Empirical model results
All independent variables that were considered in the model (ram length, extent, thickness, volume, ram volume and the percent of volume contained in the ram) were positively skewed (skewness: 0.2–2.23). These were shifted towards normality as much as possible using various transformations, becoming the logarithm of ram volume (W: 0.99, p: 1.69 × 10−2) and ram extent (W: 0.98, p: 2.15. × 10−3), the square root of total ice volume (W: 0.94, p: 3.14 × 10−9), and the cube roots of mMPS (W: 0.97, p: 7.17 × 10−5) and the percent of volume in the ram (W: 0.98, p: 3.58 × 10−4).
The VIF test indicated that the total ice volume and volume of the ram were excessively multi-collinear with other variables and were eliminated from the model. The pruned multiple linear regression model used to predict mMPS indicated ram length (rl, p < 2.0 × 10−16), ram extent (re, p = 4.49 × 10−10), an interaction term between ram length and ram percent (r l · r p, p = 2.87 × 10−3) and an interaction term between ram extent and ram thickness (r e · rt, p = 1.74 × 10−5) were statistically significant variables. This model (Eq. 3) explained 55% (R2 = 0.55) of the variation in mMPS, with ram length alone explaining 68% of the modelled variation (Table 4). Model validation metrics extracted from the k-folds bootstrapped model were considered acceptable (R2 = 0.47, RMSE = 258 kPa, MAE = 202 kPa). mMPS can, therefore, be modelled as:
\begin{align}mMPS &= \big( 0.82{r_l} + 7.31log\left( {{r_e}} \right) - 0.15{r_l} \cdot \sqrt[3]{{{r_p}}} - 0.04log\left( {{r_e}} \right) \nonumber\\
&\quad \cdot \ {r_t} + 30 \big)^3\end{align}Final multiple regression model, indicating ice island predictor variables with the most explanatory power and statistical significance
The model coefficients (fit to the transformed variables), their standard error, p-value and the percentage of total model variation explained are listed.
where rl is the length of the ram (m), re is the ram extent (m), rt is the ram thickness (m) and rp is the percent of total volume contained within the ram. The final two negative terms represent the interactions between variables from the model between ram extent and thickness, which has a light damping effect on stress.
4.6. FEA model validation and deflection analysis
Table 5 shows the comparison between the FEA and analytical solutions for the validation beam and simplified mesh simulations. For the beam, we found very good agreement between maximum stress, location of maximum stress and the location where deflection begins relative to the horizontal axis, with moderate agreement of the vertical deflection values themselves. For the simplified meshes, agreement was moderate across all stress and deflection variables. Note that the analytical solutions we compared to are based on a semi-infinite beam, and so overall length of the ice beam does not affect stress in this model. Overall, the best match between numerical and analytical solutions was seen with the medium-sized simplified mesh, with an mMPS of 743 kPa versus 661 kPa at 670 m and 505 m away from the ram, respectively. Deflection values were 2.6 versus 1.1 m at 500 and 1000 m, respectively.
3-D numerical finite element analysis (FEA) model output compared to analytical solutions. The maximum principal stress value (σmax) and its location (distance from the ram (Lmax)), plus the maximum deflection from horizontal (yx) and the location (distance from the ram) where the deflection begins.
The validation beam, created by Sazidy and others (Reference Sazidy, Crocker and Mueller2019) is followed by three simplified meshes. The models generally agree, noting that length is not accounted for in the analytical method, and stress and deflection may converge as the width of the mesh is reduced.
* The analytical model did not output physical values for the short mesh, as it is not long enough. The numerical solution for this mesh indicated tilting (+12.8 m displacement at one end and −8.8 m at the other end) with no bending.
The smallest simplified mesh was not subject to any observable bending behaviour and instead deflected 12.8 m up on the ram side of the mesh and 8.8 m down on the opposite side. The dimensions of this mesh are consistent with that of a very large iceberg, and this tilting behaviour is typical of icebergs with irregular profiles.
The medium simplified mesh experienced some bending from the horizontal plane, while simultaneously exhibiting some tilt. However, it is likely an ice island of these dimensions would fracture before bending given the relatively high stress at its steady-state (743 kPa). The largest simplified mesh remained predominantly parallel to the water surface, with bending beginning about halfway across the mesh (760 m), and deflection up to 2 m at the ram edge. See Figures 7 and A1 for examples of deflection profiles and a comparison to analytical solutions. These results are also consistent with our interpretation of the empirical model, which indicates a negative relationship between mMPS and the size of rams relative to the total ice volume.
The influence of ice island size (perpendicular to the ram) on vertical displacement and stress caused by ram buoyancy. (a) A vertical cross section of a long mesh representing a simplified ice island (1960 m-long at the waterline with a 40 m ram, 280 m-wide and 80 m-thick), with an upward displacement of 38.5 cm at the ram. (b) Maximum principal stress (MPS) for the same mesh, with a maximum value of 651 kPa at the bottom, 290 m from the end of the ram. (c) Upward displacement in a medium-sized mesh (960 m-long at the waterline with a 40 m ram, 280 m-wide and 80 m-thick), which is 35 cm at the ram. (d) The MPS for the same mesh reaching 640 kPa at the bottom, 290 m from the end of the ram. (e) Vertical displacement in a short mesh (100 m-long at the waterline with a 40 m ram, 280 m-wide and 80 m-thick), where the entire mesh tilts by 0.2° (upward displacement of 76 cm at the ram, downward displacement of 27 cm at the other end of the mesh) instead of bending in response to uneven hydrostatic pressure. (f) The MPS for the same mesh reached only 63 kPa at the bottom, 70 m from the end of the ram, which is far lower than the other simplified meshes and in the middle of the bottom surface. All cross sections are plotted with a displacement scale factor of 50, which allows displacement to be visible at this scale. A 30 s simulation is shown. Note that the stress and displacement are not at equilibrium and erosion was not turned on (the long and medium meshes would have fractured, if that were the case). For maximum stress and deflection values, please see Table 5.
5. Discussion
Using FEA, this study successfully recreated calving events observed in the CI2D3 Database. However, without information about the keel shape, the presence of a ram and its actual dimensions remain unknown. In this sample of ice islands, assumptions were made regarding the ram shape, length and extent to determine if the footloose mechanism was a plausible cause of observed calving and to provide a framework for future 3-D modelling studies. These assumptions should be adjusted as more data about rams is collected, which should provide a better understanding of how and where they develop on ice islands. At present, the only ice island known to have calved with a ram was ‘PII-B-1’ in July 2012, with a ram length of 19 m measured just days before it calved (Wagner and others, Reference Wagner2014).
In spite of the paucity of quantitative ram data, the effect that ram buoyancy has on stress in simplified iceberg shapes has been modelled since the 1980s (Diemand and others, Reference Diemand, Nixon and Lever1987). However, the Sazidy and others (Reference Sazidy, Crocker and Mueller2019) model can predict stresses in 3-D which means that previously unexplored aspects of the footloose mechanism can now be analyzed for the first time, such as the stress distributions throughout the realistic ice island shapes, and the effect that ram extent has on the stress magnitude. The significance of these analyses are summarized below. Crawford and others (Reference Crawford, Crocker, Smith, Mueller and Wagner2024) found good agreement between our modelled maximum stress values and those they modelled using a 1-D beam model. As well, both stress and deflection from the FEA compare well against analytical solutions (Fig. 5; Sazidy and others (Reference Sazidy, Crocker and Mueller2019).
5.1. Ram extent and location
In the simulations, uniform rams surrounding the perimeter of ice islands caused high tensile stress to accumulate at the middle of the bottom surface for smaller ice islands (Figs. 5 and A2) and in a ring around the centre of larger ice islands (Fig. A3). When the element erosion method was activated, the smaller ice islands calved into several similar-sized fragments, while the larger one broke into several fragments surrounding the central area that was not under stress. These patterns were not observed in the calving plots of the selected CI2D3 Database ice islands, which suggests that uniform rams are not the cause of the local footloose-style calving that was the focus of this study. Ram formation may occur around the entire perimeter of ice islands, but localized footloose calving likely results from a preferential ram growth in specific locations. When the ice islands were fit with isolated rams, they caused stress patterns and fractures more consistent with those seen in the calving plots. The location where stresses were highest (and failure would have first occurred) up to several hundred metres from the calving edge was consistent with predictions from Wagner and others (Reference Wagner2014). When the element erosion method was activated, fractures propagated along the bottom surface of the ice island, consistent with where detachment presumably occurred, and extended in the upwards direction of the ice island.
5.2. Effect of a reduced ram extent
When ram extents were reduced by half, ice island simulations generally resulted in lower mMPS (except for those with 20 m long rams), with stress accumulating in an area still consistent with where detachment was observed to occur. Wagner and others (Reference Wagner2014) indicated that failure will occur when the ram length exceeds ∼60–80 m, for an 80 m-thick ice island and 70–80 m for a 100 m-thick ice island with a yield strength of 500 kPa. In this study, ice islands with a full ram extent failed with much shorter ram lengths (e.g. 20–40 m) regardless of thickness. However, when the ram extent was reduced by half, longer ram lengths were required to exceed the yield strength threshold, although the lengths are still far shorter than the ram lengths Wagner and others (Reference Wagner2014) predict cause calving. We suspect that, if the extent of the rams in our 3-D ice island models were reduced even further, the ram length required to cause calving would converge with the 1-D model solution from (Wagner and others, Reference Wagner2014). 1-D and 2-D models, which ignore ram extent, may under-predict stress, which is corroborated by our MLR, wherein the ram extent is the second-most important predictor of mMPS (see also Section 5.4), and also contributes significantly to the model with respect to its interaction with ram thickness. Therefore, we conclude that ram extent is a critical variable for modelling ice island stresses and should be considered in future deterioration models given that our simulations show a shorter ram length may be sufficient to induce calving if the ram extent is relatively large (Huth and others, Reference Huth, Adcroft and Sergienko2022).
5.3. Variations in yield strength
The skill of the FEA model relies on various physical parameters such as setting a single, fixed value for yield strength, which is reported to be highly variable in empirical studies (Robe (Reference Robe and Colbeck1980) and Vaughan (Reference Vaughan1995) report 0.5 MPa and Gagnon and Gammon (Reference Gagnon and Gammon1995) report 0.73–1.7 MPa). It should be considered that the ice samples selected for the yield strength testing are generally free of fractures and impurities which otherwise would result in considerably lower yield strength values (Wagner and others, Reference Wagner2014). The physical size of ice samples collected for the yield strength testing are also limited, and thus the scale effect of stress throughout ice masses as large as the ice islands are not well understood. A 500 kPa yield strength threshold was used in this study to examine the location and patterns of fracture in ice islands, based on support from previous calving studies (Robe, Reference Robe and Colbeck1980; Vaughan, Reference Vaughan1995; Wagner and others, Reference Wagner2014) and was selected for comparability with Wagner and others (Reference Wagner2014). Importantly, we also examined the model output without any ice erosion (i.e. unlimited yield strength), so that results could be considered in the context of different failure thresholds caused by, for example, crevasses and other physical flaws, and to create an empirical model that was unconstrained by our specific choice of failure threshold.
5.4. Empirical model parameters
Ram length was the single most important determinant of stress magnitude in the empirical model, accounting for 68% of the model variation in mMPS (Table 4). However, only three ram lengths were evaluated in this study and the effect of intermediate-sized ram lengths (e.g. 10 and 30 m) should be explored so as to better understand the relationship between ram length and stress.
Ram extent, the product of ram extent and ram thickness and the product of ram length and the percentage of total ice volume contained within the ram, were all statistically significant interactions with a discernible effect on stress magnitude (explaining 16%, 12% and 4% of total variation in mMPS, respectively). Previous work by Wagner and others (Reference Wagner2014) and Diemand and others (Reference Diemand, Nixon and Lever1987) did not explore ram extent in their analytical models, and state that failures tend to occur with longer ram lengths than we find here. Again, we attribute this to the 1-D nature of ice islands in their studies, which may require longer ram lengths to induce a similar amount of stress compared to similar ice islands that take into consideration the effect of ram extent. This is supported by the empirical model which shows extent to be the second biggest contributor of stress build-up, and is also suggested by the difference in stress magnitude in the simplified meshes, which reach higher stresses overall, compared to the relatively thin validation beam from Sazidy and others (Reference Sazidy, Crocker and Mueller2019).
The interaction between the percent of total ice volume contained in the ram and ram length is a significant model term which is negatively correlated with mMPS (Table 4). It seems counter-intuitive that an increase in ram length would decrease mMPS, but this interaction term also relates to the relative size of the ram to the ice island. We attribute this negative relationship to the tilting of small ice islands from the buoyancy of the ram instead of the build-up of bending stress in the keel that is seen in ice islands that are too large to tilt. Ram thickness-ram extent is likewise a significant negative interaction term which suggests that thicker and more extensive rams reduce overall stress build-up. Since thicker rams are also associated with increased ice island thickness, this may be attributed to the enhanced ice strength of thicker ice which may be important enough to offset larger ram extents.
The multiple regression model had a moderate explanatory power (R2 = 0.55), which suggests that the large-scale fracture prediction is possible with these parameterizations. Unexplained variance of mMPS (45%) may be accounted for by factors outside the scope of this study, but can be examined in greater detail in the future. Some of these factors may include highly irregular morphological variations like sharp corners, ice peninsulas or embayments seen in some meshes, which were noted to cause large spikes in mMPS compared to more symmetrically-shaped ice islands in simulations.
5.5. Bending and tilting of ice islands
It is worth noting when the element erosion is turned off and the simulation is run to a steady-state (∼150–200 s), the magnitude of stress and deflection, as well as their positions along the mesh move further away from the ram than if the simulations cease closer to the time fracture occurs. This highlights the importance of selecting realistic yield strengths for examining likely failure location and deflection values. The same set of simplified ice meshes summarized in Table 5 were plotted in cross section (Fig. 7a,c,e) alongside the MPS (Fig. 7b,d,f). All three simplified meshes showed upwards vertical displacement at the ram (Fig. 7a,c,e), however, the long and medium mesh bent upward at the ram, while the short mesh tilted in the water in response to the ram buoyancy, with no apparent bending. The MPS pattern matched the displacement where high tensile stress (>500 kPa) was found at 290 m from the end of the ram for the beams that bent, while stress was an order of magnitude lower for the small mesh.
These simplified simulations corroborate our interpretation of the negative interaction term in our empirical model and can be compared to 1-D model outputs from Wagner and others (Reference Wagner2014) along with 2-D simulations of other analogous calving front scenarios in ice shelves (Mosbeux and others, Reference Mosbeux, Wagner, Becker and Fricker2020). The long and medium mesh simulation took 11.5 and 13.7 s, respectively, to exceed the stress threshold of 500 kPa at 230 m from the edge of the ram (190 m from the waterline edge). If element erosion was turned on, this is where a fracture would have formed and stress would have been released. The distance between the waterline edge and the crack (∼190 m) is comparable to lw in Wagner and others (Reference Wagner2014). The deflection and MPS of the long and medium meshes followed the pattern observed in Mosbeux and others (Reference Mosbeux, Wagner, Becker and Fricker2020; their figure 8) with a rampart and moat at the upper surface and tensile stress concentrated at the underside of the ice within a few hundred metres of the edge of the ram.
5.6. Model limitations
Fractures simulated in the FEA model of ice islands with localized rams were mostly in good agreement with the location of fractures from the calving plots. However, when the element erosion method was activated the simulated fractures would not penetrate the topmost layer of elements which prevented a complete separation of fragments, and the ice island reaching a new state of hydrostatic equilibrium. This lack of complete separation is unphysical and was taken to be an FEA model artefact caused by the default integration scheme in the element type adopted in the model (LS DYNA element formulation I). Despite the failure of complete separation of fragments, the fact that the model reproduces realistic calving location implies that the synthesized ram dimensions and location may have been similar to those of the actual ram of these ice islands. We also focused part of our analysis on the location of the highest stress (mMPS), where the elements eroded and considered the ice island to have fractured completely along that plane when the mMPS exceeded the strength threshold. While our model was simple and computationally efficient, computational limitations caused some analyses to take a long time. This could be vastly improved with better computing resources, or by parallel processing some of the FEA computations.
5.7. Future work
This study identified some of the characteristics of rams that can induce the calving of an ice island. However, quantitative (or even qualitative) data of ice island ram dimensions are rare in the literature given the impracticality of obtaining them. In the absence of these, it is difficult to validate results from any footloose mechanism model, but we believe using a back study approach to infer possible underwater shapes and sizes of rams is a good start for exploring this method of deterioration at scale. In the future, we recommend field studies to obtain complete ice island dimensions including ram measurements through the use of multi-beam sonar to help validate our findings further and better understand their underwater shapes. This would reveal how and where rams develop along the keel of ice islands, it would establish a range of realistic ram dimensions, and could be used to validate FEA fracture models along with the stress thresholds they may use.
Though this study determined that thicker rams have a negative effect on mMPS, more variations in thickness should be considered in future models to provide a more comprehensive understanding of its effect on MPS. For example, with only two thickness values modelled, our empirical model assumed a simple linear relationship between the two values, but the relationship over a continuous range of thicknesses may be more complex. Thickness was an important variable in the (Wagner and others, Reference Wagner2014) model, and it is particularly important for modelling calving events of Antarctic ice islands, which can be over three to four times the thickness of Arctic ones (England and others, Reference England, Wagner and Eisenman2020).
This study also explored some of the effects of manipulating the ram extent but these were limited in scope. Future studies should examine how the variability of ram extent influences stress magnitude. For example, how calved fragments are affected by ram extents that are double the length of the calving edge, and if there is a minimum ram extent that can still produce a calving event.
This back analysis of ice island calving aimed to determine if specific ram shapes and sizes could have accounted for observed fracture patterns. Through this experimental approach, we inferred that the footloose mechanism was associated with specific calving events. Admittedly, this is not conclusive evidence and further work could be done to validate the model, especially with more field data as was mentioned above. To fully tune the FEA model, we recommend modelling ice island shapes that remained intact and were not observed to calve to ensure that the model does not predict high stresses under these scenarios.
The FEA relates solely to stress build-up based on ice morphology and buoyancy. In the future it may be feasible to combine this approach with models that predict other modes of ice island deterioration, such as surface ablation, edge wasting and wave notching. These types of deterioration models, which are driven with meteorological and oceanographic data (e.g. temperature, wind and current speed, etc.) may form a hybrid model to predict not only when a calving episode will occur but where on the ice island it will occur as well.
6. Conclusion
The logistics of in situ studies of ice islands are exceedingly difficult to manage due to the expense, remoteness and the planning that it takes to visit them. There are dangers inherent in working on or near ice islands since their breakup is unpredictable, making it challenging to observe this ephemeral process as it unfolds. The CI2D3 Database makes it relatively easy to analyze the drift and areal deterioration of many ice islands across time and space (Crawford, Crocker, and others, Reference Crawford, Crocker, Mueller, Desjardins, Saper and Carrieres2018) including for modelling melt (Crawford, Mueller, Desjardins, and others Reference Crawford, Mueller, Desjardins and Myers2018; Crawford, Mueller, and Joyal, Reference Crawford, Mueller and Joyal2018) and fracture (Zeinali-Torbati and others, Reference Zeinali-Torbati, Turnbull, Taylor and Mueller2021). This study determined plausible shapes and dimensions of hypothetical rams that could have resulted in some of these calving events. Since the underwater dimensions of these ice islands were unknown, ram shapes were synthesized using permutations of various ram lengths, extents and thicknesses to determine their influence on stress distribution and magnitude and calving location.
Our results indicated that a ram of 20 m is not typically long enough to cause calving in most cases but that a 40 m ram is. 60 m-long rams are likely to fracture and should rarely occur in nature assuming our strength threshold of 500 kPa is appropriate. Empirical relationships between ice island dimensions, ram dimensions and stress magnitude exist such that the highest stress in the ice island (mMPS) can be predicted moderately well. This empirical modelling and deflection validation approach highlights some of the important factors that should be considered in future examinations of realistic ice island shapes that fracture in 3-D (as opposed to idealized 1- and 2-D studies).
Wagner and others (Reference Wagner2014) found that for an ice island of 80 m in thickness, a ram length of approximately 60–70 m is needed to induce calving. However, using our FEA, nearly all ice islands (99%) calved with ram lengths ≤60 m, and some of these failed with ram lengths of 20 m. Since both studies adopted the same yield strength threshold, this could be attributed to the fact that (Wagner and others, Reference Wagner2014) used a 1-D model, which ignores the stresses that are associated with ice island planar shape and ram extent. Our study determined the extent of the ram is an important ancillary factor in ice island calving, which implies that 1-D models may overestimate the ram length needed to induce calving. If footloose-type calving can occur at shorter ram lengths than assumed up to now (Wagner and others, Reference Wagner2014), this suggests that the footloose process may be more common and, therefore, more important for ice island/iceberg deterioration than previously thought.
Analyses of large-scale deterioration calving events are still quite limited (Diemand and others, Reference Diemand, Nixon and Lever1987; Scambos and others, Reference Scambos, Sergienko, Sargent, MacAyeal and Fastook2005; Wagner and others, Reference Wagner2014; England and others, Reference England, Wagner and Eisenman2020; Zeinali-Torbati and others, Reference Zeinali-Torbati, Turnbull, Taylor and Mueller2021). Despite the exploratory nature of this study, it has identified some of the factors that promote calving events, and has provided a foundation for further investigation. Future research should be conducted with ram dimensions measured in situ to corroborate our findings and to learn more about the footloose calving process. Ice island shape should also be investigated as a factor in footloose-style calving events as rams simulated along an irregularly-shaped ice island edge in this subset seemed to induce higher stresses.
Climate warming is projected to reduce Arctic sea ice extent (Overland and Wang, Reference Overland and Wang2013), but may also increase the production of ice fragments in the future (Bevis and others, Reference Bevis2019). A reduction in sea ice may lead to increased marine activity, which may result in more interactions with glacial ice hazards in the coming years (Dalton and others, Reference Dalton2025). A better understanding of the processes involved in ice island deterioration will enhance operational risk management (Fuglem and Jordaan, Reference Fuglem, Jordaan, Copland and Mueller2017). Our research into the process of footloose calving can contribute to the development of an improved deterioration model that captures episodic calving in a robust way (e.g. Crawford and others, Reference Crawford, Crocker, Smith, Mueller and Wagner2024). This can, in turn, be used to improve drift trajectory models since they require waterline length estimates and often incorporate deterioration models to simulate the evolution of this key parameter over time (Kubat and others, Reference Kubat, Savage and Carrieres2005).
Competing interests
The authors declare none.
Acknowledgements
This research was made possible by funding from an Early Career Researcher Award (Ontario Ministry of Research and Innovation), ArcticNet, Network of Centres of Excellence of Canada, Canada Foundation for Innovation, the Ontario Research Fund, Northern Scientific Training Program, Amundsen Science and the Ina Hutchison Award in Geography at the Department of Geography and Environmental Studies (Carleton University). We thank Tetiana Bukhinska and Julia Albert for the ice island artwork. Finally, we thank Dr. Till Wagner and Dr. Doug Benn for their reviews of a prior version of this manuscript.
Appendix A
Various examples of deflection profiles along the horizontal axis of (a) the validation beam first examined by Sazidy and others (2019) and (b) our additional simplified ice meshes as described in Section 4.6. Note that bending is most obvious in the largest ice island mesh, with the other two exhibiting different degrees of tilting behaviour. Tilting is very pronounced in the smallest mesh that has dimensions consistent with that of a large iceberg, not an ice island.
Simulations of Maximum Principal Stress (MPS; kPa) at the underside of ice island PJJ with an assumed thickness of 80 m. (a) PJJ with a 40 m uniform ram at simulation step 201(20 s) with erosion turned off. The contiguous ram around the entire perimeter of the ice island concentrates stress in the middle of the ice island as in Figure 7. (b) The same simulation at step with element erosion turned on, simulating the fracture of the ice island into seven fragments of similar size. (c) PJJ with a 40 m long ram isolated along the full extent of the calving edge (blue arrow) at simulation step 201. Erosion (fracture) was turned off for this simulation showing an accumulated stress of 633 kPa concentrated along where the fracture was observed in the CI2D3 database. (d) The same simulation at the same step but with element erosion turned on creating a fragment that was similar to what we observed in the CI2D3 Database. (e) Vertical displacement along a 1576 m cross section of ice island PJJ simulated just prior to fracture at step 114 (11.3 s). Displacement in cm is shown in colour with a maximum upward displacement of 27 cm at the ram. Fractures first appeared near the asterisk, approximately 340 from the outer edge of the ram. The position of the cross section is shown with a dashed line in (d). Note that the displacement scaled by a factor of 50 relative to the scale of the cross section for visualization.
Simulations of Maximum Principal Stress (MPS; kPa) at the underside of ice island CHG with an assumed thickness of 80 m. (a) CHG with a 40 m uniform ram at simulation step 152 (15 s) with element erosion turned off. The contiguous ram around the entire perimeter of the ice island concentrates stress in a ring around the middle of the ice island but the middle of the ice island is not under stress. (b) The same simulation at step 152 with element erosion turned on showing the formation of fragments surrounding the central area where stress was low. Since this pattern of calving has not been observed in the CI2D3 Database, it is unlikely that rams form evenly around ice islands in nature. (c) CHG with a 40 m long ram isolated along the full extent of the calving edge (blue arrow) at simulation step 202 (20 s). Element erosion was turned off for this simulation but a stress of 548 kPa concentrates along where the fracture was observed in the CI2D3 database and in (d). (d) The same simulation at step 202 but with element erosion turned on with a fracture aligned with the stress in (c). (e) Vertical displacement along a 2333 m cross section of ice island CHG simulated just prior to fracture at step 154 (15.4 s). Displacement in cm is shown in colour with a maximum upward displacement of 24 cm at the ram. Fractures first appeared near the asterisk, approximately 180 from the outer edge of the ram. The position of the cross section is shown with a dashed line in (d). Note that the displacement scaled by a factor of 50 relative to the scale of the cross section for visualization.
Open access
