Phase-field models of floe fracture in sea ice

We develop a phase-field model of brittle fracture to model fracture in sea ice floes. Phase fields allow for a variational formulation of fracture by using an energy functional that combines a linear elastic energy with a term modeling the energetic cost of fracture. We study the fracture strength of ice floes with stochastic thickness variations under boundary forcings or displacements. Our approach models refrozen cracks or other linear ice impurities with stochastic models for thickness profiles. We find that the orientation of thickness variations is an important factor for the strength of ice floes, and we study the distribution of critical stresses leading to fracture. Potential applications to discrete element method (DEM) simulations and field data from the ICEX 2018 campaign are discussed.

1 Introduction

The fracture of sea ice at intermediate scales, 100 m to 10 km, impacts ice transport at the basin scale on the order of hundreds to thousands of kilometers. However, current models do not account for intermediate-scale fractures whose behavior may respond to recent shifts in climate trends. A positive correlation between ice transport and deformation (Lewis and Hutchings, 2019; Rampal et al., 2009), coupled with the loss of thicker multi-year ice (Meier et al., 2014) and falling ice concentrations, suggests that mechanically weaker ice is fracturing faster, leading to earlier melt and increased advection. At these scales, kilometer-long fractures form across floes, contributing to lead formation and mechanically weaker sea ice with more open water. Advection is enhanced as the weaker ice pack yields to winds. Moreover, darker, open water lowers the effective albedo and increases melt. Incorporating such fracture processes either directly into discrete element methods (DEMs) or effectively into continuum models has been an important, long-standing challenge to improving sea ice predictions (Blockley et al., 2020; Dempsey, 2000; Weiss and Dansereau, 2017).

Modeling intermediate-scale fractures is challenging, and most smaller-scale observations are limited to the lab scale on the order of 10 m (Timco and Weeks, 2010; Dempsey et al., 2018; Schulson and Duval, 2009; Coon et al., 1998). In particular, available sea ice elastic moduli currently do not account for scale effects and effective parameters due to refrozen leads and ridges.

Previous studies have connected fracture angles and intermediate-scale models (Wilchinsky et al., 2011; Ringeisen et al., 2021; Plante and Tremblay, 2021; Dansereau et al., 2019; Hibler III and Schulson, 2000). One of the first detailed measurements of intermediate-scale deformation were made during the ICEX 2018 field campaign (Parno et al., 2022), yet, simultaneous measurements of forcings from the ocean and atmosphere and initial ice stresses are not available. In response, we aim to develop a model that can produce large ensembles of fracture events in realistic settings of floe geometry, material parameters, and forcing.

There are multiple candidates for fracture models in sea ice with various strengths and drawbacks. Lu et al. (2015) used an extended finite-element method (XFEM) that simulates fracture modeled by linear elastic fracture mechanics (LEFMs). The formulation of this model and the interpretation of lab-scale measurements of fracture parameters (Dempsey et al., 2018) are systematically based on LEFMs. The authors present models that reproduce realistic crack propagation. A drawback of this approach is that a crack tip must be inserted, but in general one does not know where a crack nucleates. Particle or discrete element methods (Tuhkuri and Polojärvi, 2018; Hopkins and Thorndike, 2006; Kulchitsky et al., 2017; Jirásek and Bazant, 1995) model fracture as a process that emerges from the failure of elastic bonds between particles or jointed edges between polygons. They are widely used in modeling ice–structure interactions. Computationally, fracture profiles from these methods are limited by particle or element geometry. In particular, when there is a need to explore fracture profiles over floe geometries, such experiments are computationally prohibitive. Montiel and Squire (2017) model floes as thin elastic disks that fracture at critical stresses. Their fracture profiles are linear and must be parameterized. While this model could be extended to other modes of fracture and floe shapes, the additional step of parameterizing the space of one-dimensional cracks on a two-dimensional floe is challenging. Finally, one could assume fracture is scale invariant, and apply models from larger scales. There is extensive literature (Rampal et al., 2019; Hutter et al., 2019; Bouchat et al., 2022) to support such an approach. However, those models aggregate intermediate-scale processes that emerge from the lab scale. In other words, one would need to boldly conjecture that self-similarity of fracture processes holds from thousands of kilometers down to tens of meters.

Phase-field models of fracture (Bourdin et al., 2008) are an appealing choice for modeling floe-scale fractures. They are based on an energetic formulation of deformation and fracture. They can directly incorporate lab-scale observations because they are based on Griffith's theory of brittle fracture (Griffith, 1921), the foundation of LEFMs. Phase fields regularize cracks and allow them to be described by the variations in a scalar phase field s. Phase-field models have been used to explore crack profiles in different geometries, materials, and physical settings. In particular, researchers have used phase fields to simulate processes such as fracture nucleation, propagation, branching, and fragmentation (Ambati et al., 2015) in a wide range of materials, including concrete, steel, and biological tissue. The literature on these models includes an extensive review of numerical implementations (Wu et al., 2020) and extensions to a variety of elasticity constitutive equations and types of failure. We aim to build on these results and on lab-scale ice measurements to model intermediate-scale sea ice fracture.

In this work, we investigate floe-scale fracture with phase-field models. We simulate fracture under boundary force or displacement conditions for a distribution of stochastic ice thickness fields, which model refrozen ice cracks or ice ridges. Our experiments show that the critical forcing at which fracture occurs depends on the geometry of the thickness anomalies and their relation to the forcing to which the floe is subjected. Additionally, we discuss measurement data from the ICEX 2018 expedition in fracture simulations and prospects for incorporating physical fracture in discrete element models.

The outline of the paper is as follows. In Sect. 2, we discuss the phase-field model and numerical implementation. In Sect. 3, we present our experiments on stochastic weaknesses. Finally, in Sect. 4, we consider future research directions.

2 Methods

In this section, we look at how phase-field models of brittle fracture can be used to model kilometer-scale fractures. We define a phase-field formulation of brittle fracture with linear elasticity in Sect. 2.1. In Sect. 2.2, we review a staggered algorithm to solve the phase-field model formulation when ice floes are subject to displacement or force boundary conditions. We review sea ice parameters from the literature and explain how we adapt them to our framework, and we discuss the implementation of the approach in the open-source finite-element method (FEM) library FEniCS (Alnaes et al., 2015) in Sect. 2.3 and 2.4.

2.1 Phase-field model of brittle fracture

We describe ice floe deformation using a linear constitutive relation. The floe geometry is modeled by a two-dimensional domain Ω∈ℝ², on which we consider displacement fields $\mathit{u}=(u_{\mathrm{1}},u_{\mathrm{2}})$. We assume that unbroken ice satisfies plane stress linear elasticity equations; i.e., the stress tensor induced by the strain tensor $\mathit{\epsilon }:=(\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^{\text{T}})/\mathrm{2}$ is

$$\begin{array}{}\text{(1)}& \mathit{\sigma }=\left[\mathit{\lambda }\text{tr}\right(\mathit{\epsilon })\mathbf{I}+\mathrm{2}\mathit{\mu }\mathit{\epsilon }],\end{array}$$

where λ and μ are the Lamé parameters, tr(⋅) denotes the trace operator, and $\mathbf{I}\in {\mathbb{R}}^{\mathrm{2}\times \mathrm{2}}$ is the identity matrix. The definition 1 results in a corresponding elastic energy density $\mathrm{\Psi }\left(\mathit{u}\right):\mathrm{\Omega }\to \mathbb{R}$ defined as

$$\begin{array}{}\text{(2)}& \mathrm{\Psi }\left(\mathit{u}\right):={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\sigma }:\mathit{\epsilon }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathit{\lambda }\left(\text{tr}\right(\mathit{\epsilon }){)}^{\mathrm{2}}+\mathit{\mu }\text{tr}({\mathit{\epsilon }}^{\mathrm{2}})\right].\end{array}$$

Constitutive relations other than Eq. (1) can be used, but in this study we focus on a linear rheology and the corresponding well-studied fracture process.

We use a phase-field model of quasi-static brittle fracture developed by Bourdin et al. (2000). Their work builds on Francfort and Marigo (1998), who developed an energetic formulation of Griffith's theory of brittle fracture. Francfort and Marigo (1998) extended Griffith's formulation by adding a surface or fracture energy that describes the cost to create a one-dimensional crack set Γ. In this approach, displacements are allowed to be discontinuous across Γ. Moreover, the surface energy is the arc length of Γ times the fracture toughness G, the energetic cost per unit length to fracture material.

Determining a crack set Γ that minimizes the total energy is a numerically challenging problem (known as the free-discontinuity problem; Farrell and Maurini, 2017), as the occurrence of fracture is part of the solution and thus a priori unknown, making it difficult to resolve using a numerical discretization. As a remedy, Bourdin et al. (2000) propose an approach that replaces the crack set with a continuous, scalar phase field $s:\mathrm{\Omega }\to \mathbb{R}$, describing a diffusive crack profile. The phase field takes values in the interval [0,1]. Fracture occurs where s=0, and ice is unbroken where s=1. Between the two values, s transitions continuously. A length-scale parameter ℓ controls the spread and smoothness of s near a crack. Bourdin et al. (2000) show that as ℓ→0 in the case of linear elasticity, the 0 set of s converges to a minimizing crack Γ.

Combined with other constitutive relations, phase-field models have also been used to capture different modes of fracture (Ambati et al., 2015). The total energy ℰ underlying these phase-field models is as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathcal{E}(\mathit{u},s)& :=\underset{\text{elastic energy}}{\underbrace{\underset{\mathrm{\Omega }}{\int }(s^{\mathrm{2}}+\mathit{\eta })\mathrm{\Psi }\left(\mathit{u}\right)\mathrm{d}x}}\\ & +\underset{\text{surface energy}}{\underbrace{G\underset{\mathrm{\Omega }}{\int }\left({\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathrm{\ell }}}(\mathrm{1}-s{)}^{\mathrm{2}}+\mathrm{\ell }|\mathrm{\nabla }s{|}^{\mathrm{2}}\right)\mathrm{d}x}},\end{array}\end{array}$$

where $\mathrm{0}<\mathit{\eta }\ll \mathrm{1}$ is a dimensionless residual elasticity parameter, added to prevent loss of a unique minimizing u and ellipticity. Note that ℰ is convex in u and convex in s but not jointly convex in (u,s). The two terms in Eq. (3) represent the elastic energy and the surface energy, and minimizing ℰ balances these two energy contributions. Points where s<1 reduce the elastic energy but incur a cost in the surface energy. Such points approximate discontinuous displacements u and thus model full or partial diffusive cracks. If ℓ is decreased, the surface energy puts a smaller penalty on $\left|\mathrm{\nabla }s\right|$, and thus minimizers tend to have thinner crack profiles, i.e., regions where s changes from s=1 to regions where s≈0 (corresponding to undamaged and damaged regions and cracks, respectively). Note that in practice, the phase-field function s may also take small negative values or values slightly larger than 1, which does not create any difficulties.

2.2 Staggered minimization algorithm

Following Bourdin et al. (2000), we use a staggered algorithm to solve the minimization:

$$\begin{array}{}\text{(4)}& \begin{array}{c}min\\ \mathit{u}\in U,s\in S\end{array}\mathcal{E}(\mathit{u},s),\end{array}$$

where U and S are spaces of sufficiently regular functions that satisfy appropriate Dirichlet and Neumann boundary conditions on the boundary ∂Ω. Typical Dirichlet conditions are

$$\begin{array}{}\text{(5a)}& {\displaystyle }\mathit{u}={\mathit{u}}_{\mathrm{0}}\text{on}\partial {\mathrm{\Omega }}_D\subset \partial \mathrm{\Omega },\text{(5b)}& {\displaystyle }s=\mathrm{1}\text{on}\partial \mathrm{\Omega }.\end{array}$$

In Eq. (5a), u₀ is a given displacement on a part of the boundary $\partial {\mathrm{\Omega }}_D\subset \partial \mathrm{\Omega }$. The condition (5b) is the assumption of undamaged material on the boundary and is common in phase-field models. On the remaining boundaries, natural (i.e., Neumann) conditions are assumed for u.

The algorithm to solve Eq. (4) alternates between minimizing with respect to u while keeping s fixed and minimizing with respect to s while keeping u fixed. Since ℰ is quadratic in u and s, these minimizations can be done exactly. To be more precise, for fixed s, the displacement u that minimizes ℰ must satisfy ${\mathit{\delta }}_{\mathit{u}}\mathcal{E}(\mathit{u},s)=\mathrm{0}$; i.e., the variation with respect to u vanishes. Using integration by parts, this results in the linear elasticity equation,

$$\begin{array}{}\text{(6a)}& {\displaystyle }\mathrm{\nabla }\cdot \left[\right(s^{\mathrm{2}}+\mathit{\eta }\left)\mathit{\sigma }\right]=\mathbf{0}\text{on}\mathrm{\Omega },\text{(6b)}& {\displaystyle }\mathit{u}={\mathit{u}}_{\mathrm{0}}\text{on}\partial {\mathrm{\Omega }}_D,\text{(6c)}& {\displaystyle }\mathit{\sigma }\mathit{n}=\mathbf{0}\text{on}\partial \mathrm{\Omega }\setminus \partial {\mathrm{\Omega }}_D,\end{array}$$

where n is the unit normal on the boundary. The condition (6c) results from partial integration. When u is held fixed, minimizing with respect to the phase fields s shows that the minimizer must satisfy ${\mathit{\delta }}_{\mathrm{s}}\mathcal{E}(\mathit{u},s)=\mathrm{0}$. Integration by parts shows that this is equivalent to

$$\begin{array}{}\text{(7a)}& {\displaystyle }\mathrm{\ell }\mathrm{\Delta }s-s\left({\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathrm{\ell }}}+{\displaystyle \frac{\mathrm{\Psi }\left(\mathit{u}\right)}{G}}\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathrm{\ell }}}\text{on}\mathrm{\Omega },\text{(7b)}& {\displaystyle }s=\mathrm{1}\text{on}\partial \mathrm{\Omega }.\end{array}$$

This equation shows that the elastic energy density drives the phase-field solution. In particular, if Ψ(u)≡0, then s≡1; i.e., the ice floe is undamaged. The staggered minimization algorithm is summarized as Algorithm 1. It proceeds by alternately solving Eq. (6) for fixed phase field s and solving Eq. (7) for fixed displacement u. The algorithm terminates either if the maximal update of the phase-field function is less than a given threshold ϵ or after a maximum number of iterations N.

Algorithm 1 Staggered minimization for phase-field model of fracture.

require N>0 and $\mathrm{1}\gg \mathit{ϵ}>\mathrm{0}$

set n=1, ${\mathit{ϵ}}^{*}=\mathrm{1}$, and s₀=1

while n<N and $\mathit{ϵ}<{\mathit{ϵ}}^{*}$ do

Solve Eq. (6) for u

Solve Eq. (7) for s

${\mathit{ϵ}}^{*}\leftarrow {max}_{\mathrm{\Omega }}|s-s_{\mathrm{0}}|$

s₀←s

$n\leftarrow n+\mathrm{1}$

end while

2.3 Floe-scale parameters

We introduce the ice thickness h(x)>0 into our model by assuming that the elastic strength is proportional to the ice thickness. Measurements of sea ice elastic moduli are typically not made at the kilometer scale (Timco and Weeks, 2010). In real ice, scale effects such as ice creep come into play as we transition from the lab scale to the size of a floe. In our model, we limit scale effects to those imparted by variations in ice thickness. These variations in thickness lead to variations in elastic strength, which may support concentrating stress. Hence, such variations contribute to the initiation of fracture. For this reason, relative differences in elastic strength are more important than absolute values.

We use reference Lamé parameter values from available sea ice measurements, namely, Young's modulus $E=\mathrm{9}\times {\mathrm{10}}^{\mathrm{9}}$ Pa and Poisson's ratio ν=0.3 (Timco and Weeks, 2010; Schulson and Duval, 2009; Dempsey et al., 1999). The corresponding Lamé parameters vary spatially due to ice impurities or the ice height h(x); i.e., λ and μ are given in terms of E, ν, and h(⋅) as

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathit{\lambda }\left(x\right)& =h\left(x\right){\displaystyle \frac{E\mathit{\nu }}{(\mathrm{1}+\mathit{\nu })(\mathrm{1}-\mathrm{2}\mathit{\nu })}},\\ \mathit{\mu }\left(x\right)& =h\left(x\right){\displaystyle \frac{E}{\mathrm{2}(\mathrm{1}+\mathit{\nu })}}.\end{array}\end{array}$$

For the fracture toughness we choose G=10 J m⁻² under the same scale assumptions. Meter-scale measurements of fracture energy (or apparent fracture toughness in the sea ice literature (Dempsey, 1991)) vary between 8 and 12 J m⁻². To our knowledge, no floe-scale measurements of G are available.

2.4 Discretization and solvers

Our implementation is based on the open-source finite-element software FEniCS (Alnaes et al., 2015). We use linear finite elements on triangular meshes to discretize u and the phase-field function s. The scale of fractures occurring in the model depends on ℓ, and the discretization mesh needs to be fine enough to resolve these fractures. As can be seen in Eq. (3), ℓ multiplies the norm of the gradient, and thus larger ℓ results in smaller gradients and therefore smoother phase-field functions s. The required mesh resolution can be compared with ℓ with a one-dimensional solution to Eq. (7) under simplified assumptions, as in Kuhn and Müller (2010) and Wu et al. (2020). This one-dimensional problem assumes the system is in a steady state, neglects elastic energy, and models a crack at a single point. We compare the transition between an undamaged and a damaged state so that it is well resolved by our mesh points. Mesh sizes are chosen so that the diffusive regions near cracks contain at least four or five nodes. The matrix systems occurring upon discretization of Eqs. (6) and (7) are solved using the direct solver MUMPS (Amestoy et al., 2019).

3 Fracture behavior of floes with random linear weaknesses

In this section, we use the phase-field formulation described in Sect. 2 to study the effect of varying ice thickness profiles on the fracture behavior of ice floes. First, in Sect. 3.1, we present a stochastic model to insert regions of thinner ice into our models, introducing regions for nucleation of fractures. Then, in Sect. 3.2 and 3.3, we study the statistical behavior of fractures arising in these models when the ice is subject to boundary displacements or boundary forces.

3.1 Experiment setup

For the following experiments, we use a 1 m thick, 1 km×1 km ice floe domain. Other floe geometries could have been chosen, but for simplicity we restrict our study to a square domain. Since the occurrence of fracture is inextricably linked to geometry or stress localization induced by material imperfections, we introduce random ice thickness variations into our ice floes. These features model ice impurities and other imperfections, e.g., those arising when an ice floe has gone through a fracture and refreezing cycle. Other possibilities for modeling ice thickness variations include random models of ice height (Bowen et al., 2018) or surface generation methods from random field models such as Gaussian fields (Rasmussen and Williams, 2006). Here, our focus is on fractures that occur along linear features which are not captured by those methods. We model these features directly with a simpler model.

Our stochastic model is similar to a model of random linear features, also called discrete fracture networks (DFNs) (Min et al., 2004). In these models, a configuration of K line features (here, we fix K=10) is parameterized by centers c_i, angles θ_i, and line lengths b_i, where $i=\mathrm{1},\mathrm{\dots },K$. We extend the DFN model by adding an ice height H_i. Angles, centers, and heights are drawn from uniform distributions: ${\mathit{\theta }}_i\in [\mathrm{0},\mathrm{2}\mathit{\pi }],{\mathit{c}}_i\in \mathrm{\Omega }$, and $H_{\mathrm{i}}\in [\mathrm{0},\mathrm{1}]$. Random lengths b_i are uniformly drawn from [500,1500 m]. Each line feature is spatially smoothed using a Gaussian-like mollifier into a height function h_i(x). The mollifying function is displayed in Fig. 1. It is parameterized by the minimum height H_i, a width d_m, and a Gaussian spread parameter σ>0, the latter two of which we choose as constants d_m=5 and σ=5. For a single line feature on a floe of 1 m height, the height field is given by

$$\begin{array}{}\text{(9)}& h_i\left(\mathit{x}\right)=\mathrm{1}-(\mathrm{1}-H_{\mathrm{i}})\exp \left(-{\displaystyle \frac{\mathit{\xi }(\mathit{x}{)}^{\mathrm{2}}}{\mathrm{2}\mathit{\sigma }}}\right),\end{array}$$

where

$$\mathit{\xi }\left(\mathit{x}\right)=max(\mathrm{0},d(\mathit{x})-d_{\mathrm{m}}).$$

Figure 1 Diagram of the mollifying function 9 used to define the ice height h_i(x) in the DFN model.

Here, d(x) is the distance of a point x to the line segment. This mollification prevents height discontinuities that may be difficult to resolve in computations and that might result in unphysical stresses. We set $h\left(x\right)={min}_i{h}_i\left(x\right)$; i.e., the ice thickness is the minimal thickness over the thickness fields generated from the different line segments. The resulting height field has the full height of 1 m away from the line features, is H_i along the lines, and takes the minimum height at line intersections.

Random realizations of this stochastic model can be seen in the top row in Fig. 3.

In our experiments, we assume the following boundary conditions on the right side ∂Ω_r of Ω and on the horizontal (i.e., top and bottom) sides ∂Ω_h of Ω:

$$\begin{array}{}\text{(10a)}& {\displaystyle }\mathit{u}=\mathbf{0}\text{on}\partial {\mathrm{\Omega }}_{\mathrm{r}},\text{(10b)}& {\displaystyle }\mathit{\sigma }\mathit{n}=\mathbf{0}\text{on}\partial {\mathrm{\Omega }}_{\mathrm{h}}.\end{array}$$

That is, the ice floe is held fixed on the right boundary ∂Ω_r and satisfies stress-free boundary conditions on the horizontal boundaries. On the left side ∂Ω_l, we either prescribe displacements in the normal or tangential direction (see Sect. 3.2) or apply a normal force until fracture occurs (see Sect. 3.3). An illustration of this setup is shown in Fig. 2.

In all experiments below, we use ℓ=1 and discretize the ice floe domain using 200×200 squares, each of which is split into two triangles. This results in a resolution of 5 m and about five mesh points to resolve the diffusive phase-field function in a direction orthogonal to a fracture. In all experiments in this work, we use the numerical parameters $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{6}}$ and N=5000 in Algorithm 1. With these parameters, the staggered algorithm typically terminated after between 1000 and 2000 iterations.

Figure 2 Illustration of experimental setup. The floe is fixed along the right edge, stress-free along the top and bottom edges, and is subject to normal stress or displacement along the left edge. The blue color indicates the floe thickness, where lighter regions are thinner.

3.2 Fracture arising from boundary displacement

We first study fractures that occur as a result of a prescribed tensile displacement. In these numerical experiments, we augment the displacement boundary conditions (10) with the following condition on the left boundary:

$$\begin{array}{}\text{(11)}& \mathit{u}\cdot \mathit{n}=-D,\mathit{\sigma }\mathit{n}\cdot \mathit{t}=\mathrm{0}\text{on}\partial {\mathrm{\Omega }}_{\mathrm{l}},\end{array}$$

where we choose $D:=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ m, and t is the tangential unit vector. That is, we impose displacement Dirichlet conditions on the normal direction and stress-free conditions on the tangential direction. In Fig. 3, we show examples of resulting displacement fields and fractures for various samples from our stochastic ice thickness model. We observe that fractures, if they occur, tend to follow linear features of thinner ice. Interestingly, the second crack profile from the left created a closed loop. Whether the ice floes fracture completely, for the given boundary displacement, depends on the height profile. We observe that we are in a regime where complete fractures may or may not occur.

A staggered algorithm or any method that computes critical points of Eq. (3) may terminate in local minima or saddle points when boundary conditions are too large and stress localization is absent (Bourdin et al., 2000, 2008; Amor et al., 2009). Our displacement experiments indicate that we are near the threshold of fracture, and the inserted weaknesses initiate cracks by localizing stress. We cannot ensure that Algorithm 1 finds a global minimum, but by being close to the critical stress needed for fracture to occur and by introducing inhomogeneities, we avoid two known issues that may lead to spurious minima.

In separate calculations, we have performed compression experiments by enforcing a positive normal displacement in Eq. (11). These experiments led to qualitatively similar fracture patterns to those in the tension experiments displayed in Fig. 3, so we do not discuss them further.

Figure 3 Top row: random realizations of ice thickness; brighter features correspond to thinner ice. Bottom row: corresponding fractures under tensile stress, where the phase field s=0 (red), may arise due to the given displacement of the left boundary in the leftward direction. Black arrows illustrate the displacement field u over the ice floe. Fracture does not occur for all ice thickness fields.

Next, we consider a follow-up experiment on the ice thickness realizations with shear displacements. Different from Eq. (11), we now use a shear displacement condition given by

$$\begin{array}{}\text{(12)}& \mathit{u}\cdot \mathit{t}=\stackrel{\mathrm{‾}}{D},\mathit{\sigma }\mathit{n}\cdot \mathit{n}=\mathrm{0}\text{on}\partial {\mathrm{\Omega }}_{\mathrm{l}},\end{array}$$

where $\stackrel{\mathrm{‾}}{D}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{2}}$. All other boundary conditions are as above. We observe fracture patterns of higher complexity, including fractures into multiple pieces and fractures with multiple sharp turns; see Fig. 4 for example fracture profiles and displacement fields. These results demonstrate the ability of the phase-field model to generate fracture with an intricate spatial structure under diverse loading and weakness scenarios.

Figure 4 The same as Fig. 3 but for shear experiments corresponding to the condition 12 on the left boundary.

3.3 Critical stress distribution

In our next set of experiments, we combine the boundary conditions (10) with a normal pulling force condition on the left edge ∂Ω_l of the domain; i.e.,

$$\begin{array}{}\text{(13)}& \mathit{\sigma }\mathit{n}\cdot \mathit{n}=-L,\mathit{\sigma }\mathit{n}\cdot \mathit{t}=\mathrm{0}\text{on}\partial {\mathrm{\Omega }}_{\mathrm{l}}.\end{array}$$

We increase the boundary force L≥0 until the floe completely fractures vertically. Specifically, for each floe with a realization from the stochastic ice thickness model discussed in Sect. 3.1, we increase L in steps of ΔL:=5 kN and minimize the energy in Eq. (3) using the staggered algorithm. We terminate this incremental loading schedule when we detect complete vertical fracture, and we call the corresponding value of L the critical stress value. To detect complete horizontal fracture, we check if $min\left(s\right)<s^{*}={\mathrm{10}}^{-\mathrm{3}}$ and $\Vert \mathit{u}{\Vert }_{\mathrm{\infty }}>\mathrm{1}$ m are satisfied, where $\Vert \mathit{u}{\Vert }_{\mathrm{\infty }}$ is the maximal displacement vector norm. The first condition verifies that fracture occurs at all, while the second condition detects a large displacement that occurs when the floe completely breaks vertically, as then the broken-off part is not subject to any Dirichlet conditions and is only connected to the clamped left part of the floe due to the residual stiffness $\mathrm{0}<\mathit{\eta }\ll \mathrm{1}$ in Eq. (3).

Figure 5 Histogram of critical stresses leading to vertical fracture using 1000 random samples of linear weakness fields.

Using this approach, we study the distribution of critical stresses for the various ice thickness fields. In Fig. 5, we show a histogram of the critical stress at which fracture occurs over the distribution of ice thickness fields. As can be seen, all floes fracture at or below 60 kN. Some ice floes fracture at almost an order of magnitude lower than normal forces, depending on their thickness anomalies.

Figure 6 Each dot in the scatterplot in panel (a) corresponds to a random thickness field, and the color indicates its critical stress of fracture. Shown on the x axis is the average thickness of each ice floe and on the y axis a measure for the average orientation of the linear ice features. We also show representative thickness fields corresponding to four of the samples in the scatterplot (top-left and top-right panels). The scatterplots in panels (b) and (c) show weak correlations, approximately 0.39 and −0.18, between the critical stress and average thickness and between the critical stress and α, respectively.

In Fig. 6a, we study correlations between the critical normal stress, the average ice thickness, and the orientation between the pulling force and the linear features of thinner ice in our stochastic floe model. Note that the average ice thickness (x axis) correlates with the ice strength, as can be seen by the larger number of red-shaded dots on the right of Fig. 6a, indicating larger critical stresses. This correlation is also visualized in Fig. 6b, which shows a scatterplot of critical stress vs. average thickness. The correlation is rather weak, with a corresponding correlation coefficient of 0.39, though it is statistically significant with a significance level p equal to 0, with machine precision based on a t test.

Next, we study if the angle between the force direction and the linear features of thinner ice affects the floe strength. We define the average absolute value of the sine of the angles of the linear thinner ice regions as

$$\begin{array}{}\text{(14)}& \mathit{\alpha }:={\displaystyle \frac{\mathrm{1}}{K}}\sum _{i=\mathrm{1}}^K\left|\sin \right({\mathit{\theta }}_i\left)\right|.\end{array}$$

Figure 6c shows a scatterplot of critical stress vs. α. We observe a mild tendency of ice floes to endure a larger critical stress when α is smaller, i.e., when the average angle to the force direction is smaller. The reason for this is that the displacement field largely follows the horizontal boundary force, and thus the vertical line features result in stronger stress concentration. The correlation coefficient is found to be −0.18, which is statistically significant at level $p=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{8}}$.

To further study which weaknesses are the most important for floe fracture, we compare floes with all K=10 linear features of thinner ice, as described in Sect. 3.1, with floes that only retain the most vertical feature of thinner ice. In Fig. 7, the critical stresses for the floes with all linear features are plotted against the height of the most orthogonal feature as blue-shaded dots. The shading indicates by how much the critical stress increases when the floe only contains the most vertical feature. The fact that lighter colors appear near the front of the maximal critical stress at the given height of most of the orthogonal feature indicates that the thickness of the most vertical line feature limits the stress the floe can withstand. This holds despite the fact that the most vertical feature does not generally go through the entire ice floe. We did not find any striking trends when we compared the critical stress against the orientation of the most vertical line or the height of the thinnest crack. Further, the weak correlations in Fig. 6 show that the average thickness and α do not contain significant information on critical stress.

Additionally, we study the critical stress for ice floes with a single vertical feature of thinner ice that spans the entire floe and goes through the middle. In Fig. 7, these critical stresses are plotted as a red curve against the thickness of the feature. We find that generally the critical stress of floes with the vertical feature of thinner ice is an upper bound to the critical stresses of the floes with all K=10 weaknesses.

Figure 7 Height H_i of the most vertical line (x axis) against the critical stress (y axis). Colors indicate the difference in critical stress when the floe only contains the most vertical feature versus all K=10 features. Dark-blue-shaded dots are random thickness fields where the critical stresses between the most vertical line configuration and the full network are the largest. The red line indicates the critical stress against the depth of a single vertical-spanning line.

4 Outlook and challenges

Here, we discuss the potential use of phase-field models for individual floe fractures for incorporation into large-scale DEMs and for analyzing measurements from the 2018 ICEX field campaign. The following two subsections begin by outlining the applications, followed by discussing technical challenges and potential paths forward.

4.1 Fracture in large-scale DEMs

Time-evolving DEM simulations of interacting ice floes are likely to be more accurate if they include realistic intermediate-scale fractures. Ideally, such fracture calculations should not drastically increase the computational cost of the time steps in a DEM. This creates critical hurdles to implementing fracture models for floes. Below we discuss potential paths forward in order to incorporate a physics-based fracture model, such as the phase-field model proposed here, into a DEM.

Existing DEM implementations stimulate sea ice as tightly bounded particles (Tuhkuri and Polojärvi, 2018; Jirásek and Bazant, 1995), polygons (Manucharyan and Montemuro, 2022), or polygons joined along edges (Hopkins and Thorndike, 2006). Contact forces between floes or with land masses result in stresses that may lead to fracture. Stress computations can be based on assuming constant stress over an entire floe or on more detailed floe stress fields computed, e.g., using FEM simulations. Based on these collision-caused stresses, a simulation model must decide if a floe fractures (on the floe scale, fractures are assumed to be instantaneous) and how, i.e., in how many pieces and in which directions it fractures.

Computational cost prohibits running the phase-field fracture algorithm in each time step and for each floe. In a typical DEM, each floe element requires estimating stresses due to collisions, which, at most, amounts to a partial differential equation (PDE) solve. Running the staggered phase-field solver on each element would require many PDE solves for each floe, making the straightforward use of the phase-field fracture model as part of a DEM infeasible. Nevertheless, we imagine two possible approaches to combining phase-field fracture models with a DEM: (1) speeding up the computation using a fracture dictionary generated by a large number of individual ice floe experiments, potentially interpolating between fracture experiments using machine learning, and (2) by detecting the instances when a phase-field fracture computation is crucial in order to decide whether and how a floe should fracture. In Sect. 3, we produce a sample study of fracture profiles in a fixed geometry and with fixed physical forcing, experimenting with the first approach. However, the occurrence and type of fracture likely depend strongly on floe geometry and involved forces, which would require a massive sample size. In the second approach, we envision that fracture nucleation of a single or multiple elements can be inferred from a large number of experiments. A cheap algorithm could be used to detect when a fracture computation for a floe is necessary in order to initiate the staggered phase-field solver or other algorithm.

A particular challenge would be to model leads across multiple floes, which can be observed in data. Such a behavior could for instance be achieved if one floe's fracture increases the stresses in the neighboring floe, which, as a result, also fractures. This likely requires a fracture model where fractures are initiated through stress localization and occur in the lead direction. In Sect. 3 we indicate that the phase field could use contact forces to produce crack profiles in preferred directions, depending on stresses and ice impurities.

4.2 Towards predicting fracture from ICEX 2018 data

Motivated by the high-frequency displacement data available from the ICEX 2018 expedition, we are interested in whether such data can be used to identify fractures when or even before they occur. The ICEX expedition ran from 8 to 21 March 2018 at Ice Camp Skate in Beaufort Sea, roughly 230 km north of Prudhoe Bay, Alaska. Researchers spread 24 reflectors on the ice within 1 km of the camp. The observation area contained both first-year and multi-year ice. A high-precision robotic system measured the location of the reflectors roughly every 1–3 min. Measurements were taken both before and after a crack appeared within the observation area. For a more detailed description of the data and its acquisition, we refer the reader to Parno et al. (2022).

There are several challenges to using such data to reproduce and estimate the time and shape of the physical crack. First, using position time series data in a quasi-static fracture model as proposed in this paper is not straightforward. Reflector positions need to be converted to displacements, and we do not have information on existing large-scale background stresses or displacements. Second, key data are missing. The experiments on stressed ice floes in Sect. 3 show that boundary stresses and ice impurities largely govern fracture. Thus, the low-resolution ice thickness information and the fact that only a few displacements and even fewer stress observations are available may be limiting. Narrow weaknesses could remain undetected, and stress conditions could change sharply.

Parno et al. (2022) estimate boundary displacements from the point displacement measurements under the assumption of a linear elastic model for ice deformation. Using Bayesian inference, they conclude that the linear elastic model is unlikely to explain the observations. They argue that this is due to the occurrence of a fracture that is not captured by the linear elastic mode. One could aim at replacing the linear elastic displacement model assumed in Parno et al. (2022) with the phase-field fracture model to fit observations. To do so, one could build on theoretical results regarding differentiability of objectives governed by the phase-field model (Neitzel et al., 2017), which is useful to find the best-fitting parameters. However, using the phase-field model to infer the ice state from observations would make the inference problem substantially nonlinear and thus challenging. Another issue, which could also be incorporated in the inference problem, would be to select values for the elastic moduli and fracture toughness from Sect. 2.3 that are appropriate for the spatial scales of the measurements. We believe that addressing these challenges is an interesting avenue for future work that would allow for coupling phase-field models with inference methods to predict sea ice fracture from observational data.