A statistical fracture model for Antarctic ice shelves and glaciers

Antarctica and Greenland hold enough ice to raise sea level by more than 65 m if both ice sheets were to melt completely. Predicting future ice sheet mass balance depends on our ability to model these ice sheets, which is limited by our current understanding of several key physical processes, such as iceberg calving. Large-scale ice flow models either ignore this process or represent it crudely. To model fractured zones, an important component of many calving models, continuum damage mechanics as well as linear fracture mechanics are commonly used. However, these methods have a large number of uncertainties when applied across the entire Antarctic continent because the models were typically tuned to match processes seen on particular ice shelves. Here we present an alternative, statistics-based method to model the most probable zones of the location of fractures and demonstrate our approach on all main ice shelf regions in Antarctica, including the Antarctic Peninsula. We can predict the location of observed fractures with an average success rate of 84 % for grounded ice and 61 % for floating ice and a mean overestimation error rate of 26 % and 20 %, respectively. We found that Antarctic ice shelves can be classified into groups based on the factors that control fracture location.

1 Introduction

In recent years, increased positive-temperature anomalies have been observed in Antarctica (Jansen et al., 2007; Johanson and Fu, 2007; Steig et al., 2009; Vaughan et al., 2003), and future climate change in this area may be even more pronounced (Vaughan et al., 2003). This may cause the state of the Antarctic Ice Sheet to change significantly and could lead to a release of fresh water currently stored in the ice sheet; West Antarctica alone could contribute up to 4.3 m to global sea level (Fretwell et al., 2013). Thus, understanding the factors that control the mass balance of the Antarctic Ice Sheet is crucial if we want to better understand the future impact of climate change and contribution of Antarctic ice mass loss to global sea level rise (SLR). Increased calving from the major ice shelves between 1998 and 2003 led to growing concern about the ice sheet stability (Shepherd et al., 2012).

Overall, research on crevasse propagation started as early as 1955, and calving parameterization has been under development for the last 20 years. It has been shown that the increased ice mass loss from Antarctica is caused primarily by an increased number of calving events in the last 2 decades, which has led to significant ice front retreat (e.g. the collapse of the Larsen B Ice Shelf (IS) and break-off of a part of the Larsen C Ice Shelf; Borstad et al., 2013; Gudmundsson, 2013; Hogg and Gudmundsson, 2017; Jacobs et al., 1992; Katz and Worster, 2010; Mercer, 1978). A number of studies have shown that increased calving can lead to destabilization of ice shelves and thus to a loss of the supporting mechanism (known as the “buttressing effect” or “back stress”) they provide to inland ice in Antarctica (Borstad et al., 2013; De Angelis and Skvarca, 2003; Dupont and Alley, 2005; Goldberg et al., 2009; Gudmundsson, 2013; Jezek, 1984; Katz and Worster, 2010) . This support can be crucial for the overall stability of the West Antarctic Ice Sheet as strong basal melting and reduced ice shelf buttressing can make the ice sheet unstable (Miles et al., 2013).

Developing a reliable calving law requires knowledge of where these fractures are located and how they evolve. Fractures in the Antarctic Ice Sheet and ice shelves are visible in satellite imagery and can occur more often than every 50 m. Because of the size of Antarctica, the only feasible means of creating a database of fractured zones is through the analysis of satellite imagery and altimetry. However, fractures can be covered by snow and/or not be visible because of poor resolution of the available imagery. It is for these reasons that inverse methods are often used (Borstad et al., 2013, 2016). Because of the incomplete knowledge of the location of fractured zones in Antarctica, models for representing fracturing in ice shelves/glaciers remain poorly constrained. The modelling of initiation and propagation of fractures is an active area of research; however no models to date have been able to successfully model this process on an ice sheet scale. Existing studies have focused on Greenland (Cook et al., 2014; Krug et al., 2014; Nick et al., 2010; Sugiyama et al., 2015), Svalbard (Chapuis and Tetzlaff, 2014) and the Antarctic Peninsula (Bassis and Walker, 2012; Otero et al., 2010), and a method that can universally describe calving at any ice shelf in Antarctica has not yet been developed.

The aim of this study is to construct an empirical model that can predict the locations of fractures. We focus on modelling of crevasses (surface fractures less than 200 m wide) on the surface of the Antarctic Ice Sheet and surrounding ice shelves. Our model for predicting fractured zones is based on a probabilistic approach, where we utilize a logistic regression algorithm (LRA) to find a relationship that enables the prediction of fracture locations. Our approach accounts for many potential parameters – including geometry, mechanical properties and flow regime (predictor parameters) – and is based on a combination of modelling and remote sensing. We use a data set of fracture observations, built by careful manual selection of the locations of visible fractures in satellite images, to build a model that can identify fractured regions on most of the Antarctic ice shelves as well as grounded ice regions around ice shelves in Antarctica, including the Antarctic Peninsula. We compare the ability of our model to match observations of fractures from satellite imagery versus the predictive ability of the damage-based method of Borstad et al. (2013).

In Sect. 2 we describe previous studies focused on fracture and calving modelling. In Sect. 3 we provide information about the data sets used to construct our model. In Sect. 4 we construct our probability model. In Sect. 5 we group glaciers/ice shelves based on common characteristics; then in Sect. 6 we describe each group. Our best results were achieved combining the LRA with Bayesian as well as Jensen–Shannon divergence (JSD) theory described in Sect. 4.

Nye (1955)Weertman (1973)Smith (1976)Vaughan (1993)Motyka (1997)van der Veen (1998a, b)Pralong et al. (2003)Pralong and Funk (2005)Benn et al. (2007a, b)Nick et al. (2010); Otero et al. (2010)Borstad et al. (2012)Levermann et al. (2012)Duddu and Waisman (2012)Duddu and Waisman (2013)Bassis and Jacobs (2013)Åström et al. (2013)Nick et al. (2013)Cook et al. (2014)Albrecht and Levermann (2014)Krug et al. (2014)Morlighem et al. (2016)

Table 1 Development of calving parameterizations.

2 Background

2.1 Current state of calving computations in ice sheet models

A number of calving parameterizations have been developed and implemented in some software packages, but none of them includes the propagation of fractures both vertically and horizontally. Most of the available parameterizations are specific to a particular case and set of predictors (Pralong and Funk, 2005), or calibrated for a particular location (Krug et al., 2014), and therefore cannot be applied generically to any ice shelf.

The majority of calving models are one- or two-dimensional and are built on simplified physics (Alley et al., 2008) using either stress-related criteria or parameterizations based on crevasse penetration depth (Benn et al., 2017). The Community Ice Sheet Model (CISM) assumes that calving takes place when the water depth reaches a certain value (Price et al., 2014). This water depth calving model uses flotation criteria to estimate the location of the glacier terminus. It allows calving to be linked to glacier dynamics as well as surface melting when applied to marine-terminating glaciers in Greenland (Nick et al., 2010). However, it cannot describe calving at floating ice shelves in Antarctica since the floating part is simply removed from the CISM model (the water depth relationship requires a glacier to calve once it floats). The Parallel Ice Sheet Model (PISM) uses a calving algorithm based on along- and across-flow strain rate (Levermann et al., 2012), which is based on the correlation between the calving rate and the first-order approximation of local ice flow spreading rate (it includes spreading rates of the first order and assigns all higher-order terms to zero). This idea is based on the observations of the increase of calving rate with along-flow ice shelf spreading rates and the spreading rates perpendicular to the calving front. However, it considers only large-scale behaviour and does not take into account the formation and propagation of crevasses. A second method implemented in PISM involves calculating a calving rate based on the critical ice thickness, which is mainly used to model calving of marine-terminating glaciers rather than floating ice shelves (due to different physics governing calving between the grounded ice and floating ice). More explicit simulations have been performed using discrete element models describing short-term calving events (Åström et al., 2013).

A number of other approaches such as calving laws proposed by Pralong and Funk (2005) and Duddu and Waisman (2012) have advanced our understanding of calving, and the main existing studies are presented in Table 1. These methods can include hydrofracture and other modes of failure but have largely been applied to grounded calving margins, in contrast to the methods by Borstad et al. (2016) that are calibrated to remote sensing data and have been applied to ice shelves but not grounded calving margins. Moreover, most of the mentioned methods might not be applicable in a generalized large-scale case.

To date, linear elastic fracture mechanics (LEFM) (van der Veen, 1998a) and continuum damage mechanics (CDM) (Kachanov, 1958) are the most common methods used to model fractured zones, which is important information for modelling calving. A method based on CDM proposed by Borstad et al. (2016) can describe both the formation and evolution of fractures in a fully viscous continuum damage model, although coupling with an elastic damage model might be more appropriate for representing short-timescale evolution of fractures. Krug et al. (2014) built an alternative scheme by combining LEFM and CDM and found that they could match the observed evolution of a tidewater glacier in Greenland. This method is more complex compared to earlier approaches as it allows for both viscous and elastic behaviour and is able to reproduce development of small crevasses over a long period of time. The LEFM-based approach consists of calculating a stress intensity factor around fractures and assuming that they propagate until the factor falls below a certain critical value. To apply this method to any ice shelf, the modelled fractured zones need to be in good agreement with the observed surface fractures. Therefore, modelling the location of the fractured zones is an important basis for the subsequent estimation of fracture depth as well as calving, and it must be described in the ice sheet models accurately. The Elmer/Ice model (Gagliardini et al., 2013) combines CDM and LEFM to model calving; however, most of the experiments with Elmer/Ice calving has only been applied to Greenland glaciers (Krug et al., 2014), which have a different calving mechanism from the floating ice shelves in Antarctica (van der Veen, 2002). The Ice Sheet System Model (ISSM) also has a calving algorithm only for marine-terminating glaciers, which is derived from a tensile von Mises yield criterion (Morlighem et al., 2016).

The continuum damage mechanics approach, based on the method suggested by Kachanov (1958), includes estimation of damaged zones where the ice is softened due to the presence of fractures. The continuum damage mechanics approach has been successfully applied to a wide range of engineering problems and to model damage at individually selected ice shelves such as the Ross, Filchner-Ronne, Amery (Bassis and Ma, 2015), Larsen C (Borstad et al., 2013) and Larsen B (Borstad et al., 2016) ice shelves. Damage is a scalar variable used to determine failure of ice and the nucleation of fractures, usually when the damage predictor reaches a certain value. There are two different methods for calculating damage: methods applied to invert for damage and methods used to model damage advection in ice sheet models. Borstad et al. (2012) suggested a direct inversion for damage using a cost function. Later, Borstad et al. (2013) proposed a method to calculate damage as a post-processing routine after inverting for the ice viscosity. Inversion of ice rheology is performed following Larour et al. (2005); then damage is calculated from the inversion of velocities at the ice shelves, which is based on minimizing the cost function that quantifies the misfit between the observed and modelled surface velocities. A method for modelling of advection of damage was suggested by Krug et al. (2014) and Albrecht and Levermann (2014) using an advection scheme and a source function. The main limitations of this method are the choice of what should be used as a source function as well as the number of decisive parameters that define the damage evolution (damage threshold, initiation threshold and the enhancement factor). The source function is the controlling factor in the damage advection, and Pralong et al. (2003) as well as Krug et al. (2014) proposed a source function definition. Both of the approaches work well for particular regions, but control predictors in this model have been derived using data from only one specific glacier in Greenland (Pralong and Funk, 2005) and through the use of small-scale laboratory experiments (Duddu and Waisman, 2012) and have not yet been generalized to be applicable to all ice shelves/glaciers. Another limitation of these models is that they do not yet account for factors such as ice fabric and impurities (Borstad et al., 2012).

Recently, Borstad et al. (2016) proposed a framework where, instead of computing a damage source term as is usually done, damage is part of a generalized constitutive relationship. This approach has a number of advantages as it allows the calculation of mechanical ice weakening and the prediction of the degradation of ice shelves. This can significantly improve the accuracy of identifying zones where the ice is weakened by illuminating the uncertainties related to the source function. The main weakness of the approach lies in determining the constant parameters that define damage, because the validity of the parameter values can only be tested when an ice shelf undergoes pronounced mechanical changes, as did the Larsen B Ice Shelf in 2002. The parameters suggested by Borstad et al. (2016) have not been tested for other ice shelves apart from Larsen B, and so it is unknown whether the approach is valid for other locations and settings.

3 Construction of the input data to the statistical model

Our statistical model is built upon knowledge of the flow and geometry parameters collected from 35 ice shelf regions in Antarctica (see Fig. 3), each including both ice shelves and the grounded ice around 100 km upstream from the grounding line (hereafter referred to as ice shelf regions or ice shelf/glacier).

To construct the model, we needed two types of data: observations of fractures and the flow/geometry information in each fracture. To run the model for a particular ice shelf region, only information about the flow/geometry is required. We chose 45 ice shelves/glaciers where high-resolution imagery was available. For calibration of the model we constructed a data set containing information about flow regime/geometry and the locations of some observed fracture or lack thereof (described in Sect. 3.2) collected from 35 ice shelf regions. We use these observations to derive the β values in the LRA. The 35 ice shelves/glaciers were chosen at random out of the 45 analysed regions for which we had observed fractures. To validate the model results, we use an independent set of observed fractures on these 35 ice shelves/glaciers. Thus, there is no duplication of use of observations in the calibration and validation processes. Overall, there are 20 times more observations of fractures used for validation of the model than in calibrating the model. We include in the validation process fracture observations on an additional 10 ice shelves/glaciers, which provide completely independent observations against which to compare our model results.

Data about geometry and velocity used in our statistical model were taken from observations, while other flow and geometry parameters were calculated either using ISSM or independently (described in Sect. 4.2). The 450 m resolution horizontal ice velocities were taken from interferometric synthetic aperture radar (InSAR; Rignot et al., 2011a, b). The data for the geometry of the ice shelves and surrounded grounded ice (bedrock topography, ice thickness and glacier surface) were interpolated from Bedmap2 data (Fretwell et al., 2013) at 1 km spatial resolution. To calculate stresses, strain rates, back stresses, friction coefficient and viscosity of the ice, we used ISSM (Larour et al., 2012).

3.1 Ice Sheet System Model setup

ISSM is a fully dynamic model that includes both two-dimensional (2-D) and three-dimensional (3-D) stress balance approximations. Our experiments rely on the shelfy-stream approximation (SSA) as it is computationally cheap and suitable for modelling floating ice shelves and grounded ice streams undergoing widespread basal sliding. We ran one simulation to create a stress balance solution per region (ice shelf/glacier), which allowed us to obtain the predictor parameters required for the calculation of the probability of fracturing. We used SeaRISE air temperature, snow accumulation and geothermal heat flux (Le Brocq et al., 2010) as climate forcing data. The information about the ice temperature for grounded ice is calculated as the depth-averaged temperature from Van Liefferinge and Pattyn (2013) and Pattyn (2010) (21 vertical levels). The steady-state depth-averaged ice temperature on floating ice shelves (mainly used for the calculation of damage) was calculated using surface, basal temperatures and basal melting rate according to Holland and Jenkins (1999). To calculate ice temperature, we corrected the surface temperature with a lapse rate and imposed it on the ice surface. Basal melting rates on ice shelves were taken from Depoorter et al. (2013). Basal friction under grounded ice and rheology for floating ice were calculated from an inversion of velocities (Khazendar et al., 2007), and the sliding law is the Budd sliding law (Budd et al., 1979). In the inversions we used regularization to penalize sharp gradients of the cost function, calibrated using an L-curve analysis (Morlighem et al., 2013). We set boundary conditions as follows: the upper surface is considered stress-free, and friction is applied at the ice–bedrock interface. At the inflow boundary we applied Dirichlet conditions, and at the ice front we applied Neuman boundary conditions. The position of the grounding line position is calculated using a flotation criterion. We used the outputs of the simulations for each of the 45 regions to evaluate the predictor factors in our probabilistic method (see Sect. 4.1).

For all the simulations, we used a multi-resolution mesh approach for the chosen domains in East and West Antarctica as well as the Antarctic Peninsula. This method was chosen due to the fact that, on the one hand, using a 50 or 100 m mesh resolution created a significant increase in the computational time of the model but that, on the other hand, it was important to have a fine-resolution mesh in order to model surface fractures, as the distance between them is normally around 50–100 m. In order to have a fine resolution together with smaller computational time, we first calculated all the main predictors on a 200 m resolution mesh (to achieve a faster computational speed) and then interpolated the values to nodes on a 100 m resolution mesh (to use in our fracture model resolved at 100 m). All further computations and analyses were performed on this finer mesh.

3.2 Observing fractures using satellite images

We focus only on predicting the location of surface crevasses without modelling rifts, since the processes that cause rift opening might differ from processes that allow surface crevasses to stay open. In fact rifts might be formed due to presence of basal fractures, tidal deformation (Bromirski et al., 2010), ocean swell (Bassis and Walker, 2012) and stay open due to presence of mélange (a mixture of icebergs and sea ice) (Bassis et al., 2005; Fricker et al., 2005; Larour et al., 2004; Rignot and MacAyeal, 1998). To model this, we would need to include information about ocean temperature, sea ice and seismic activity, all of which is outside the scope of this paper. Moreover, rifts form when cracks propagate through the entire ice thickness and, therefore, the ice becomes effectively discontinuous. We, therefore, do not include rifts and focus only on surface crevasses.

In order to obtain the observations of the location of fractures on the ice sheet surface we used satellite images taken from Google Earth Pro, where images of the Antarctic Ice Sheet were available at different spatial resolutions. However, to be able to see surface fractures, we limited our choice to only images with a horizontal resolution smaller than 10 m for the period between 2011 and 2015. We included only regions with at least one high-resolution satellite image and where it was relatively easy to identify surface fractures.

The visual images of the ice surface include many features, and it is important to distinguish the surface fractures from other patterns such as surface troughs (formed as a result of presence of bottom crevasses or subglacial channels). It has been suggested by Luckman et al. (2012) that wide features in the images with a large spacing between them (e.g. >1 km) are more likely to be troughs associated with basal crevasses. Modelling of basal fractures is outside the scope of this paper. Thus, such large linear features that are visible in the satellite images are surface troughs and should not be interpreted as surface fractures that our model fails to predict.

To construct a set of observed fractures, we manually selected fractured locations as well as non-fractured zones that we could identify in the satellite images. Most of the identified non-fractured regions are located in blue-ice regions, which are areas with low snow accumulation or where the snow has been removed by the wind. In such areas we can clearly see where the ice is not damaged. It is important to note that in some locations the resolution of images was not always sufficiently high to clearly see every fracture. Moreover, some surface fractures may be covered in snow and, therefore, are not identified by our analysis.

For 35 ice shelf regions we constructed two different types of data sets: “calibration” and “evaluation”, for building the statistical model and for studying the output of the model, respectively. For the other 10 ice shelves we only constructed the evaluation data set as we did not use these 10 regions to construct our model. In the calibration data set we select a subset of observed fractures, being a representative sample of locations where fractures are found on 35 ice shelves/glaciers. The statistical approach requires training in a large number of ice shelf regions with different characteristics and a variety of observations. Therefore, we used the calibration data set to build the LRA model. This improves the reliability of the model, as the diversity in sampling provides a better estimation of correlation coefficients (called β coefficients in LRA). We assign a value of 1 to fractured nodes and 0 to non-fractured nodes (due to the fact that the LRA that we use in our approach, described later in Sect. 4.1, uses a categorical input).

We form the evaluation data set to test how well our new approach predicts fractures for each ice shelf region individually. The evaluation set for each glacier/ice shelf is much larger than the number of fractures selected from each of the regions for calibration as we did not need every observed fracture to construct the model (as previously mentioned, it was the variety and not the number of data points that was required for a successful construction of our model). Although we did not need to select all the fractures on the ice sheet surface to build the calibration data set, to construct the evaluation data set we made a concerted effort to select the majority of the visible surface fractures in each of the ice shelf regions. It is possible that some fractures were missed due to the large spatial extent of the experiments. Moreover, we do not present every fracture in the figures in this paper in order to make the figures legible. In addition, we perform validation experiments with another 10 ice shelf regions to test how well the LRA works for a randomly selected ice shelf/glacier that was not a part of the construction of the model.

It is important to note that the evaluation data sets are not just discrete values (0 and 1) but are rather a continuous field representing the probability of observing a fracture in a location. In a node where we could see a fracture, we assigned the probability of observing a fracture to 1. Nodes around the observed fracture are more likely than not to be fractured. It is important to mention that the spacing of crevasses is often linked to their depth. A single crevasse can penetrate much deeper than a crevasse in a set of closely spaced crevasses. However, in this study we do not focus on estimating either depth or spacing of crevasses. Therefore, we then set the probability of observing a fracture to simply decrease from 1 to 0.55, decreasing with increased distance from the observed fracture (within 500 m radius). On the other hand, when a non-fractured node was found within a region with high-resolution imagery, we assigned the probability of fracturing in this node to 0.05. Within a 500 m radius of the non-fractured node we allowed the probability to increase linearly from 0.05 to 0.4. In all other nodes we set the probability of observing a fracture to 0.5. The last assumption is due to the fact that if there are no fractures visible in the area of poor resolution of the image it is equally likely for the node to be fractured or non-fractured. This allows us to account for uncertainties of the observations, since it is not always possible to determine whether there are no fractures or whether fractures are just not visible. We do not include any information about the depth and spacing of the crevasses.

4 Methods

We used statistics-based methods as an alternative to physics-based approaches in order to gain insights into the location of fractured zones in ice shelves and glaciers. In the well-known damage-based approach, the damage variable varies from 0 to 1, representing the fraction of a volume that is fractured, with 0 being not fractured and 1 being fully fractured. Instead of using the damage-based method, we use the LRA, which provides us with the probability of fracturing (also varying from 0 to 1). We then apply this method to derive fracture likelihood functions for both floating ice shelves and the grounded ice for any ice shelf region. To construct the likelihood function, we need to find coefficients that describe the relationships between predictor variables and what we want to predict (in our case it is surface fractures not including rifts). Thus, in order to create a statistical model, we use our calibration data set of observations of surface fractures and non-fractures as well as information about the flow regime at the locations of each observation (predictor parameters).

Our main goal is to determine the most likely location of surface fractures. We do not focus on identifying the location of their initiation, since it is not possible to know whether the observed fractures were formed where observed or advected to that position after having formed upstream. We tried to select observed fractures where there were no other fractures visible upstream, meaning that the observed fractures would identify the initiation zones, but this may not always be possible. The model will, therefore, predict the locations not only of initiation of fractures but also of some zones to which fractures have advected. For this reason, we do not distinguish between the high-advection (advection from upstream) and low-advection (because of local stresses) cycles (Colgan et al., 2016). Although we do not directly model advection, the statistical model predicts the presence of both initiated and advected fractures without distinguishing one from the other. The question that arises then is, how do we know that the flow regime conditions that caused opening of the fracture are the conditions at the point where the fracture is observed and not the conditions upstream from the observed fracture (in the case of advection)? However, even if an observed fracture was not formed at a particular location but was advected with the ice flow, it is still visible in the satellite image. The fact that fractures can be seen indicates that there are factors at that location that act to permit the fractures to exist, whether they formed in that particular location or remained open after being advected from upstream (since another combination of factors could close the fracture).

This section is structured in the following order. First, we present our method (logistic regression algorithm) used for predicting the formation of fractures. Second, we describe the predictor factors (predictors) we include in this method. Then, two methods used for optimizing a set of predictor factors are presented (Bayesian-based algorithm and Jensen–Shannon divergence). Finally, we present the damage calculation used for a qualitative comparison with our results.

Table 2 Predictor factors (predictors).

4.1 The logistic regression algorithm

Logistic regression is a statistical technique generally used to classify data based on values of input fields. The method is similar to linear regression but takes a categorical target field (in our case nodes which are fractured or non-fractured) instead of a numerical series. The logistic function allows us to calculate the likelihood of an event as a function of different predictor factors (see Table 2 for the predictors used in our model). Taking any range of data, it produces values from 0 to 1, and thus it can be used to represent the probability of fracturing (Hosmer Jr. and Lemeshow, 2004).

To apply the logistic regression algorithm, we constructed a logistic function P_j (Eq. 1) that describes the probability of a certain node being fractured as a function of the predictor factors x_i. This function is not designed to provide any information about the depth of a fracture, just its spatial location.

$$\begin{array}{ll}\text{(1)}& {\displaystyle }& {\displaystyle }{P}_j=\mathrm{Prob}(X=\mathrm{1}|x_{ij})={\displaystyle }& {\displaystyle \frac{\exp ({\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}\cdot x_{\mathrm{1}j}+{\mathit{\beta }}_{\mathrm{2}}\cdot x_{\mathrm{2}j}+{\mathit{\beta }}_{\mathrm{3}}\cdot x_{\mathrm{3}j}+\mathrm{\dots })}{\mathrm{1}+\exp ({\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}\cdot x_{\mathrm{1}j}+{\mathit{\beta }}_{\mathrm{2}}\cdot x_{\mathrm{2}j}+{\mathit{\beta }}_{\mathrm{3}}\cdot x_{\mathrm{3}j}+\mathrm{\dots })}},\end{array}$$

where j is the node number, x_ij is the value of predictor x_i for node j and β are correlation coefficients.

The unknown coefficients β_i are found by maximizing the likelihood function L (Eq. 2):

$$\begin{array}{}\text{(2)}& L\left({\mathit{\beta }}_j\right)=\prod _{j=\mathrm{1}}^n{P}_j^{{\mathit{\delta }}_j}(\mathrm{1}-P_j{)}^{\mathrm{1}-{\mathit{\delta }}_j},\end{array}$$

where n is the number of observations and δ is the Kronecker symbol:

$$\begin{array}{}\text{(3)}& \mathit{\delta }=\left\{\begin{array}{ll}\mathrm{1},& \mathrm{when}\mathrm{there}\mathrm{is}\mathrm{a}\mathrm{surface}\mathrm{fracture}\mathrm{visible}\\ & \mathrm{on}\mathrm{a}\mathrm{satellite}\mathrm{image}\\ \mathrm{0},& \mathrm{otherwise}.\end{array}\right.\end{array}$$

Once the values of β are found, we can find the probability of a node being fractured by substituting a chosen set of predictors into Eq. (1).

4.2 Predictor parameters

We started with a set of 19 predictors, x_i. Some of them are known to influence fracturing (stresses, strain rates, ice rheology), while others we considered to be potentially important (various geometrical properties, proximity to the ice front and the grounding line, etc). Temperature and accumulation were not included in the list of predictors due to the incompatibility of their spatial resolution with the relatively fine 100 m mesh we used to model fractures. They might be important for the formation and propagation of fractures, as warmer temperatures can increase the number of fractures due to the effect of meltwater (Mobasher et al., 2016; Weertman, 1973; van der Veen, 1998b), but a better-resolution climate data set would be needed to assess this.

All the experiments and sets of parameters used in LRA were constructed separately for floating and grounded ice. This is due to the fact that some parameters that were used for prediction of fractures on grounded ice are not applicable for predicting fractures on floating ice and vice versa (for example, friction and bed slope are irrelevant on floating ice, whereas back stress cannot be applied to grounded ice).

The calculation of some predictors was performed using methods already implemented in ISSM (e.g. stresses, strain rates, friction coefficient). Other predictors (e.g. calculation of curvature, distances to ice front, grounding line, proximity to glacier edges and nunataks) are not produced by ISSM and were calculated independently. Here we describe the methods we used to calculate each predictor parameter as well as a brief description as to why each parameter may have an impact on the location of fractures:

1. Principal values of the deviatoric stress and effective stress. Both stress and strain rates are calculated using the observed surface velocities. Following the Shallow-ice approximation, the deviatoric stress is

   $$\begin{array}{}\text{(4)}& {\mathit{\sigma }}^{\prime }=\left[\begin{array}{ccc}{{\mathit{\sigma }}^{\prime }}_{xx}& {{\mathit{\sigma }}^{\prime }}_{xy}& \mathrm{0}\\ {{\mathit{\sigma }}^{\prime }}_{xy}& {{\mathit{\sigma }}^{\prime }}_{yy}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -{{\mathit{\sigma }}^{\prime }}_{xx}-{{\mathit{\sigma }}^{\prime }}_{yy}\end{array}\right],\end{array}$$

   where ${{\mathit{\sigma }}^{\prime }}_{ij}$ are the deviatoric stress components.

   The deviatoric stress values have a direct effect on the opening and closing of crevasses; the sign of the first principal stress component determines whether it is compressive (negative) or tensile (positive).

   Effective deviatoric stress is calculated as

   $$\begin{array}{}\text{(5)}& {{\mathit{\sigma }}^{\prime }}_{\mathrm{e}}=\sqrt{{{\mathit{\sigma }}^{\prime }}_{xx}^{\mathrm{2}}+{{\mathit{\sigma }}^{\prime }}_{yy}^{\mathrm{2}}+{{\mathit{\sigma }}^{\prime }}_{xy}^{\mathrm{2}}+{{\mathit{\sigma }}^{\prime }}_{xx}{{\mathit{\sigma }}^{\prime }}_{yy}}.\end{array}$$

   Von Mises stress is calculated as

   $$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{\mathrm{vm}}=\sqrt{{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\sum _{i,j}{{\mathit{\sigma }}^{\prime }}_{ij}{{\mathit{\sigma }}^{\prime }}_{ij}}=\sqrt{\mathrm{3}}B{\dot{\mathit{\epsilon }}}_{\mathrm{e}}^{\mathrm{1}/n},\end{array}$$

   where B and n are the creep parameter and the creep exponent, respectively, and ${\dot{\mathit{\epsilon }}}_{\mathrm{e}}$ is the effective strain rate (Eq. 8).

2. Back stress. Back stress is calculated in ISSM using velocity inversion results as described in Borstad et al. (2013).

3. Effective strain rate. The effective strain rate ${\dot{\mathit{\epsilon }}}_{\mathrm{e}}$ is included in our analysis because it is known that crevasse initiation is linked to strain rates (Campbell et al., 2013). If the strain rate in the horizontal plane is sufficiently high, crevasses can propagate to greater depth (Benn and Evans, 2010). In addition, stresses can trigger brittle fracturing; however, to take into account a gradual viscoelastic effect that can lead to fracture formation, strain rates are included in our model.

   The principal strain rates are calculated using the observed velocities as eigenvalues of the matrix

   $$\begin{array}{ll}\text{(7)}& {\displaystyle }& {\displaystyle }\dot{\mathit{\epsilon }}={\displaystyle }& {\displaystyle }\left[\begin{array}{ccc}{\displaystyle \frac{\partial u}{\partial x}}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)& \mathrm{0}\\ {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)& {\displaystyle \frac{\partial v}{\partial y}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -{\displaystyle \frac{\partial u}{\partial x}}-{\displaystyle \frac{\partial v}{\partial y}}\end{array}\right],\end{array}$$

   where u and v are horizontal components of the observed surface velocity.

   Using again the shallow-ice approximation, vertical shear is neglected, and the effective strain rate is approximated as

   $$\begin{array}{}\text{(8)}& {\dot{\mathit{\epsilon }}}_{\mathrm{e}}=\sqrt{{\dot{\mathit{\epsilon }}}_{xx}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{yy}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{xy}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{xx}{\dot{\mathit{\epsilon }}}_{yy}},\end{array}$$

   where ${\dot{\mathit{\epsilon }}}_{ij}$ are the strain rate components (since in 2-D we neglect ${\dot{\mathit{\epsilon }}}_{xz}$ and ${\dot{\mathit{\epsilon }}}_{yz}$ and use incompressibility ${\dot{\mathit{\epsilon }}}_{zz}=-{\dot{\mathit{\epsilon }}}_{yy}-{\dot{\mathit{\epsilon }}}_{xx}$).

4. Horizontal strain rate change. Changing geometry or boundary conditions can cause changes in strain rate and can also cause fractures (e.g. a glacier flowing over a convex slope or icefall: the change in bed slope causes a change in strain rate and also causes fractures. However, it is important to mention that there might be cases when change in strain rate is not the cause but the consequence of the presence of fractures. However, the aim of our study is to identify where fractures are present without attempting to fully describe the process by which they are formed. Thus, we use the change in strain rate as a predictor precisely because it tells us where we can expect to find fractures. This predictor allows us to discover new regions where crevasses are present even if they were not seen in the imagery. A lack of observed fractures but high strain rate might mean that fractures may not be visible but should still be present (e.g. if fractures are covered in snow or not visible due to bad resolution of the satellite images) or might be due to an error in the velocity field. This parameter is calculated as a maximum difference between the strain rate at the location of observations and strain rate upstream.

5. Velocity. Velocities are included in the model because the increase of the ice flow velocities can result in more fracturing. The observations of the surface velocities are taken from InSAR.

6. Friction (removed from the analysis, described in Sect. 4.2.1). Low friction levels at the base of glaciers will lead to a higher sensitivity to membrane stresses, which can lead to more crevassing in tensile mode. We obtain this parameter from the inversion of surface velocities in ISSM.

7. Rheology B (stiffness of ice). In addition, we include the viscosity parameter B in Glen's flow law because, when ice stiffness increases and ice crystals cannot creep fast enough, fracture may occur. Therefore, this parameter is added as a predictor. Adding temperature directly into the analysis did not improve the prediction results, which might be due to the uncertainties in the temperature estimation. Ice stiffness is obtained from the inversion of observed velocities (implemented in ISSM).

8. Ice thickness. Ice thickness was included due to the fact that fracture formation in thicker glaciers/ice shelves might differ from fracturing in thin glaciers/ice shelves. This parameter was taken directly from the Bedmap2 model.

9. Distance to the ice front and the grounding line. We can see in the satellite images that more fractures are present at a certain distance from the ice front as well as near the grounding line. We found that the relation between the presence of fractures and distance to the ice front as well as the distance to the grounding line is non-linear (Fig. 4). For most ice shelves/glaciers we can see more fractures 3–5 km as well as 10–13 km away from the front and a slightly smaller number of fractures closer than 3 km to the front or between 5 and 10 km. Therefore, instead of using d_IF and d_GL (distance to the ice front and the grounding line in kilometres) as predictor variables, we construct dummy variables: DM_IF and DM_GL, respectively, which represent two-column arrays in the following form:

   $$\begin{array}{}\text{(9)}& {\mathrm{DM}}_{\mathrm{IF}}=\left\{\begin{array}{ll}(\mathrm{1},\mathrm{1}),& \mathrm{when}\mathrm{3}\mathrm{km}\le d_{\mathrm{IF}}<\mathrm{5}\mathrm{km}\\ & \mathrm{or}\mathrm{10}\mathrm{km}\le d_{\mathrm{IF}}<\mathrm{13}\mathrm{km}\\ (\mathrm{1},\mathrm{0}),& \mathrm{when}\mathrm{5}\mathrm{km}\le d_{\mathrm{IF}}<\mathrm{10}\mathrm{km}\\ (\mathrm{0},\mathrm{1}),& \mathrm{when}{d}_{\mathrm{IF}}<\mathrm{3}\mathrm{km}\\ (\mathrm{0},\mathrm{0}),& \mathrm{else}\end{array}\right.,\end{array}$$
   $$\begin{array}{}\text{(10)}& {\mathrm{DM}}_{\mathrm{GL}}=\left\{\begin{array}{ll}(\mathrm{1},\mathrm{1}),& \mathrm{when}\mathrm{5}\mathrm{km}\le d_{\mathrm{GL}}<\mathrm{15}\mathrm{km}\\ (\mathrm{1},\mathrm{0}),& \mathrm{when}{d}_{\mathrm{GL}}<\mathrm{5}\mathrm{km}\\ (\mathrm{0},\mathrm{1}),& \mathrm{when}\mathrm{15}\mathrm{km}\le d_{\mathrm{GL}}<\mathrm{20}\mathrm{km}\\ (\mathrm{0},\mathrm{0}),& \mathrm{else}\end{array}\right..\end{array}$$

10. Proximity to glacier edges. Generally the lateral friction along the glacier boundary is not considered in ice sheet models when stress is calculated. The stress field alone can predict transverse, longitudinal and radial splaying crevasses. They are all formed due to opening stress and are normally considered in existing damage modelling methods. However, the prediction of crevasses near the edges of glaciers requires a parameterization of the lateral drag. Thus, we include the proximity to edges of glaciers and to nunataks as a predictor in our model.

11. Surface change. We include this predictor in the analysis due to the fact that fracturing can be caused by an increase in stress due to an abrupt change in surface elevation. The change of surface elevation is calculated using the topography data from Bedmap2.

12. Curvature. It is clearly visible in satellite images that more fractures occur around horizontal bends in glaciers. Therefore, the curvature of the glacier channel was included as a predictor, calculated as

    $$\begin{array}{}\text{(11)}& \mathit{\alpha }=\arccos \left({\displaystyle \frac{\mathit{v}\left(P\right)\cdot \mathit{v}\left(E\right)}{\left|\mathit{v}\right(P\left)\right|\cdot \left|\mathit{v}\right(E\left)\right|}}\right),\end{array}$$

    where v(P) is the ice horizontal velocity at the point of observations and v(E) is the velocity D metres away from the point. The distance D is based on the velocity magnitude v(P), because if the velocity is high, we need to increase D so that two subsequent points capture the geometry of the bend of a glacier. Thus, if v(P) is greater than 2000 m yr⁻¹, we assign D=3v(P); otherwise we assign D=6v(P). These values are not arbitrary: this assignment is used in the model only to have enough data points to see the local curvature of a glacier, and it does not affect the calculation of the curvature itself.

13. Bed and surface slopes. Calculation of the stress field already considers surface velocity, surface slope and a curvature of a glacier channel (Larour et al., 2012), but for our method we estimate the effect of each parameter on fracturing separately, because including only stress was not sufficient to explain fracture formation (described in Sect. 4.2.1). Thus, bed and surface slopes are included in the model since shear stress increases on a steeper slope and can lead to fracturing (e.g. ice fall is an extreme case). The calculation of these parameters is performed by computing the derivative of the bed and surface using the finite-element nodal function in ISSM (a linear function calculated for each node using the information about the horizontal coordinates).

Figure 1 Comparison between the success of identifying fractures using LRA (purple for grounded ice, green for floating ice) and using test set 1: effective deviatoric stress (effective dev. stress); test set 2: principal deviatoric stress 1 and 2; and test set 3: von Mises stress (blue for grounded ice, yellow for floating ice). Left column represents grounded ice (a, c, e), and right column show the results for floating ice (b, d, f).

Finally, due to the fact that all predictor parameters have different units, as well as significantly different magnitudes, we normalize each predictor used in Eq. (1) as follows:

$$\begin{array}{}\text{(12)}& x_i^{*}={\displaystyle \frac{x_i-{\mathit{\mu }}_i}{{\mathit{\sigma }}_i}},\end{array}$$

where μ_i and σ_i are the mean and standard deviation of the predictor variables, respectively.

4.2.1 Test run with a small set of parameters

Including stress variables to predict fractures is intuitive as they are one of the major indicator of ice being fractured or non-fractured. Other variables such as geometry correlate to stress variables, but we found that it is important to include them in the model because the results are inferior if the parameters are not included. This might be caused by limitations of the predictor parameter values produced by the ice sheet model or the simplification of the Stokes equation.

In order to show that our fracture model works better when including both physics-based and geometry-based predictors, we ran three additional experiments. In the first test run we included only the effective deviatoric stress as a predictor and found that, although it produces reasonable results matching the observations, the success of identifying fractures is about 20 % lower than the results of the model with the chosen optimal set of the predictors (Fig. 1a, b). The second set of experiments contained only the principal deviatoric stresses and produced results that did not agree well with the observations (Fig. 1 c, d). Similarly, the results of the third test, which included only von Mises stress as the predictor, did not agree with the observed fractures (success rate not exceeding 50 %; Fig. 1 e, f). This shows that stress derived from velocity observations through a model relying on SSA is not sufficient to model fractured zones and that a more complex set of predictor parameters is required.

Moreover, including both friction and strain rate might be ambiguous since less friction can lead to larger strain rates. By looking at the predictor data sets, we found that the optimal choice of parameters for each group includes either friction or strain rate but never both at the same time. We ran an experiment replacing strain rates by friction and found that the prediction success for some glacier decreases by only about 3 %. We therefore kept only strain rate as a predictor parameter and discarded friction.

4.3 Optimization the choice of predictors

To construct the probability function for each glacier, we sought a set (or subset) of the predictor factors required to include in the LRA. We started with a first guess (calculated using LRA and a potential choice of the predictor parameters) and then improved it based on three methods: random walk, Bayesian and Jensen–Shannon divergence algorithm.

4.3.1 Random walk

For each of the 45 ice shelf regions we performed a 100 000-step run with random sets of predictors used at each step (the number and selection of predictors were chosen at random at each step). We defined a potentially good model to be the one with a success of identifying fractures larger than 70 % and the error rate of overestimation not exceeding 15 % (however, when a good-fit model was not found after 2000 steps, we looked for a model with a 65 % success rate and 20 % error rate). Once a good fit was found, we saved it as a potentially good set and continued running the model with different sets of predictor factors for the remaining number of steps to search for a better model. At the end of each run the algorithm provided us with a mean set of factors for a best-fitting model.

4.3.2 Bayesian

To test the behaviour of the models with different sets of parameters and, thus, to choose with more precision an optimal set of predictors from the full set, we performed a non-linear Bayesian inversion.

To find the likelihood function for a Bayesian inversion, we need to add the probabilities of fracturing for all nodes. The area of each ice shelf region is ∼10⁸ km², which, when adding for all nodes, leads to a very large sum of all modelled probabilities (using LRA) for non-fractured nodes and therefore extremely large likelihoods (note that these values should not be confused with 0 and 1 values set for observed fractures only). In order to achieve a more realistic magnitude of the likelihood function and to account for the fact that in general there are more observed non-fractured nodes, we re-calculated the estimated LRA probabilities by scaling them between 0.55 and 1. To do this, we assigned all probabilities below 0.55 to the value of non-fractured nodes (value of zero) and scaled the remaining values to the range from 0 to 1.

In addition to defining a likelihood function, Bayesian inversion requires an input of prior model and prior scores. For a prior model, we took a calculated fracture probability p^∗ from the best-fitting model obtained earlier using the random-walk search. We assumed that the prior probability density function (PDF) was uniformly distributed between 1 and 17 (U[1,17], because the maximum number of predictor factors in a set is equal to 17).

Finding an expression for a likelihood function, L_i, for our model was problematic. We tested a number of different commonly used expressions, such as

$$\begin{array}{}\text{(13)}& \left\{\begin{array}{ll}{L}_i& =\sum _i\log \left(f_i\right),\mathrm{where}\\ f_i& =(\mathrm{1}-p_i{)}^{\mathrm{1}-d_i}+p_i^{d_i},\end{array}\right.\end{array}$$

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{ll}{L}_i& =\sum _i\log \left(f_i\right),\mathrm{where}\\ f_i& =(\mathrm{1}-p_i)\cdot (\mathrm{1}-d_i)+p_i\cdot d_i,\end{array}\right.\end{array}$$

where d_i and p_i represent observed and modelled fractures on an ice shelf/glacier, respectively, and f_i is the probability of agreement between two predictions in a cell i.

However, all of them produced very large likelihoods which increased dramatically with a small percentage change in the probability density function. The value of the likelihood function increased up to an order of 10⁵, which was unrealistic and made the inversions unstable. Therefore, it was crucial to choose a better representative likelihood function.

In order to construct the function, first we assumed that the measure R of the total agreement between two models (the sum of all probabilities) followed a Gaussian distribution with a mean E (Eq. 15) and a standard deviation σ (Eq. 16):

$$\begin{array}{}\text{(15)}& E\left(f_i^{\mathrm{pred}}\right)=\sum _{i=\mathrm{1}}^N{f}_i,\end{array}$$

$$\begin{array}{}\text{(16)}& \mathit{\sigma }=\left|E\left(f_i^{\mathrm{obs}}\right)-E\left(f_i^{\mathrm{best}}\right)\right|,\end{array}$$

where the two expected values for both data and the chosen model, $E\left(f_i^{\mathrm{obs}}\right)$ and $E\left(f_i^{\mathrm{best}}\right)$, respectively, are defined as

$$\begin{array}{}\text{(17)}& E\left(f_i^{\mathrm{obs}}\right)=\sum _{i=\mathrm{1}}^N{d}_i^{\mathrm{2}}+(\mathrm{1}-d_i{)}^{\mathrm{2}},\end{array}$$

$$\begin{array}{}\text{(18)}& E\left(f_i^{\mathrm{best}}\right)=\sum _{i=\mathrm{1}}^N{p}^{\ast }\cdot d_i+\left(\mathrm{1}-p^{\ast }\right)\cdot (\mathrm{1}-d_i),\end{array}$$

where p^∗ is the best-fit probability.

Second, our idea was to calculate the likelihood L_i as an exponential function of the misfit ϕ_d(m) between the data and the model, assuming that either the data (observed fractures) or the analysed model contains an error (Eq. 19):

$$\begin{array}{}\text{(19)}& L_i=e^{-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\varphi }}_d\left(m\right)},\mathrm{where}\end{array}$$

$$\begin{array}{}\text{(20)}& {\mathit{\varphi }}_d\left(m\right)={\displaystyle \frac{{\left(E\left(f_i^{\mathrm{obs}}\right)-E\left(f_d^{\mathrm{pred}}\right)\right)}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}}}.\end{array}$$

We performed a Bayesian analysis for 500 steps, then narrowed down the selection and accepted only those models that had likelihoods greater than 90 % of the best likelihood. Each step included two criteria: if a new likelihood was greater than the prior likelihood or was greater than a certain percentage (taken at random at each step) of the old likelihood, we accepted the model. This allowed us to identify the most commonly chosen sets of parameters.

4.4 Glaciers classification and Jensen–Shannon divergence (JSD)

In order to select a set of predictors for a general case and to find whether it is possible to identify a set that can be used for any ice shelf region, we started with the construction of a binary array for each ice shelf/glacier, where the number of rows represents the number of well-fitting models for an ice shelf/glacier and the number of columns represents each of the predictor factors.

We then found the average occurrence of each predictor:

$$\begin{array}{}\text{(21)}& A_i={\displaystyle \frac{\mathrm{1}}{N}}\sum _{j=\mathrm{1}}^N{k}_j,\end{array}$$

where $i\in [\mathrm{1},\mathrm{17}]$ is the predictor index and N is the number of well-fitting models. k_j=1 when the predictor is included in the well-fitting model j and 0 otherwise.

We could then determine how often a certain predictor was included in the good-fit models. If a predictor was selected more than 50 % of the time, then it was assigned as a candidate for the best-fitting model. Thus, we obtained a 45×17 array (45 glaciers vs. 17 predictors) that consisted of 1 when the predictor was included in the best-fit model and 0 otherwise.

Then, we classified the glaciers in groups. There were a large number of possible combinations to select such groups. Therefore, we constructed a test that assessed every possible combination and calculated a percentage of similarity between glaciers in a group (Eq. 22).

$$\begin{array}{}\text{(22)}& S={\displaystyle \frac{M}{\mathrm{17}}}\times \mathrm{100},\end{array}$$

where M is the number of matches between sets of predictors for two glaciers and S is a group number.

Finally, we found that we could categorize all 45 glaciers into four different groups, with group 1 having glaciers/ice shelves that can be more easily combined and group 4 being a narrower group of specific glaciers/ice shelves that cannot be placed in any of the other three groups.

The JSD method (Dagan et al., 1997) can be used as a tool to measure the distance between two distributions and can provide a value that we can use to assign a particular glacier/ice shelf to one of the groups. The JSD formula is widely used in statistics to measure a divergence of one probability distribution from another. We applied JSD to identify the similarity between the best probability for each glacier and a probability calculated by placing the glacier in a certain group.

The Kullback–Leibler divergence (Kullback and Leibler, 1951) is defined as

$$\begin{array}{}\text{(23)}& \mathrm{JSD}\left(P_{\mathrm{1}}\right|\left|P_{\mathrm{2}}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}D\left(P_{\mathrm{1}}\right|\left|M\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}D\left(P_{\mathrm{2}}\right|\left|M\right),\end{array}$$

where two conditional PDFs, $D\left(P_{\mathrm{1}}\right|\left|M\right)$ and $D\left(P_{\mathrm{2}}\right|\left|M\right)$, are defined as

$$\begin{array}{}\text{(24)}& D\left(P_{\mathrm{1}}\right|\left|M\right)=\sum \left(P_{\mathrm{1}}\log \left({\displaystyle \frac{P_{\mathrm{1}}}{M}}\right)\right)\end{array}$$

and

$$\begin{array}{}\text{(25)}& D\left(P_{\mathrm{2}}\right|\left|M\right)=\sum \left(P_{\mathrm{2}}\log \left({\displaystyle \frac{P_{\mathrm{2}}}{M}}\right)\right),\end{array}$$

where $M=\frac{P_{\mathrm{1}}+P_{\mathrm{2}}}{\mathrm{2}}$ and P₁ and P₂ are the new probability (when assigning a glacier/ice shelf to a new group) and the old probability (the best-fit model), respectively. Both probabilities P₁ and P₂ have to be normalized before applying Eq. (23).

4.5 Calculation of damage

Here we utilize the damage-based model as an independent method in order to compare it with our statistics-based method. We do not compare our probability-based model with the damage model directly; rather, we evaluate their respective ability to predict the formation of fractures in ice. For this we compare calculated damage with the observations of fractures and identify areas where it can and cannot accurately predict the presence of fractures.

In this study we use the damage inversion method proposed by Borstad et al. (2013) to identify regions where fractures are initiated. Damage in this context has no vertical coordinate but comes from a linear mapping of the depth-integrated shallow-shelf equations. This approach is performed in two steps. First, the inversion of rheology based on the misfit between observed and modelled velocity is performed on floating ice. Then, damage is calculated as

$$\begin{array}{}\text{(26)}& D=\mathrm{1}-{\displaystyle \frac{B_{\mathrm{I}}}{B_{\mathrm{T}}}},\end{array}$$

where B_I and B_T are viscosity parameters calculated from an inversion and initialization of viscosity based on temperature analysis, respectively.

It is important to keep in mind that the inversion only infers damage in areas where fractures (crevasses or rifts) are being actively formed and, thus, creating a jump in strain rate/velocity. Many rifts are formed at one point in time and then only intermittently propagate. If the velocity observations do not show a discrete jump across a fracture, then there is nothing for the inversion to pick up in terms of damage. It only finds fractures that are actively enhancing the flow, and it is not meant to locate every fracture.

Estimation of B_T is the source of the main uncertainty in damage calculations due to the lack of ice temperature data, which can be crucial in affecting the accuracy of the viscosity parameter (Bassis and Ma, 2015). Thus, the errors in assumed temperature may affect the inferred value for damage.

Fractures that have been advected can be identified by damage, but this is not always the case, due to the fact that the inverse method for calculating damage will only find damage where there are fractures that give rise to velocity gradients. Damage will capture some fractures that were formed upstream and advected to a region with different stress conditions only if the fracture enhances the flow and creates a local velocity gradient. Thus, we first calculate flow lines for each observed fracture. If upstream from the fracture the damage is larger than 50 %, we assume that the damage calculation may be correct and that the observed fracture was formed upstream. If there is no damage initiated at the point or damage upstream from the observed fracture, we assign the observation point as not captured by the damage method and consider this as a failure of damage to identify the fracture (which can be due to the fact that the fracture in observation point does not cause a local gradient in strain rate).

Physics-based methods, such as LEFM and CDM, are necessary when modelling fractures in Antarctica. We do not intend to substitute these methods; rather, we seek a method that can improve on some aspects and cases when physics-based models do not predict well the formation of fractures. In particular it is possible that some fractures are initiated upstream from the grounding line rather than on floating ice. It is therefore important to be able to predict the formation of fractures in both cases. Damage is calculated only on floating ice based on model inversions using ISSM (Larour et al., 2012) because it is not possible to distinguish between basal friction and damage on grounded ice, as they have similar effects on the ice velocity. Thus, the main motivation of this study is not to replace the damage approach, which in fact provides a strong physical background for ice sheet modelling, but to find an alternative method that can be applied to both ice shelves and grounded ice, can work for a large set of glaciers/ice shelves and does not depend on temperature observations or threshold parameters.

5 Results

We applied the LRA method combined with the random-walk method to 45 ice shelf regions, including both ice shelves and surrounding grounded ice (the corresponding names and locations can be found in Table 6 and Fig. 3, respectively), and found a best-fitting model for 44 of them. The fracturing of the remaining ice shelf cannot be described using the predictors we have, which produce unacceptably large or small probabilities.

In total, for each ice shelf/glacier the random-walk analysis gave a number of possible sets of predictors that can produce a well-fitting model. We combined all of these possible sets for each glacier to see which predictors are always present in the well-fitting model and which ones are never included. The results of the random walk and the Bayesian inversion agreed well. Most of the essential predictors for each particular glacier selected in the Bayesian approach were also chosen when performing the random walk. In most cases, the Bayesian analysis showed equal importance of most of the predictors although effective strain rate and velocity had a slightly higher rate of selection. There was no universal set of factors that could be used to model all ice shelf regions. However, subsets of glaciers had some similarities in terms of the predictors that had to be included in order to achieve a well-fitting model.

To estimate how well our probability model and the damage model identify observed fractures, we calculated the percentage of success and error for each ice shelf/glacier model. First, we found the number of cases when there is a modelled fracture in the vicinity of an observed fracture (within 100 m radius). Then, we divided this number by the total number of observed fractures to find the percentage of success. To find the percentage of failure, we calculated how many times there is a modelled fracture when there are no observed fractures within a 100 m radius. We divide this number by the total number of non-fractured nodes to find the failure percentage.

Table 3 Formed groups of ice shelf regions.

Table 4 The value of the calculated weight coefficient β in LRA for each of the included parameters on grounded ice (rounded to two decimal places).

It is important to mention that β does not define wether a parameter is included in the model or not. The β coefficients are calculated only for the parameters that were already included in the model as they are a part of the model itself; it is the weight defining the calculation of the probabilities of fracturing (see Sect. 4.1). No value indicates that the predictor is not included in the model.

Table 5 Same as Table 4 but for floating ice.

Figure 2 Success and error percentages for LRA for grounded ice are shown in (a). Results for the floating ice applying LRA vs. damage method is shown in (b). For glacier 28 on floating ice (b) the success rate for both amage and LRA is 100 %; therefore, to make both bars visible, LRA is shown to have a slightly lower value than damage for this glacier.

Thus, we categorized the 45 glaciers/ice shelves into four groups, requiring that the deviation from the best-fit models did not exceed 5 %. Next, we performed a test to assess whether these selected sets were the optimal choice, by estimating the deviation from the best solution using the Jensen–Shannon divergence algorithm. We assigned each glacier to a particular group based on its minimum value of the deviation from the best-fitting model in JSD analysis. In so doing, we slightly modified the members of each group that we had previously created. For example, glacier 27 belonged to group 1 previously, and it fit well with only a slight change of the best-fit score. However the JSD showed that, if we move this glacier to group 2, the deviation from the best-fit decreases from 0.01 to 0.003. However, we had to take into account the fact that the JSD algorithm measures the total distance to the best-fit probability and, thus, can decrease the overestimation error while, at the same time, significantly decreasing the success rate. This took place for six glaciers/ice shelves with the numbers 10, 13, 15, 11, 30 and 32 (see Table 3). Therefore, since these six glaciers were of a similar type to the glaciers in group 1 and their JSD was similar for group1 and group 2 (e.g. JSD = 0.02 in group 2, and JSD = 0.0205 in group 1), they were assigned to group 1 to avoid a decrease in the success rate of identifying fractures.

Finally, to reach an optimal agreement between our model and the observations of fractures, we assigned each glacier to a particular group, and the set of factors for each group are presented in Tables 3, 4 and 5 (18 predictors, as friction was discarded). We found that the sets of predictors for each group varied significantly; however surface velocity was included in the grounded ice set for all groups. For the floating ice the analysis showed that surface velocity was not a determining factor in predicting fractures; instead effective strain rate and deviatoric principal stress values were present in each predictor set.

While the success rate of identifying fractures on floating ice was lower than for grounded ice, we were still able to identify the main fracture patterns, and the success rate was high for the majority of ice shelves (see Fig. 2a). Our method is able to identify up to 99 % of the exact location of fractures on grounded ice with an average of 84 % (Fig. 2a) and 61 % for floating ice (Fig. 2b). The mean overestimation error rate for grounded ice and floating ice was 26 % and 20 %, respectively. There are many cases where our method agrees with the results produced by the damage-based approach. However, in almost all cases the success rate of LRA on floating ice was higher than that of the damage-based method, with the exception of two glaciers. Overall, for all four groups where we could not achieve a high score using LRA, the damage-based method did not produce a high success score either (see Fig. 2b).

5.1 Group 1

This was the largest group of glaciers, and the best-fit model includes as many as 10 predictors for grounded ice and seven predictors for floating ice. The analysis of the estimated coefficients in LRA showed that predictors with the highest weights in our model for this group of glaciers were effective strain rate, proximity to glacier edges and nunataks, and the surface elevation change. We present the modelled probability of fractures in Figs. 6b and 5c as well as comparison with the damage-based results in Figs. 8a and 11c.

Figure 3 The location of each of the ice shelf regions in our analysis.

Figure 4 The number of observed fractures versus distance from the ice front.

Table 6 A list of analysed ice shelf regions. Fracture observations for the model calibration process were taken from regions marked with *.

AP: Antarctic Peninsula; WA: West Antarctica; EA: East Antarctica.

The main pattern of surface fractures is well represented for this group. On grounded ice the success rate of identifying fractures is larger than 88 %, with a quarter of glaciers at almost 100 %. The failure related to overestimation of fractures is 27 %. On floating ice the success amounted to 55 %, and the failure was equal to 15 % on average. For Vanderford IS (see Fig. 6b) the overall pattern is well represented (success rate of 91 % and 70 % for grounded and floating ice, respectively), even though high-resolution images were not available for this glacier. The overestimation error is mainly related to the region that is far from the ice front and has a relatively high accumulation rate, possibly obscuring the fractures in the imagery. On floating ice the probability of fracturing is relatively smaller, mainly showing a higher chance of fracturing closer to the groundling line. Conversely, Drygalski Ice Shelf has a larger number of high-resolution areas and, as a possible result, less overestimation of fracturing (see Fig. 5c). We can see that the “definitely non-fractured nodes” (selected in blue-ice areas) are successfully represented in our model. For this glacier, none of the observed non-fractured nodes was assigned a high probability of fracturing, with the modelled probability being as low as 0.1. Moreover, in the regions with a large number of observed fractures, the probability is as high as 0.9, and it is slightly lower in the areas with a smaller number of observed fractures (between 0.6 and 0.8). Observed fractures not captured by our model were not captured by the damage-based model either.

The modelling results for the Cook Ice Shelf are shown in Fig. 8b. There are distinct fractures visible towards the front and in the central part of the ice shelf that are not captured by either approach, which we interpret as showing that most of the fractures are formed further upstream near the groundling line. In general, the probability and the damage-based models show good agreement on the floating ice near the grounding line. However, damage does not reach 50 % at the majority of the locations. Moreover, at many locations where rifts are visible it shows 0 damage.

Figure 5 Modelled probability of a fracture for group 2, Rayner Thyner IS (a) and Edward VII IS (b), and group 1: Drygalski IS (c). Observed surface fractures are shown in green, and observed non-fractured ice is marked with blue crosses. Cyan boxes represent regions where high-resolution images were available. Black solid line shows the location of the grounding line.

Figure 6 Modelled probability of a fracture for group 2, Larsen A IS (a), and group1: Vanderford IS (b). Labels are the same as in Fig. 5.

The modelling results for Larsen B IS are illustrated in Fig. 11c. It is clear that the nodes where damage is high have a high probability of fracturing due to the fact that we added damage as one of the predictor parameters to this glacier. It can be also seen that there are two lines of high probability of fracturing that coincide with the location of the large rifts that can be seen in the satellite images.

The results for Nansen IS (Fig. 9c) as well as for Pine Island (Fig. 7a) agree well with observations even though the data from these two glaciers were not included in the calibration data set used to construct the LRA model. For Pine Island, we observe fractures in the central part of the shelf that were not captured by the model, but our model predicted high probabilities of fracture upstream where the ice is grounded. Thus, these fractures are likely advected from the grounding line out onto the floating ice shelf.

5.2 Group 2

The model for the second group of glaciers has the best fit when the bed slope is excluded. Effective strain rate and surface slopes were found to be the most important predictors in the model for this group.

For this group the LRA method predicts fractures with a 70 %–90 % success rate on grounded ice and finds about 67 % of observed fractures on floating ice with an overestimation of 25 % and 27 %, respectively. In most cases the model represents the non-fractured nodes with high precision, except for the slight overestimation at the front of the ice shelf. Similar situations are observed for most glaciers in this group: the area of floating ice is relatively small; thus the main prediction is performed for grounded ice. For, example, for Edward VII IS and Rayner Thyner IS (Fig. 5a and b, respectively) the modelled probability captures most of the fractured as well as non-fractured nodes on grounded ice with the exception of a few very small regions that are outside of the high-resolution image areas.

Interesting results were found for Larsen A IS (see Fig. 6a), showing a very good agreement between our model and the observations. We can see that the nodes observed to be non-fractured and located within the high-resolution regions are predicted by the LRA method to have lower probabilities of fracturing. Both the LRA and damage-based results for Shirase IS (Fig. 12a) are similar on floating ice, although the LRA method captures a slightly higher number of fractures. We infer that most of these fractures are formed further upstream on grounded ice.

Figure 7 Modelled probability for Pine Island (group 1) (a), modelled probability for Amery IS (group 4) (b) and inferred damage for Amery (c). Labels are the same as in Fig. 5. Red dots represent observed fractures that were considered to have advected downstream from areas identified as fracture formation zones by the damage model and, thus, were not used for evaluation of the damage-based results.

5.3 Group 3

This group includes four ice shelf regions, namely Totten IS, Nivl IS, Dibble IS and Holmes IS. These glaciers were very sensitive to the choice of predictor factors, and the JSD process could not assign them to either of the two aforementioned groups. The mean success rate for this group was around 93 %, with an overestimation rate just above 23 % on grounded ice and 56 % and 23 % success and error rate on floating ice, respectively. Potentially, in the model for this group we could include the proximity to the ice front since it produces slightly better results for three of the four glaciers. However, it lowers the success rate for Nivl IS significantly. Thus, in order to achieve a set that would give a well-fitting model for all of the glacier, we exclude back stress and the ice front proximity from the list of predictor factors for this group.

Figure 8 Modelled probability (a) and damage (b) at Mertz IS (group 1). Labels are the same as in Fig. 7.

In terms of results for Holmes IS (Fig. 9a), there was good agreement between damage and LRA models as they both were able to predict the main fracture pattern. However, the LRA predicted the observed fracture pattern slightly better, especially at the front of the ice shelf. Similarly, for Dibble IS (Fig. 10c), both methods produced similar patterns, although the LRA method had a higher overestimation error rate. Nevertheless, the LRA method was able to more precisely estimate fractures at the western part of the ice shelf.

Finally, we present the result for the Totten IS (see Fig. 12c). The images we used for this glacier were very hard to interpret due to the presence of many features on the ice surface as well as the low resolution of the imagery. Interestingly, for the Totten IS, only a certain set of factors can produce a good fit to observed fractures. For this glacier, including back stresses in the model produces a slightly smaller number of modelled fractures. Both the LRA model and the damage model capture most of the fractures we could observe on the floating ice, displaying similar distribution patterns. Although the overestimation is relatively high on grounded ice, it is not possible to say whether there are no fractures there and our model shows overestimation for this glacier or whether fractures are not visible due to the low-resolution images and our model provides a correct prediction.

Figure 9 Modelled probability of a fracture vs. modelled damage for Holmes IS (group 3) (a, b) and Nansen IS (group 1) (c, d). Labels are the same as in Figs. 5 and 8.

5.4 Group 4

This group includes Larsen C, Amery, George IV and Borchgrevnik IS. The average success rate for this group amounted to 66 % and 56 % for floating and grounded ice, respectively, while the average error rate amounted to 15 % and 20 %. The most important predictor factors in this group for floating ice are effective strain rate, surface change and ice thickness, while for prediction of fractures on grounded ice curvature and surface velocity need to be included in the model. For all of the ice shelves/glaciers in this group we found that including the ice front and grounding line proximity distorts the model, increasing significantly the error due to overestimation of fractures. For Borchgrevnik IS it also led to a drop of the success rate of fracture prediction. In addition, Larsen C and Amery ice shelves can be grouped together but cannot be included in any of the groups mentioned above. For these glaciers only a small number of predictors needed to be included in the model. The Bayesian analysis also confirmed the sensitivity of the Amery fracture model to this set of predictor factors. The LRA model for the Amery IS (Fig. 7b) was able to capture most of the fracture pattern on the grounded ice and near the edges of the floating ice and demonstrated a similar pattern to the damage model. For George IV (Fig. 10a) we observe a very good agreement with observations both on floating and grounded ice. There is a high probability of fracturing at most locations where the high damage zones are predicted. Overall, for both ice shelves we can see that the majority of fractures are formed further upstream of the grounding line. Plots for other glaciers/ice shelves are shown in Figs. S1–S8 in the Supplement.

6 Glacier characteristics

6.1 Group 1

The ice shelves/glaciers in this group have a number of characteristics that distinguish them from other ice shelves/glaciers. Most of them are relatively wide with a large floating area. The floating part is not restricted by any channel walls, and the width of the shelf is similar to its length. All glaciers in this groups are relatively static, with less curvature or significant surface elevation changes. However, there were two exceptions: Abbot IS and Drygalski IS have slightly different characteristics. First, Abbot is a wide glacier that exhibits most of the properties of group 1. However it has a large number of glaciers that restrict its outflow towards the ocean and, therefore, has similarities with the glaciers from group 4. This observation is in good agreement with the JSD results that showed that Abbot IS could also be assigned to group 4 as the change in JSD distance in this case would be very small. Second, the JSD results showed that Drygalski IS could be as well placed in group 2 or group 3. This ice shelf has some characteristics similar to group 2 (large number of nunataks) and group 3 (a very long floating tongue). Therefore, we suggest that some glaciers have mixed features of group 1–group 2 (such as Vanderford) or group 1–group 3 (Ekström, Tracy-Tremenchus, Rennik); however they still have more characteristics of group 1 and produce better-fitting results when assigned to this group.

6.2 Group 2

This group includes a relatively smaller number of ice shelves/glaciers. All of the ice shelves/glaciers have a large number of nunataks and smaller ice thickness as well as many small narrow channels and fast ice streams. They are mostly located on the Antarctic Peninsula or near the Transantarctic Mountains. All of the ice streams are relatively steep, which may explain why it is necessary to include surface slopes in order to achieve a well-fitting model.

6.3 Group 3

Group 3 glaciers were found to have many similar features. Most of the ice shelf regions in this group have one relatively long glacier that flows inside an embayment. For most of them the ice shelf is much longer than it is wide, and they all have a very low glacier channel curvature. The surface velocities of these ice shelves/glaciers are relatively high, which explains why changes in strain rate and surface velocity are the most important predictors for this group.

Interestingly, although the average back stress for Totten IS is one of the highest out of all 45 ice shelf regions, including it in the model does not significantly change the fracture probabilities. Thus, apparently, even though predictors may have large magnitudes, they can make just minor contributions to the constructed probability, and other predictors dominate the fracturing process for the Totten IS. The effective strain rate is also one of the highest for Totten, but we found that it is not this predictor that most contributes to fracturing; rather it is the effective strain rate change. Thus, sudden changes in the flow regime of the glacier would be the most likely cause to promote an increase of the number of fractures.

6.4 Group 4

The JSD analysis has shown that Borchgrevnik IS could also be assigned to group 1, but it produced slightly better results being placed in group 4. On the other hand, Amery and Larsen C need to be strictly assigned to group 4 only. George IV IS and Amery IS have similar characteristics as they are both narrow and long (in fact much longer than any other ice shelves of this type in Antarctica) and are located inside an embayment. Although Larsen C IS is not inside an embayment, it is a significantly long and narrow ice shelf stretching around the coast. Borchgrevnik IS also has similar features to the Amery and George IV ice shelves as it is of a similar shape and is located inside a narrow channel. However, it does not have exactly the same characteristics as the other ice shelf regions in this group as it is much shorter, which could be why JSD showed that it could also be placed in group 1.

On the Amery IS (see Fig. 7b) most of the fracturing occurs upstream of the grounding line. There were a number of fractures (right-hand side of Fig. 7b, close to the ice front of the glacier) that could not be represented by our models. However, the uncertainty of the fracture observations in this area is high due to the difficulty distinguishing them from the surface features caused by basal crevasses.

Adding the proximity to the grounding line and the ice front as predictors did not produce a good fit for the Borchgrevnik IS because of the specific shape of this region. The distance between the ice front and the grounding line is very small relatively to other glaciers in our analysis.

6.5 Discussion

We found that, in general, the most important predictor factors for modelling surface fractures on grounded ice for all analysed glaciers were the surface velocity and the surface change (maximum difference between the surface elevation within 500 m radius), which is in agreement with the theory of possible mechanism of fracture formation (Colgan et al., 2016). Interestingly, the required parameters on floating ice were different from grounded ice, with effective strain rate and principal stress being the most important. Previous analysis was based on damage accounts for effective stresses, thickness and viscosity, but it did not include such predictors as proximity to glacier edges, nunataks and the grounding line as well as the curvature of a channel, which helped to improve the modelling of fractures on most ice shelves in our analysis.

We do not claim that all the predictors that were chosen in the final set for each group represent the exact fracture mechanisms for each glacier. For some ice shelf regions, sets containing different predictors can lead to results close to the best-fitting model. However, for some cases, such as Amery and Totten ice shelves, the number of well-fitting models is very limited. For example, by including the effective strain rate and proximity to the ice front in the analysis, we can achieve a better fit to the observations. Therefore, we conclude that some factors have a very strong effect on fracturing, while others are only minor for some glaciers. Ultimately, we seek only to be able to develop a model that can identify correctly the geographical location of fractures, not necessarily explain why they are there.

The main uncertainty of our method is related to the overestimation of the number of fractures. It could be argued that we predict fractures at the locations where no observations of fractures are detected due to the fact that they are not visible due to snow accumulation or coarse resolution. A possible solution to this could be to supplement satellite images with radar and seismic measurements (Delaney and Arcone, 2005; Navarro et al., 2005) or to acquire higher-resolution satellite images that can help to identify location of fractures even if they are hidden under the snow surface. However, for our method it would require a large set of observational data. In our paper we assume that 20 %–25 % rate of overestimation is reasonable, since most of our results show overestimation when predicting fractures in the areas around observed fractures. Our assumption is based on the fact that the ice regime conditions are similar within a 500 m radius, not implying any direct influence of the old fractures on the new fractures (Colgan et al., 2016). In addition, the area of high probability of a fracture is larger than the number of observed fractures mainly due to the fact that it is not possible to select all of the fractures in the satellite images manually. The observed fractures in our evaluation data set capture the main areas of fracturing assuming that surrounded nodes are likely to be fractured as well.

Figure 10 Modelled probability of a fracture vs. modelled damage for George IV IS (group 4) (a, b) and Dibble IS (group 3) (c, d). Labels are the same as in Figs. 5 and 8.

Figure 11 Modelled probability of a fracture vs. modelled damage for Moscow University IS (group 1) (a, b) and Larsen B IS (group 1) (c, d). Labels are the same as in Figs. 5 and 8.

Figure 12 Modelled probability of a fracture vs. modelled damage for Shirase IS (group 2) (a, b) and Totten IS (group 3) (c, d). Labels are the same as in Figs. 5 and 8.

A significant overestimation of predicted fracturing can be seen at the front of Vanderford IS (see Fig. 6b), suggesting that there is a very high chance of developing surface fractures in that area. It may be that fractures exist there that are just not visible in the satellite images. In fact, that region has a relatively high snow accumulation rate reaching 1 m yr⁻¹. After closer inspection of the satellite image areas where we see the overestimation error, we can recognize the presence of surface fractures, and we may have under-identified existing fractures in that area. Due to a very large number of fractures around Antarctica, the manual selection of all the fractured data points is a time-intensive process. Sampling with a higher spatial density would require an automated algorithm. However, the low-spatial-density sampling does not influence the results of the fracture probability (as mentioned previously, only diversity in sampling is important); it affects only the observations data set that we use to estimate the success of our model. Thus, under-identification of fractures can lead to an apparent overestimation of the error.

We looked at various properties of Ronne IS, for which we could not find a good approximation using any of the 17 predictors. We found that the Ronne IS has the lowest elevation change as well as the principal stress components. We do not have enough samples to cover values that are non-typical for the majority of glaciers, which may explain why we could not find a good-fit model for this ice shelf, neither with LRA nor using the Bayesian analysis. Thus, we conclude that our probabilistic model is not appropriate in this case.

The damage-based approach sometimes produces areas of high damage downstream of the observed fractures (Fig. 10d). If we assume that there are some fractures that are not visible in the images, it is still unclear why the modelled ice is not damaged upstream where observed fractures are present. If we could see damaged ice upstream from the observed fractures, we could assume that the ice was damaged and transported downstream where we can see fractures. However, in many results based on the damage approach we could not identify this type of behaviour. Even after correcting the observations of fractures by back-integrating the flow of ice, we still can find fractures that do not have zones of high damage upstream of the fractures. In addition, the method based on damage inversion predicts only damage on floating ice, whereas fractures are often formed upstream of the grounding line. In our probability-based model this type of behaviour is accounted for, and the model in most cases overestimates fractures only in the vicinity of or upstream of observed fractures. In Fig. 9a we show modelled probability of fracturing for Holmes IS. The overestimation rate is quite small, and it all occurs upstream of the observed fractures, which might be due to fractures that were formed further upstream being transported downstream where they are visible in the satellite images.

7 Conclusions

Most previous large-scale modelling of surface fractures has focused on applying zero-stress models (propagation of crevasses to the depth where the overburden pressure and the tensile stress are equal; Nye, 1955), linear elastic fracture mechanics and continuum damage mechanics. We have shown that, using the suggested nominal parameters, the damage-based approach does not fully reproduce the location of fracturing for any ice shelf region. In this study, we constructed a probability-based method to model surface fractures and generated improved predictions of fractures when physics-based models did not predict well the location of surface fractures. From this different perspective, we can construct an alternative method to predict the location of fractured zones not only on floating ice but also on grounded ice.

We found that the logistic regression algorithm, combined with other statistical methods, can significantly improve the prediction of fractured zones for the Antarctic ice shelves/glaciers and can lead to the identification of up to 99 % of observed surface crevasses for some ice shelf regions, with an average of 70 % for all ice shelf regions. Our approach has a number of uncertainties and leads to some overestimation of the number of fractures in comparison to the observations, but the rate is not significantly higher than the overestimation error found when using the damage-based method. However, the damage-based method did not predict location of many fractures either upstream or downstream from the observed locations, suggesting an underestimation when applying the damage method. The probabilistic results suggest that our statistics-based methods are more reliable in identifying fractures and rifts at the locations where the damage method does not predict them (which is related not to a failure in the damage method but to the fact that it is not constructed to do so). There are also uncertainties in the damage-based method related to the surface temperatures in Antarctica, which may be poorly represented with available observations. It is possible that the damage-based method needs to be tuned for every ice shelf separately, but this is beyond the scope of this study.

We classified the Antarctic ice shelf regions into four groups, where ice shelves/glaciers in each group have similar characteristics and each group has a set of predictors that can be used to predict the location of fractures. Although there were ice shelves/glaciers of specific shapes and having specific regimes that are more difficult to describe applying the general set of factors suggested in this study, overall our method provides a tool that can be used in the analysis of fracturing for most of the ice shelves/glaciers in Antarctica.

Our model is easy to implement and can be effectively used as a basis for modelling of fracture propagation and the first step in implementing a calving parameterization in ice sheet models. This statistics-based method can help to expand our current knowledge of the crevasses as well as improve mapping of potential hazards. Our results can be used to identify potential regions with snow-covered crevasses that may pose hazards for navigation in Antarctica and, thus, complement field campaigns.