# A viscoelastic phase-field model for iceberg calving

Jakub Stocek  Robert J. Arthern  Oliver J. Marsh

(…)

###### Abstract

Iceberg calving accounts for around half of the ice lost annually from
Antarctica, but realistic representation of fracture and calving in large-
scale ice sheet models remains a major unsolved problem in glaciology. We
present a new phase-field viscoelastic model for fracture that simulates the
slow deformation of ice and the distribution and evolution of cracks. Cracks
nucleate and propagate in response to the evolving stress field, and are
influenced by water pressure below sea level. The model incorporates
nonlinear-viscous rheology, linear-elastic rheology, and a phase-field
variational formulation, which allows simulation of complex fracture
phenomena. We show that this approach is capable of simulating the physical
process of calving. Numerical experiments supported by a simplified model
suggest that calving rate will scale with the fourth power of ice thickness
for a floating ice front that has no variation across flow. The equations make
no assumptions about the style of calving, so they would also simulate
numerous more realistic settings in Antarctica for which material parameters
and three-dimensional effects can be expected to influence the calving rate.

##  1 Introduction

Calving of tabular icebergs from ice shelves accounts for approximately half
the ice lost from Antarctica each year (Greene et al., 2022). Failure of
calving ice cliffs is also hypothesised to make the Antarctic Ice Sheet
vulnerable to an instability that could drive rapid sea level rise (DeConto &
Pollard, 2016; Pattyn et al., 2018). Despite its central importance to the
problem of predicting global sea level, the fracture process that drives
calving, comprising both initiation and evolution of cracks, is complex and
not yet well understood. Rather than simulating the physical process of
fracture in detail, most models of calving and cliff retreat rates have been
based on heuristic arguments or limited observations scaled up to the whole of
Antarctica (DeConto & Pollard, 2016). There is no guarantee that current
methods will accurately capture the sea level contributions when conditions
deviate from present-day observations, so a focus on physically-based
modelling of fracture in glaciological settings is needed.

Glacial ice can be modelled as a viscoelastic material with Maxwell rheology
and nonlinear viscosity dependent on the strain rate (Glen, 1955). It is more
than twenty years since Meier (1997) emphasised that “iceberg calving is
largely a problem in fracture mechanics coupled to ice dynamics". This remark
was later highlighted by (Benn et al., 2007) as crucial to understanding
calving. Despite this, there have been relatively few attempts to engage fully
with the implications of this statement. Since Griffith (1921) it has been
clear that brittle fracture inherently involves the transfer of stored elastic
energy into surface energy, yet most large-scale ice sheet models neglect
elasticity in their formulation, on the assumption that timescales of interest
far exceed the Maxwell timescale at which elastic behaviour transitions to
viscous behaviour. Similarly, the literature in fracture mechanics tends to
concentrate on the elastic problem in isolation, because the timescale under
which many engineering samples undergo catastrophic failure can be considered
short enough to neglect viscous deformation. To accept Meier’s statement at
face value is to acknowledge that calving is a coupled viscoelastic fracture
problem, and requires a system of equations that involves ice viscosity,
elasticity and brittle fracture. This is the framework that we adopt here.
Following (Miehe et al., 2010a) we use phase-field approach to simulate the
brittle failure of ice.

Previously, fracture and crevasse propagation in ice have been addressed using
a combination of analytical and computational techniques. Approximations for
crevasse depths in idealised geometries have been described using a zero-
stress model advanced by Nye (1955), Meier (1958), and Weertman (1973). Under
this paradigm, crevasses are assumed to form under any tensile stress, no
matter how small, but to stop propagating at the ‘Nye depth’, where the
compressive component of stress generated by the weight of overlying ice
exceeds the tensile component, generated by longitudinal stretching. Further
extensions have been made to incorporate water filled crevasses (Benn et al.,
2007).

Computational modelling of damage evolution in ice shelves and ice sheets has
been investigated in (Bassis & Walker, 2012; Clerc et al., 2019; Lhermitte et
al., 2020; Mosbeux et al., 2023). Models based on linear elastic fracture
mechanics have been used to estimate crevasse depths on the assumption that
ice behaves elastically on short time scales (van der Veen, 1998; Lipovsky,
2020; Zarrinderakht et al., 2022). These models can make useful predictions,
but are most applicable to idealised geometries with simple boundary
conditions, and explicitly specified initial flaws. The heuristic criteria
needed to predict the onset of crack nucleation can make such models difficult
to parameterise. Implementing numerical methods that represent bifurcation and
coalescence of cracks also remains a challenge within the framework of linear
elastic fracture mechanics. Together, these factors have discouraged
widespread use of linear elastic fracture mechanics for large-scale ice sheet
modelling.

One can overcome many of the drawbacks of linear elastic fracture mechanics
with diffusive crack modelling. This uses variational approaches that are
based on energy minimisation (Francfort & Marigo, 1998). Introducing a phase
field for fracture, one effectively has a variable that interpolates between
the solid material and fracture induced voids in a sufficiently smooth manner.
This alleviates issues related to the complex crack topology. The crucial
difference for numerical treatment is the fact that the regularised problem
does not require an explicit treatment of the crack configuration. All
computations can be executed on a fixed mesh using standard finite-element
techniques. Further, regularised functionals based on phase-field formulations
are $\Gamma$–convergent to the sharp crack topology functionals for vanishing
length-scale regularisation parameter (Ambrosio & Tortorelli, 1990). In
practice, this means that the correct transfers of energy between stored
elastic energy and surface energy are approached as the regularistaion
lengthscale decreases.

Recently, phase-field models for fracture have gained a large following due to
their ability to predict complex cracking phenomena such as crack branching
and coalescence, or crack nucleation. A considerable amount of research has
been focused on brittle fracture in elastic solids (Bourdin et al., 2000;
Miehe et al., 2010a, b). New phase-field models have been developed for
dynamic fracturing, fluid-driven fracture propagation (Mikelić et al., 2015c,
b, a), as well as thermo-, visco-, elasto-, plastic materials (Miehe, 2011,
2012; Miehe et al., 2016, 2017; Miehe & Schänzel, 2014).

In this work we present a phase-field formulation of fracture for Maxwell
viscoelastic materials. This model is capable of capturing the creep of
glacial ice as well as an instantaneous elastic deformation. Phase-field
approaches have been used to predict hydrofracture in compressible elastic
glacial ice (Sun et al., 2021) as well as incompressible viscous materials
(Clayton et al., 2022). A model for viscoelastic materials with phase-field
fracture has been also discussed in (Shen et al., 2019). They focus on short-
term evolution of the fracture network in ice. In contrast, our goal is to
describe a model capable of long time evolution, crucial for modelling calving
from ice shelves.

This article is organised as follows: Section 2 discusses the viscoelastic
phase-field fracture model for ice shelf dynamics. For convenience, we divide
the presentation into three subsections. Subsection 2.1 introduces notation
for the model domain and Subsection 2.2 presents ice rheology and the
underlying equations. In Subsection 2.3 we review relevant information from
phase-field modelling of fracture, incorporate them into the rheology, and
present a non-dimensional form of the full system. In Section 3 we discuss
implementation details. Numerical experiments are presented in Section 4,
accompanied with discussion of the results in Section 5. We present concluding
remarks in Section 6. The article is accompanied by appendices and
supplementary materials that address technical details.

##  2 Model description

In this section, we describe the relevant equations governing the evolution of
ice sheets and fracture phase-field equations. We begin by describing model
geometry and its associated notation. Next we discuss the rheological
properties of ice flow and introduce their energy storage and dissipation
potential functions. We then briefly review fracture phase-field equations and
include the additional variables into the system of equations that govern ice
sheet evolution. We incorporate a hydrostatic fluid pressure condition into
the cracks. Finally, we present the nondimensionalised system of equations
that govern the evolution of the viscoelastic material and the fracture phase-
field equations.

###  2.1 Geometry

Figure 1: Illustration of geometry and associated notation.

Throughout this article we consider an intact body $\Omega_{B}$ with
inclusions $\mathcal{C}_{i}$. We introduce an extended domain
$\Omega=\overline{\Omega_{B}\cup\bigcup_{i}\mathcal{C}_{i}}$ which combines
both the intact body $\Omega_{B}$ and the inclusions $\mathcal{C}_{i}$. Denote
the boundary of $\mathcal{C}_{i}$ as
$\partial\mathcal{C}_{i}=\partial_{E}\mathcal{C}_{i}\cup\partial_{I}\mathcal{C}_{i}$,
where we decompose the boundary into exterior and interior parts with respect
to $\Omega$, respectively. We further denote boundaries of $\Omega_{B}$ and
$\Omega$ by
$\partial\Omega_{B}=\bigcup_{i}\partial_{I}\mathcal{C}_{i}\cup\partial_{E}\Omega_{B}$
and
$\partial\Omega=\bigcup_{i}\partial_{E}\mathcal{C}_{i}\cup\partial_{E}\Omega_{B}$,
respectively.

We further decompose the boundary
$\partial\Omega_{B}=\partial_{D}\Omega_{B}\cup\partial_{N}\Omega_{B}$
according to the boundary conditions. Here $\partial_{D}\Omega_{B},\
\partial_{N}\Omega_{B}$ correspond respectively to Dirichlet and Neumann,
boundary conditions. Similarly we decompose the boundary of the extended
domain $\partial\Omega=\partial_{D}\Omega\cup\partial_{N}\Omega$.

An example of a domain with three inclusions is illustrated in Figure 1. Here,
$\partial_{E}\Omega_{B}$ is the bold part of the boundary,
$\partial_{E}\mathcal{C}_{i}$ is the dashed part of the boundary, and
$\partial_{I}\mathcal{C}_{i}$ is the solid fine line of the boundary. Even
though the inclusions in Figure 1 are represented as having a volume, in case
of sharp cracks we will have to treat them as a lower dimensional object. In
that case we treat $\mathcal{C}_{i}$ and $\partial\mathcal{C}_{i}$
interchangeably. This detail is alleviated by the phase field formulation.

We also indicate the height of the domain as $H$ and the height of the water
level measured from the base as $H_{w}$. In the case of a freely floating
iceberg we note that $\partial_{D}\Omega_{B}=\emptyset$.

###  2.2 Ice Rheology

First we will focus only on the behaviour of the intact material $\Omega_{B}$
and neglect any inclusions in the material. We aim to discuss the rheological
relations governing the evolution of polycrystalline ice. Behaviour of ice
sheets and ice shelves is typically represented by a Maxwell visco-elastic
model with nonlinear stress–dependent viscosity referred to as Glen’s flow law
(Glen, 1955; Cuffey & Paterson, 2010). In the Maxwell model it is assumed that
the total strain in a material body $\Omega_{B}$ can be decomposed additively
into the elastic and viscous strains, whereas the total Cauchy stress acting
on $\Omega_{B}$ is equal to both the elastic and viscous stresses:

$\displaystyle\bm{\varepsilon}$
$\displaystyle=\bm{\varepsilon}_{e}+\bm{\varepsilon}_{v},$
$\displaystyle\bm{\bm{\sigma}}$
$\displaystyle=\bm{\bm{\sigma}}_{e}=\bm{\bm{\sigma}}_{v}.$

The elastic stress $\bm{\bm{\sigma}}_{e}$ is typically (Greve & Blatter, 2009)
related to the elastic strain $\bm{\varepsilon}_{e}$ through Hooke’s law

$\displaystyle\bm{\bm{\sigma}}_{e}=\lambda\operatorname{\mathrm{tr}}(\bm{\varepsilon}_{e})\mathds{1}+2\mu\bm{\varepsilon}_{e}.$
(1)

Here $\mathds{1}$ is the second order identity tensor and $\lambda,\mu$ are
the first and second Lamé parameters that are related to the Young’s modulus
$E$ and Poisson ratio $\nu$:

$\displaystyle\lambda=\frac{E\nu}{(1+\nu)(1-2\nu)},\quad$
$\displaystyle\mu=\frac{E}{2(1+\nu)}.$

Alternatively, we may introduce the effective pressure $p$ that is, in the
compressible case, related to the volumetric component of the elastic stress
$p=-\left(\lambda+\frac{2\mu}{3}\right)\operatorname{\mathrm{tr}}(\bm{\varepsilon}_{e})$.
Then the elastic stress is given by

$\displaystyle\bm{\bm{\sigma}}_{e}=-p\mathds{1}+2\mu\operatorname{\mathrm{dev}}(\bm{\varepsilon}_{e}),$
(2)

where
$\operatorname{\mathrm{dev}}(\bm{\varepsilon}_{e})=\bm{\varepsilon}_{e}-\frac{1}{3}\operatorname{\mathrm{tr}}(\bm{\varepsilon}_{e})\mathds{1}$
is the deviatoric part of the elastic strain tensor. The volumetric part is
the effective pressure $p$ and deviatoric part
$\operatorname{\mathrm{dev}}(\bm{\sigma}_{v})$ is governed by the Glen’s flow
law (Glen, 1955),

$\displaystyle\bm{\bm{\sigma}}_{v}=-p\mathds{1}+\operatorname{\mathrm{dev}}(\bm{\sigma}_{v})=-p\mathds{1}+2\eta(\dot{\bm{\varepsilon}}_{v})\dot{\bm{\varepsilon}}_{v}.$
(3)

The effective pressure can be expressed as
$p=-1/3\operatorname{\mathrm{tr}}(\bm{\bm{\sigma}}_{v})$. The viscous strain
rate $\dot{\bm{\varepsilon}}_{v}$ is typically given by
$\nabla_{s}\dot{\mathbf{w}}$, the symmetric gradient of $\dot{\mathbf{w}}$,
where $\dot{\mathbf{w}}$ is the rate of change of $\mathbf{w}$, and
$\mathbf{w}$ is the irrecoverable part of the total displacement $\mathbf{u}$.
The viscous stress $\bm{\bm{\sigma}}_{v}$ is decomposed into the volumetric
and deviatoric parts. In effect, Glen’s flow law (Glen, 1955) defines a
viscosity that depends upon strain rate as follows.

$\displaystyle\eta(\dot{\bm{\varepsilon}}_{v})=\frac{1}{2}{A}^{-1/n}\left(\frac{1}{2}\dot{\bm{\varepsilon}}_{v}:\dot{\bm{\varepsilon}}_{v}\right)^{\frac{1-n}{2n}},$
(4)

with $n$ usually considered to be $n=3$ and ${A}$ is a constant dependent on
temperature via the Arrhenius law (Cuffey & Paterson, 2010). Note that the
viscosity $\eta$ can be equivalently represented as a function
$\widetilde{\eta}$ with dependency on the viscous shear stress (Greve &
Blatter, 2009),

$\displaystyle\eta(\dot{\bm{\varepsilon}}_{v})=\widetilde{\eta}(\operatorname{\mathrm{dev}}(\bm{\sigma}_{v}))=\frac{1}{2}{A}^{-1}\left(\frac{1}{2}\operatorname{\mathrm{dev}}(\bm{\sigma}_{v}):\operatorname{\mathrm{dev}}(\bm{\sigma}_{v})\right)^{\frac{1-n}{2}}.$
(5)

Note that in our notation $\dot{\bm{\varepsilon}}_{v}$ is deviatoric due to
the Glen’s flow law. In certain places it may prove beneficial to write
$\operatorname{\mathrm{dev}}(\dot{\bm{\varepsilon}}_{v})$ to highlight such a
fact.

It proves useful to define the free energy and dissipation potential functions
of the visco-elastic system. The free energy function of the system is the
recoverable elastic energy:

$\displaystyle\begin{split}\psi(\bm{\varepsilon}_{e})&=\frac{1}{2}\left(\lambda\operatorname{\mathrm{tr}}(\bm{\varepsilon}_{e})^{2}+2\mu\
\bm{\varepsilon}_{e}:\bm{\varepsilon}_{e}\right)\\\
&=\frac{1}{2}\left(\left(\lambda+\frac{2\mu}{3}\right)\operatorname{\mathrm{tr}}(\bm{\varepsilon}_{e})^{2}+2\mu\operatorname{\mathrm{dev}}(\bm{\varepsilon}_{e}):\operatorname{\mathrm{dev}}(\bm{\varepsilon}_{e})\right).\end{split}$
(6)

We note that the elastic stress is given by
$\bm{\bm{\sigma}}_{e}=\partial_{\bm{\varepsilon}_{e}}\psi(\bm{\varepsilon}_{e}).$

The dissipation potential function is given by:

$\displaystyle\phi(\dot{\bm{\varepsilon}}_{v})={\frac{2n}{n+1}}\eta(\operatorname{\mathrm{dev}}(\dot{\bm{\varepsilon}}_{v}))\operatorname{\mathrm{dev}}(\dot{\bm{\varepsilon}}_{v}):\operatorname{\mathrm{dev}}(\dot{\bm{\varepsilon}}_{v}).$
(7)

The reason for defining the free energy and dissipation of the visco-elastic
system is many fold. Primarily, it allows us to incorporate the fracture
variable into the system and easily alter the system for different rheological
assumptions. Secondly, it allows us to understand thermodynamic consistency of
the system (Miehe et al., 2010b). Finally, it serves as a basis for the
numerical implementation.

The stored energy and dissipation potential functionals of the system are the
integral of the free energy and dissipation potential functions over the
domain $\Omega_{B}$ (Maugin, 1990):

$\displaystyle\mathbf{E}$
$\displaystyle=\int_{\Omega_{B}}\psi\mathop{}\\!\mathrm{d}V,$ (8)
$\displaystyle\mathbf{D}$
$\displaystyle=\int_{\Omega_{B}}\phi\mathop{}\\!\mathrm{d}V.$ (9)

The free energy stored in $\Omega_{B}$ due to the deformation is given by
$\mathbf{E}$. $\mathbf{D}$ is related to the power dissipated within the
material during its deformation.

We assume that the macroscopic motions of the body $\Omega_{B}$ are given by
the displacement field $\mathbf{u}\in\mathbb{R}^{3}$. In the small strain
context we assume that the total strain is given by the symmetric part of the
displacement gradient

$\displaystyle\bm{\varepsilon}=\nabla_{s}\mathbf{u}.$

We further assume that the total displacement $\mathbf{u}$ can be additively
decomposed into the elastic and viscous parts

$\displaystyle\mathbf{u}=\mathbf{v}+\mathbf{w}.$

Therefore,

$\displaystyle\bm{\varepsilon}=\nabla_{s}\mathbf{u}=\nabla_{s}\mathbf{v}+\nabla_{s}\mathbf{w}=\bm{\varepsilon}_{e}+\bm{\varepsilon}_{v}.$

In order to derive the system of equations governing the behaviour of the
viscoelastic material, we choose the independent constitutive state variables
to be $\mathbf{u},\mathbf{w}$. We will make a standard assumption that the
external forces act only on the total displacements (Miehe, 2011).

We could also have chosen the total and elastic displacements ($\mathbf{u}$
and $\mathbf{v}$) as state variables (instead of $\mathbf{u}$ and
$\mathbf{w}$). This would have lead to a Stokes-like system of equations.

The constitutive functions $\mathbf{E}$ and $\mathbf{D}$, defined above, are
related respectively to energy storage and dissipation due to the deformation
of the material (Miehe, 2011). The rate of energy storage at state
$\mathbf{u},\mathbf{w}$ is the time derivative of the energy functional.

$\displaystyle\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\mathbf{E}=\int_{\Omega_{B}}\delta_{\mathbf{u}}\psi(\mathbf{u},\mathbf{w})\cdot\dot{\mathbf{u}}+\delta_{\mathbf{w}}\psi(\mathbf{u},\mathbf{w})\cdot\dot{\mathbf{w}}\mathop{}\\!\mathrm{d}V+\int_{\partial_{N}\Omega_{B}}\left[\partial_{\nabla_{s}\mathbf{u}}\psi\cdot\mathbf{n}\right]\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S,$

where we introduced the functional derivatives of the free-energy function

$\displaystyle\delta_{\mathbf{u}}\psi:=-\nabla\cdot\left[\partial_{\nabla_{s}\mathbf{u}}\psi\right],\quad\delta_{\mathbf{w}}\psi:=-\nabla\cdot\left[\partial_{\nabla_{s}\mathbf{w}}\psi\right].$

Note that the functional derivatives do not contain terms with
$\partial_{\mathbf{u}},\ \partial_{\mathbf{w}}$ due to requirements of frame
invariance (Maugin, 1990; Miehe, 2011).

The internal potential $\Pi_{\mathrm{int}}$ is then composed of the
elastically-stored and dissipated contributions as follows,

$\displaystyle\Pi_{\mathrm{int}}$
$\displaystyle=\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\mathbf{E}+\mathbf{D}.$
(10)

As such, the internal potential is determined by both the energy storage
function $\psi$ and the dissipation potential function $\phi$. The effect is
to combine the influence of the total displacement $\mathbf{u}$ and the
viscous displacement $\mathbf{w}$.

We further assume that the external load functional is given by body and
surface forces $\mathbf{f},\ \mathbf{t}$ that act only on the external
variable $\mathbf{u}$

$\displaystyle\Pi_{\mathrm{ext}}$
$\displaystyle=\int_{\Omega_{B}}\mathbf{f}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V+\int_{\partial_{N}\Omega_{B}}\mathbf{t}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S.$
(11)

We define the potential $\Pi$ as the difference between the internal potential
functional and external load functional:

$\displaystyle\begin{split}\Pi=\Pi_{\mathrm{int}}-\Pi_{\mathrm{ext}}=&\int_{\Omega_{B}}\left[\delta_{\mathbf{u}}\psi-\mathbf{f}\right]\cdot\dot{\mathbf{u}}+\delta_{\mathbf{w}}\psi\cdot\dot{\mathbf{w}}+\phi\mathop{}\\!\mathrm{d}V\\\
&+\int_{\partial_{N}\Omega_{B}}\left[\partial_{\nabla_{s}\mathbf{u}}\psi\cdot\mathbf{n}-\mathbf{t}\right]\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S\\\
&+\int_{\partial_{N}\Omega_{B}}\left[\partial_{\nabla_{s}\mathbf{w}}\psi\cdot\mathbf{n}\right]\cdot\dot{\mathbf{w}}\mathop{}\\!\mathrm{d}S.\end{split}$
(12)

On thermodynamic grounds and related principles (Miehe, 2011; Maugin, 1990) we
assume that the rates of the external and internal variables at a given state
are determined by the variational principle

$\displaystyle\\{\dot{\mathbf{u}},\
\dot{\mathbf{w}}\\}=\mathrm{Arg}\\{\inf_{\dot{\mathbf{u}},\
\dot{\mathbf{w}}}\Pi(\dot{\mathbf{u}},\ \dot{\mathbf{w}})\\}.$ (13)

Taking the variation of $\Pi$ we obtain an expression for virtual rates of the
internal and external variables which satisfy homogeneous Dirichlet boundary
conditions on $\partial_{D}\Omega_{B}$:

$\displaystyle\dot{\mathbf{u}}\in\\{\mathbf{v}|\mathbf{v}=\mathbf{0}\
\mathrm{on}\
\partial_{D}\Omega_{B}\\},\quad\dot{\mathbf{w}}\in\\{\mathbf{v}|\mathbf{v}=\mathbf{0}\
\mathrm{on}\ \partial_{D}\Omega_{B}\\}.$

Applying the fundamental lemma of the calculus of variations then results in a
coupled system of equations in a domain $\Omega_{B}$ with Neumann type
boundary conditions:

$\displaystyle-\nabla\cdot\left(\lambda\nabla\cdot\left(\mathbf{u-w}\right)\mathds{1}+2\mu\nabla_{s}\left(\mathbf{u-w}\right)\right)$
$\displaystyle=\mathbf{f}$ $\displaystyle\mathrm{in}\ \Omega_{B}$ (14)
$\displaystyle
2\eta(\operatorname{\mathrm{dev}}(\nabla_{s}\dot{\mathbf{w}}))\operatorname{\mathrm{dev}}(\nabla_{S}\dot{\mathbf{w}})-2\mu\operatorname{\mathrm{dev}}\left(\nabla_{s}\left(\mathbf{u-w}\right)\right)$
$\displaystyle=0$ $\displaystyle\mathrm{in}\ \Omega_{B}$
$\displaystyle\nabla\cdot\dot{\mathbf{w}}$ $\displaystyle=0$
$\displaystyle\mathrm{in}\ \Omega_{B}$
$\displaystyle\left(\lambda\nabla\cdot\left(\mathbf{u-w}\right)\mathds{1}+2\mu\nabla_{s}\left(\mathbf{u-w}\right)\right)\cdot\mathbf{n}$
$\displaystyle=\mathbf{t}$ $\displaystyle\mathrm{on}\ \partial_{N}\Omega_{B}.$

The first equation corresponds to the momentum balance equation of the elastic
stress. The second equation corresponds to the balance of elastic and viscous
deviatoric stresses. The third equation is the standard incompressibility
condition and the fourth equation is a traction boundary condition.

The system can be also be rewritten using (2) in the momentum balance equation

$\displaystyle-\nabla\cdot\left(-p\mathds{1}+2\mu\operatorname{\mathrm{dev}}\left(\nabla_{s}\left(\mathbf{u-w}\right)\right)\right)$
$\displaystyle=\mathbf{f}$ $\displaystyle\mathrm{in}\ \Omega_{B}$ (15)
$\displaystyle
2\eta(\operatorname{\mathrm{dev}}(\nabla_{s}\dot{\mathbf{w}}))\operatorname{\mathrm{dev}}(\nabla_{s}\dot{\mathbf{w}})-2\mu\operatorname{\mathrm{dev}}\left(\nabla_{s}\left(\mathbf{u-w}\right)\right)$
$\displaystyle=0$ $\displaystyle\mathrm{in}\ \Omega_{B}$
$\displaystyle\nabla\cdot\left(\mathbf{u-w}\right)+\left(\lambda+\frac{2\mu}{3}\right)^{-1}p$
$\displaystyle=0$ $\displaystyle\mathrm{in}\ \Omega_{B}$
$\displaystyle\nabla\cdot\dot{\mathbf{w}}$ $\displaystyle=0$
$\displaystyle\mathrm{in}\ \Omega_{B}$
$\displaystyle\left(-p\mathds{1}+2\mu\operatorname{\mathrm{dev}}\left(\nabla_{s}\left(\mathbf{u-w}\right)\right)\right)\cdot\mathbf{n}$
$\displaystyle=\mathbf{t}$ $\displaystyle\mathrm{on}\ \partial_{N}\Omega_{B}.$

In the case of nearly incompressible materials, when the Poisson ratio $\nu\to
0.5$, $\lambda$ is much larger than $\mu$. This leads to a well known volume
locking phenomenon (Babuška & Suri, 1992). This can be avoided by introduction
of the pressure variable (Braess, 2007).
To ensure stability of the solution, the inf-sup stability condition needs to
be fulfilled (Braess, 2007). This means that the discretisation of the
viscoelastic system needs to be carefully chosen, see (Braess, 2007).

We assume that the external loading functions $\mathbf{f},\ \mathbf{t}$ are
the gravitational force and depth varying hydrostatic water pressure,
respectively, given by

$\displaystyle\mathbf{f}$ $\displaystyle=\rho_{s}\mathbf{g},\ \mathrm{in}\
\Omega_{B}$ (16) $\displaystyle\mathbf{t}$ $\displaystyle=-p_{w}\mathbf{n},\
\mathrm{on}\ \partial_{E}\Omega_{B}\cup\partial_{N}\Omega_{B}$ (17)

where $\rho_{s}$ is ice density, $\mathbf{g}$ is gravitational acceleration,
$\mathbf{n}$ is the unit outward pointing normal vector, and $p_{w}$ is depth
varying water pressure

$\displaystyle p_{w}=\begin{cases}\rho_{w}|\mathbf{g}|(H_{w}-z),\
&\mathrm{for}\ z<H_{w}\\\ 0,&\mathrm{for}\ z\geq H_{w}\end{cases}$ (18)

with $H_{w}=\rho_{s}/\rho_{w}H$ and $H$ being the ice-shelf thickness.
Throughout this section we focused on an intact material body without any
inclusion or cracks. Therefore, the traction boundary condition is acting only
on the exterior boundaries as no cracks exist inside of the material
$\Omega_{B}$. Water filled cracks will be incorporated into the system via the
phase field variable in Subsection 2.3.

###  2.3 Phase field fracture formulation

Before describing the details of the implementation we first highlight the
principal differences from the situation without fracture considered in the
previous section. We will briefly introduce the energetic approach to
fracture, then extend the equations from Subsection 2.2 to incorporate
fracture into the viscoelastic constitutive equations via a phase-field
variable $d$ that takes values of $d=1$ near the cracks and $d=0$ away from
the cracks. In doing so, we eliminate treatment of evolving geometry due to
fracture evolution. Instead, our equations will be defined over the whole
extended domain $\Omega$ instead of $\Omega_{B}$.

In the standard theory of brittle fracture, the drop in stored elastic energy
$G$ that occurs when a specific crack $\mathcal{C}$ is introduced into the
material, is compared to the critical energy $G_{c}$ that is needed to create
the crack. Propagation of the crack $\mathcal{C}$ occurs when $G\geq G_{c}$.
This is known as the Griffith criterion (Griffith, 1921). For pure brittle
failure $G_{c}$ is envisioned to be the surface energy that is required to
separate the crack faces, but this can also be generalised for more ductile
materials, in which additional energy must be expended to form the crack. The
resulting material parameter $G_{c}$ is commonly known as the Griffith energy
release rate and can be estimated from laboratory experiments.

In the phase-field approach, outlined in (Francfort & Marigo, 1998; Miehe et
al., 2010a; Miehe, 2011), the energy $\Gamma(\mathcal{C})$ needed to create a
crack $\mathcal{C}$ is approximated by an elliptic functional that depends
upon the critical energy release rate $G_{c}$, the phase-field variable $d$,
and its spatial gradient $\nabla d$,

$\displaystyle\Gamma_{\ell}(d)=\int_{\Omega}G_{c}\gamma(d,\nabla
d)\mathop{}\\!\mathrm{d}V.$ (19)

The crack density function $\gamma$ approximates the specific surface area of
cracks per unit volume and is defined by

$\displaystyle\gamma(d,\nabla d)=\frac{1}{2\ell}\left(d^{2}+\ell^{2}|\nabla
d|^{2}\right).$ (20)

In practice, the regularisation parameter $\ell$ acts to control the
lengthscale over which the phase-field variable $d$ varies in the neighborhood
of cracks. Larger $\ell$ corresponds to a smoother regularised transition
between fully-fractured material ($d=1$) and fully-intact material ($d=0$). It
is this regularisation that allows sharp cracks to be represented on a finite
computational mesh. Although approximate, there are theoretical reasons to
expect $\Gamma_{\ell}$ to become ever closer to the true energy $\Gamma$ as
$\ell$ is decreased (Ambrosio & Tortorelli, 1990).

To capture the release of elastic energy upon crack formation, alternative
degraded forms for the free-energy function
$\widetilde{\psi}(\bm{\varepsilon}_{e},d)$ and the dissipative potential
$\widetilde{\phi}(\dot{\bm{\varepsilon}}_{v},d)$ are used. These functions
play the same role as $\psi(\bm{\varepsilon}_{e})$ and
$\phi(\dot{\bm{\varepsilon}}_{v})$ in the previous section, but now account
for reduced ability to maintain elastic or viscous stresses within the cracked
material.

The net result is that we can define the internal potential of the system as
modified by the presence of the phase field $d$ as follows:

$\displaystyle\begin{split}\Pi_{\mathrm{int}}&=\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\widetilde{\mathbf{E}}+\widetilde{\mathbf{D}}+\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\Gamma_{\ell}\\\
&=\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\int_{\Omega}\widetilde{\psi}(\bm{\varepsilon}_{e},d)\mathop{}\\!\mathrm{d}V+\int_{\Omega}\widetilde{\phi}(\dot{\bm{\varepsilon}}_{v},d)\mathop{}\\!\mathrm{d}V+\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\int_{\Omega}G_{c}\gamma(d,\nabla
d)\mathop{}\\!\mathrm{d}V,\end{split}$ (21)

where $\widetilde{\psi}(\bm{\varepsilon}_{e},d)$ and
$\widetilde{\phi}(\dot{\bm{\varepsilon}}_{v},d)$ are the modified free-energy
density and dissipation potentials, and $\gamma(d,\nabla d)$ is the crack
density function.

In a similar fashion, we modify the external load functional, which we assume
to be given as in (11), but with an extra term that represents the effects of
pressurised water within the cracks.

$\displaystyle\Pi_{\mathrm{ext}}=\int_{\Omega}{g(d)}\mathbf{f}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V+\int_{\partial_{N}\Omega}{g(d)}\mathbf{t}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S\
{+}\int_{\Omega}p_{w}{\nabla
g(d)}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V$ (22)

In this expression, external forces have been multiplied by a degradation
function $g(d)$ that takes a value $g=1$ in fully-intact regions and $g=0$ in
fully-fractured regions. The intention is to approximate (11), in which the
external forces do not contribute within the cracked material. In line with
the literature (Miehe et al., 2010a) we choose

$\displaystyle g(d)=(1-d)^{2}.$ (23)

This is dependent on the phase-field variable $d$ and smoothly interpolates
between the undamaged state $d=0$ and a fully damaged state $d=1$. The
motivation for this choice of function is that it satisfies the following
criteria,

$\displaystyle g(0)=1,\quad g(1)=0,\quad g^{\prime}(1)=0.$

These conditions impose limits on the amount of degradation for the undamaged
and fully damaged states. The last condition ensures that the fracture force
converges to a finite value when $d\to 1$ (Miehe et al., 2010a).

The final term in (22) is motivated by an approximation that allows the
effects of water pressure on internal crack faces to be represented as a
volume integral

$\displaystyle{\int_{\cup_{i}\partial\mathcal{C}_{i}}p_{w}\mathbf{n}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S\approx-\int_{\Omega}p_{w}\nabla
g(d)\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V=\int_{\Omega}2p_{w}(1-d)\nabla
d\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V.}$ (24)

This approximation is derived in Appendix A.

Turning to the choice of $\widetilde{\psi}(\bm{\varepsilon}_{e},d)$, a
straightforward degradation of free energy $\psi(\bm{\varepsilon}_{e})$ could
be obtained simply by multiplying this functions by $g(d)$. However, this
would lead to unrealistic behaviour due to equal treatment of fracture under
tension and compression. In order to alleviate this issue we seek a
tensile–compressive decomposition so that the intact material behaves in line
with equations in Subsection 2.2, but in fully damaged parts of the material
should hold no or almost no tensile stresses.

Inspired by (Miehe et al., 2017, 2010a) we split the energy storage function
into tensile and compressive parts as follows:

$\displaystyle\widetilde{\psi}(\bm{\varepsilon}_{e},d)=g(d)\left(\psi^{+}(\bm{\varepsilon}_{e})-\psi_{crit}\right)+\left(\psi^{-}(\bm{\varepsilon}_{e})+\psi_{crit}\right).$
(25)

Here,
$\psi(\bm{\varepsilon}_{e})=\lambda/2\operatorname{\mathrm{tr}}(\bm{\varepsilon}_{e})^{2}+\mu\operatorname{\mathrm{tr}}(\bm{\varepsilon}_{e}^{2})$
corresponds to an isotropic energy function of an unbroken material and
$\psi_{crit}$ is a material parameter that acts as a crack energy threshold
(Miehe et al., 2017). Note that the crack energy threshold has no impact on
the intact material, where $d=0$. In the case of a fully damaged material,
$\psi_{crit}$ acts as a regularisation parameter that prevents a complete
degradation of the tensile energy-density function.

We consider an additive decomposition of $\psi=\psi^{+}+\psi^{-}$ where

$\displaystyle\psi^{\pm}$
$\displaystyle=\frac{1}{2}\left(\lambda+\frac{2\mu}{3}\right)\langle\operatorname{\mathrm{tr}}(\bm{\varepsilon}_{e})\rangle_{\pm}^{2}+\mu\operatorname{\mathrm{tr}}\left({\left(\operatorname{\mathrm{dev}}(\bm{\varepsilon}_{e})_{\pm}\right)^{2}}\right),$
(26)

where $\langle\cdot\rangle_{\pm}$ is the Macaulay bracket given by
$\langle\cdot\rangle_{\pm}=1/2(\cdot\pm|\cdot|)$ and $\bm{\varepsilon}_{\pm}$
is a spectral decomposition of $\bm{\varepsilon}$ into positive and negative
parts. The spectral decomposition is given by

$\displaystyle\bm{\varepsilon}_{\pm}:=\sum_{a}\langle
e_{a}\rangle_{\pm}\mathbf{m}_{a}\otimes\mathbf{m}_{a},$

where $e_{a}$ are the principal strains and $\mathbf{m}_{a}$ are the principal
directions.

Next, we define the dissipative potential
$\widetilde{\phi}(\dot{\bm{\varepsilon}}_{v},d)$. Contrary to the free energy,
we do not decompose the dissipative potential function into tensile and
compressive parts, but simply degrade the dissipative potential
$\phi(\dot{\bm{\varepsilon}}_{v})$ as follows,

$\displaystyle\widetilde{\phi}(\dot{\bm{\varepsilon}}_{v},d)=g(d)\phi(\dot{\bm{\varepsilon}}_{v})=g(d){\frac{2n}{n+1}}\eta(\operatorname{\mathrm{dev}}(\dot{\bm{\varepsilon}}_{v}))\operatorname{\mathrm{dev}}(\dot{\bm{\varepsilon}}_{v}):\operatorname{\mathrm{dev}}(\dot{\bm{\varepsilon}}_{v}).$
(27)

Equipped with the modified expressions for $\Pi_{\mathrm{int}}$ and
$\Pi_{\mathrm{ext}}$, we proceed as in Subsection 2.2, using the same
variational principle (Maugin, 1980, 1990) to derive the strong form of the
equations for the viscoelastic phase-field fracture system. Now, as well as
all admissible virtual rates $\mathbf{\dot{u},\ \dot{w}}$ we must also
consider $\Pi_{\mathrm{int}}$ and $\Pi_{\mathrm{ext}}$ to balance for all
admissible variations with respect to $\dot{d}$, the rate of change of the
phase-field variable. This provides an additional equation that must be
satisfied by the phase-field whenever $\dot{d}>0$,

$\displaystyle\frac{G_{c}}{\ell}\left(d-\ell^{2}\Delta d\right)$
$\displaystyle=2(1-d)\left(\psi^{+}-\psi_{crit}\right)$
$\displaystyle\mathrm{in}\ \Omega,$ (28) $\displaystyle\nabla
d\cdot\mathbf{n}$ $\displaystyle=0$ $\displaystyle\mathrm{on}\
\partial\Omega.$

As in most models, fracture is designed to be irreversible, with a requirement
on thermodynamic consistency (Miehe et al., 2010b). The rate of change of
$\Gamma_{\ell}(d)$ defines the dissipation of power used in creating the crack
field

$\displaystyle\dot{\Gamma}_{\ell}(d)=\int_{\Omega}G_{c}\delta_{d}\gamma(d,\nabla
d)\dot{d}\mathop{}\\!\mathrm{d}V,$

where $\delta_{d}\gamma(d,\nabla d)$ is the variational derivative of the
crack density function.

To ensure thermodynamically-consistent, irreversible fracture, we require that
the dissipated power is positive

$\displaystyle\dot{\Gamma}_{\ell}(d)\geq 0.$

This condition can be satisfied by locally imposing that

$\displaystyle\delta_{d}\gamma(d,\nabla d)$ $\displaystyle\geq 0,$ (29)
$\displaystyle\dot{d}$ $\displaystyle\geq 0.$ (30)

The second inequality is a natural constraint that assumes the phase-field
evolution is locally irreversible and does not account for any healing.

In order to satisfy the conditions (29) and (30), we follow (Miehe et al.,
2010b) and introduce a local history field of maximum tensile energy over time
or loading steps

$\displaystyle\mathcal{H}(x,t):=\max_{s}\langle\psi^{+}(\bm{\varepsilon}_{e}(x,s))-\psi_{crit}\rangle_{+}.$
(31)

The fracture phase-field evolution is then modified to be

$\displaystyle\frac{G_{c}}{\ell}\left(d-\ell^{2}\Delta d\right)$
$\displaystyle=2(1-d)\mathcal{H}$ $\displaystyle\mathrm{in}\ \Omega,$ (32)
$\displaystyle\nabla d\cdot\mathbf{n}$ $\displaystyle=0$
$\displaystyle\mathrm{on}\ \partial\Omega.$

We will use this formulation for the implementation.

For completeness, we present the full system of equations derived from the
modified expressions for $\Pi_{\mathrm{int}}$ and $\Pi_{\mathrm{ext}}$ below.
To simplify the dependence of these equations upon material parameters, we
present them in non-dimensional form. This is obtained by introducing non-
dimensional scaled quantities (labelled by asterix) defined as follows:

$\displaystyle\begin{split}&\mathbf{x}=L\mathbf{x}^{*},\quad t=\tau
t^{*},\quad\mathbf{u}=u_{c}\mathbf{u}^{*},\quad\mathbf{w}=u_{c}\mathbf{w}^{*},\quad\ell=L\ell^{*},\\\
&\lambda=\mu_{c}\lambda^{*},\quad\mu=\mu_{c}\mu^{*},\quad
p_{w}=p_{wc}p_{w}^{*},\quad\mathbf{f}=f_{c}\mathbf{f}^{*},\quad\mathbf{t}=p_{c}\mathbf{t}^{*},\\\
&\rho_{s}=\rho_{c}\rho_{s}^{*},\quad\rho_{w}=\rho_{c}\rho_{w}^{*},\quad\eta=\eta_{c}\eta^{*},\quad\bm{\bm{\bm{\sigma}}}=p_{c}\bm{\bm{\sigma}}^{*},\quad
p=p_{c}p^{*},\\\
&\mathcal{H}=\mathcal{H}_{c}\mathcal{H}^{*},\quad\psi=\mathcal{H}_{c}\psi^{*},\quad\psi^{\pm}=\mathcal{H}_{c}\psi^{\pm*},\quad\psi_{crit}=\mathcal{H}_{c}\psi_{crit}^{*}.\end{split}$
(33)

The scaling of $x,\ y,\ z$ is to be understood as the scaling of all physical
dimensions by a characteristic scale for the ice thickness $L$. Similarly,
time is scaled by $\tau$ and displacements by $u_{c}$. The other scales are
chosen in terms of $L$, $\tau$, $u_{c}$ and material constants as follows:

$\displaystyle p_{c}=\mu u_{c}L^{-1},\,\eta_{c}$
$\displaystyle=A^{-\frac{1}{n}}\left(\frac{u_{c}}{L\tau}\right)^{\frac{1-n}{n}},\,\mathcal{H}_{c}=\mu\left(\frac{u_{c}}{L}\right)^{2},\,\mu_{c}=\mu,\,f_{c}=\rho_{c}g,\,p_{wc}=\rho_{c}gL.$

Having derived the full system in exactly the same way as described previously
for the situation without fracture, we substitute these relationships, then
drop the asterix on non-dimensional quantities to provide the following non-
dimensional system of equations.

$\displaystyle-\nabla\cdot\bm{\bm{\sigma}}$
$\displaystyle=C_{1}\left[g(d)\mathbf{f}-p_{w}\nabla g(d)\right]$
$\displaystyle\mathrm{in}\ \Omega$ (34) $\displaystyle
g(d)2\eta(\operatorname{\mathrm{dev}}(\nabla_{s}\dot{\mathbf{w}}))\operatorname{\mathrm{dev}}(\nabla_{S}\dot{\mathbf{w}})-2C_{2}\operatorname{\mathrm{dev}}(\bm{\bm{\sigma}})$
$\displaystyle=0$ $\displaystyle\mathrm{in}\ \Omega$
$\displaystyle\nabla\cdot\left(\mathbf{u-w}\right)+\frac{3(1-2\nu)}{2(1+\nu)}p$
$\displaystyle=0$ $\displaystyle\mathrm{in}\ \Omega$
$\displaystyle\nabla\cdot\dot{\mathbf{w}}$ $\displaystyle=0$
$\displaystyle\mathrm{in}\ \Omega$
$\displaystyle\bm{\bm{\sigma}}\cdot\mathbf{n}$ $\displaystyle=g(d)\mathbf{t}$
$\displaystyle\mathrm{on}\ \partial_{N}\Omega$ $\displaystyle d-\ell^{2}\Delta
d$ $\displaystyle=C_{3}\ell 2(1-d)\mathcal{H}$ $\displaystyle\mathrm{in}\
\Omega$ $\displaystyle\nabla d\cdot\mathbf{n}$ $\displaystyle=0$
$\displaystyle\mathrm{on}\ \partial\Omega.$

where we have defined
$\bm{\bm{\sigma}}=g(d)\bm{\bm{\sigma}}_{e}^{+}+\bm{\bm{\sigma}}_{e}^{-}$ with
$\bm{\bm{\sigma}}_{e}^{\pm}:=\langle-p\rangle_{\pm}\mathds{1}+2\operatorname{\mathrm{dev}}\left(\nabla_{s}\left(\mathbf{u-w}\right)\right)_{\pm}$
for notational convenience. The history function $\mathcal{H}$ is defined as
in (31).

The only non-dimensional parameters that enter into the system are the Poisson
ratio $\nu$, the regularisation length $\ell^{*}$, the threshold
$\psi_{crit}^{*}$ that enters via (31), and the constants $C_{1},C_{2},C_{3}$.
These constants depend on the scales $L$, $\tau$, $u_{c}$, and other material
and physical parameters as follows:

$\displaystyle C_{1}$
$\displaystyle=\frac{L^{2}\rho_{c}|\mathbf{g}|}{u_{c}\mu},$ $\displaystyle
C_{2}$
$\displaystyle=A^{\frac{1}{n}}\left(\frac{u_{c}}{L}\right)^{1-1/n}\mu\tau^{\frac{1}{n}},$
$\displaystyle C_{3}$ $\displaystyle=\frac{\mu u_{c}^{2}}{G_{c}L}.$

These non-dimensional constants govern the operating regime of our equations.
Constant $C_{1}$ corresponds to the ratio between external and elastic
stresses, constant $C_{2}$ is the ratio between elastic and viscous stresses,
and constant $C_{3}$ is the ratio between elastic stress and fracture stress.
For increasing $C_{1}$, gravity and water pressure effects become more
important to the problem. Increasing $C_{2}$ takes the model from fully
elastic through visco-elastic to fully viscous. Increasing $C_{3}$ increases
the role of fracture from being irrelevant towards a zero-stress failure
criterion. This gives a sense of how changes in material parameters taken from
the literature (Table 1) will influence the model. As an example, for fixed
$C_{2}$, doubling the rate factor $A$ halves the characteristic timescale
$\tau$.

For our simulations we have chosen the characteristic length-scale $L=100m$,
characteristic displacement $u_{c}=10^{-2}m$, and characteristic timescale
$\tau=1a$. This leads to nondimensional constants $C_{1}=2.79,\ C_{2}=11.82,\
C_{3}=3520.75$. The nondimesional crack regularisation length-scale is chosen
as $\ell^{*}=5\times 10^{-3}$ which corresponds to a dimensional quantity of
$\ell=0.5m$.

An alternative choice of scaling, not considered further here, would be to fix
$L$ then choose scales $u_{c}$, and $\tau$ such that $C_{1}=C_{2}=1$, and
finally choose the non-dimensional regularisation length $\ell^{*}$ such that
$C_{3}\ell^{*}=1$. This would allow further simplification of the equations,
leaving just three non-dimensional parameters; the Poisson ratio $\nu$, the
regularisation length $\ell^{*}$, the threshold $\psi_{crit}^{*}$. The Poisson
ratio for ice is known approximately (Greve & Blatter, 2009), while
$\ell^{*}$, and $\psi_{crit}^{*}$ can be viewed as small parameters that
regularise the system of equations. Because $\ell^{*2}$ multiplies the highest
derivative of $d$, the problem has the character of a singular perturbation
problem.

Parameters | Values
---|---
$\rho_{s}$ | $900\ [kg/m^{3}]$
$\rho_{w}$ | $1000\ [kg/m^{3}]$
$|\mathbf{g}|$ | $9.81\ [m/s^{2}]$
$E$ | $9.33\times 10^{9}\ [N/m^{2}]$
$\nu$ | $0.325\ [-]$
$A$ | $1.2\times 10^{-25}\ [Pa^{-3}s^{-1}]$
$n$ | 3 [-]
$G_{c}$ | $1\ [N/m]$
$\psi_{crit}$ | $1\ [Pa]$
Table 1: Table of material parameters and their characteristic values. The
parameters for Young’s modulus $E$, and Poisson ratio $\nu$ are from (Greve &
Blatter, 2009), flow constant $A$ is from (Cuffey & Paterson, 2010), and
critical energy release rate $G_{c}$ is from (Goodman, 1980).

##  3 Implementation Details

In this section we briefly describe the implementation details of the
viscoelastic phase field fracture system (34). The computational domain is
decomposed into a triangular mesh and is refined in areas where cracks are
expected to propagate. In order to obtain mesh independent results, the mesh
size near the crack path needs to be chosen at least two times smaller than
the regularisation parameter $\ell$. The weak formulation is implemented using
a mixed finite element method in space, using piecewise polynomial functions
to approximate the solution. In order to ensure numerical stability of the
resulting system we use Taylor-Hood elements, where piecewise quadratic
functions are used for the displacements $\mathbf{u},\ \mathbf{w}$ and
piecewise linear functions for the pressure $p$ and phase field fracture
variable $d$ (Mang et al., 2020). In order to enforce the incompressibility
condition in (34) we include a viscous pressure-like variable $q$ that will
act as a Lagrange multiplier, and is approximated by piecewise linear
functions. This way we also ensure symmetry of the finite element system. We
approximate the time derivatives using implicit Euler time stepping method.

We note that the system (34) is nonlinear and nonconvex due to the proposed
split into tensile-compressive parts and due to the degradation function
$g(d)$. Consequently, a discretised problem of (34) would lead to a
nonsymmetric linear system. In order to alleviate this problem we propose an
extension of the staggered scheme (Miehe et al., 2010a) on an augmented system
of equations. We iteratively solve a system for the displacement-pressure, and
the fracture phase field variable using a staggered alternating minimisation
algorithm.

At each time step $k$ and at each step $i$ of the staggered iteration we
consider the phase field variable $d$ fixed and solve the following system in
order to obtain displacement and pressure variables:

$\displaystyle\begin{bmatrix}2C_{2}^{-1}/\Delta
t\bm{\tilde{K}}+2\bm{K}&-2\bm{K}&-\bm{V}^{T}&\bm{V}^{T}\\\
-2\bm{K}&2\bm{K}&\bm{0}&-\bm{V}^{T}\\\ -\bm{V}&\bm{0}&\bm{0}&\bm{0}\\\
\bm{V}&-\bm{V}&\bm{0}&-\frac{3(1-2\nu)}{2(1+\nu)}\bm{M}\\\
\end{bmatrix}\begin{bmatrix}\mathfrak{w}^{(k,i)}\\\ \mathfrak{u}^{(k,i)}\\\
\mathfrak{q}^{(k,i)}\\\
\mathfrak{p}^{(k,i)}\end{bmatrix}=\begin{bmatrix}2C_{2}^{-1}/\Delta
t\bm{\tilde{K}}\mathfrak{w}^{(k-1,I)}\\\
C_{1}\left(\bm{F}+\bm{P}\right)+\bm{T}\\\ -\bm{V}\mathfrak{w}^{(k-1,I)}\\\
\bm{0}\end{bmatrix}$ (35)

The finite element matrices and vectors are defined in Appendix B. The time
step is denoted $\Delta t$ and we take a uniform step of 1 day. The inner
iterative process continues until the system converges to a steady state and
the last iteration is denoted $I$.
We then treat the displacement and pressure variables as fixed and solve for
the fracture phase field variable ${d}$:

$\displaystyle(\bm{M}+2C_{3}\ell\hat{\bm{M}}-\ell^{2}\hat{\bm{K}})\mathfrak{d}^{(k,i)}=2C_{3}\ell\bm{\mathcal{H}}^{(k,i)}.$
(36)

The finite element matrices and vectors are again defined in Appendix B. The
history function $\bm{\mathcal{H}}^{(k,i)}$ is computed from the solution of
$\left(\mathfrak{w}^{(k,i)},\mathfrak{u}^{(k,i)},\mathfrak{q}^{(k,i)},\mathfrak{p}^{(k,i)}\right)$
as defined in (31).

Under the assumptions of small displacements and strains, linear rheology, and
load independent boundary conditions, we can usually discount differences
between the original and current configurations as the differences are
negligibly small. However, for nonlinear rheological equations and load
dependent boundary conditions as in (34) extra care needs to be taken to
obtain physically meaningful results, even under the assumptions of small
displacements and small strains (Bathe et al., 1975; Maugin, 1980).
Furthermore, evolution over long time periods may invalidate the small
displacement assumption. In such cases we need to address the differences
between the body configurations at different times.
We adopt an approach described in (Bathe et al., 1975) referred to as
material-nonlinearity-only based on the updated Lagrangian formulation, which
allows us to distinguish the body configuration at different times, and is
described in Appendix B.

##  4 Numerical Experiments

In this section we present numerical experiments to demonstrate the
capabilities of the model. Our primary focus is on tabular icebergs that float
freely in the water and are unconstrained laterally. We ignore any tidal and
drift effects.

We consider a tabular iceberg with varied thickness
$H=100,200,300,400,500,600m$ and length of $16000m$. In each run we initialise
the iceberg in a floating equilibrium with external forces switched off. We
then turn on the external forces due to gravity and water pressure. The
iceberg will initially deform elastically and then start to evolve due to
viscous creep. We initialise the icebergs either without any notches present,
with notches $5m$ in depth and one ice thickness away from the ice front, or
with notches uniformly spaced at $100m$ intervals across the whole surface of
the iceberg. Models are run for a period equivalent to $3$ years.

Figure 2 demonstrates the initial deformation and stress distribution for an
iceberg with thickness $H=300m$. Upon initiation of the model, there is a
concentration of tensile stress at the ice surface due to bending with a
maximum at a distance of just under one ice thickness from the ice front
(Reeh, 1968). Surface damage immediately appears at this stress maximum
forming a crack with depth approximately equal to the Nye depth (Nye, 1955;
Weertman, 1973).

This is a stable configuration that then evolves through creep, whereby the
ice shelf spreads, thins and bends over time. Figure 3 demonstrates the
temporal evolution of both ends of the iceberg. The initial crack length
remains constant as the presence of the crack transfers stress to the
remaining intact ice area below the crack. The stress concentration associated
with the crack tip is clearly visible just below the base of the crevasse.
During this phase of creep, the initial crack tip migrates towards the water
line. After a finite time, the crack tip crosses the water line where the
additional stress provided by water pressure provides a tensile component
sufficient to cause rapid crack propagation. This ultimate stage of failure
occurs within one time-step of the model and can be considered a type of
hydro-fracture, occurring almost instantaneously.

Figure 2: The iceberg geometry at $t=0$ days, showing the initial elastic
displacement caused by interactions between gravity and water pressure.
Displacements are exaggerated by a factor of 100.

(a)  (b)  (c)

(d)  Figure 3: Figure (a) details at $t=0$ days of iceberg geometry, principal
stress (in colour) and phase field (shading to black at $d=1$ corresponding to
fully fractured locations). Details are shown at both ends of the iceberg
(displacements are exaggerated by a factor of 100). Figures (b) and (c)
demonstrate the viscous creep of the iceberg and stress concentration near the
crack tip at intermediate times of $t=98$ and $258$ days, respectively. Figure
(d) displays both ends of the iceberg at moment of calving ($t=260$ days).

Crack location and the timing of propagation both vary strongly with the ice
thickness. The exact timing will also depend on the values of material
constants (Table 1). For simulations with $n=3$ in the Glen flow law, the time
from model initiation to calving is proportional to the ice thickness $H^{-3}$
while the distance from the crack to the ice front is proportional to $H$.
This produces an effective calving rate that is proportional to $H^{4}$.

The inclusion of initial notches in the model domain has a negligible effect
on the timing of calving (Figure 4), although naturally initiated cracks occur
slightly closer to the ice front than when notches are pre-defined at $100m$
intervals and icebergs with pre-defined notches calve slightly sooner than
naturally initiated icebergs.

(a)  (b)

(c)

Figure 4: Results for icebergs of various initial thickness (a) time to first
calving event; (b) horizontal dimension of first iceberg calved; (c) lateral
calving rate from first calving. Results are shown with and without pre-
imposed notches. Reference curves (dotted) show the predictions of a simple
model in which strain thinning allows a fixed-depth crevasse to approach the
water level, whereupon hydro-fracture occurs throughout the full depth.

##  5 Discussion

The simulations demonstrate that even brittle failure with no inherent
timescale can have its timing controlled by the slow viscous deformation of
ice. In the simulations, cracks can exist for many years in a stable
configuration before very rapidly hydro-fracturing through the full depth,
once the slow viscous thinning and bending allows their deepest points to
become submerged below the water line. At this point, a critical stress value
is reached and the crevasse propagates abruptly through the whole thickness of
the iceberg. To explain the mechanism for this rapid stage of crevasse
propagation, we note that the model assumes that crevasses flood with water
once they reach the water level. In terms of the regularised phase-field
formulation there is a softer zone ahead of the crack that can contain water
and act as an additional stress that propagates the crevasse.

To illuminate the main ideas behind the above failure mechanism, we can
consider a highly simplified model. Initial crevasse depth estimates due to
(Nye, 1955; Weertman, 1973) postulate that a crevasse propagates until a depth
$d_{0}$ that is a fraction $R_{0}$ of the iceberg thickness $h_{0}$. As noted
above, the iceberg thins and bends due to the viscous creep, while the
crevasse maintains a fixed vertical extent, so that the crevasse tip advects
towards the water line and can fill with water. Once the crevasse occupies a
certain higher fraction $R$ of the iceberg thickness $h(t)$ at a time $t$, the
crevasse will abruptly propagate towards the base of the iceberg.

Assuming no accumulation and ablation and constant density of the iceberg we
can derive a strain thinning approximation in the flow line as outlined in
(Jansen et al., 2005)

$\displaystyle\frac{\partial h}{\partial t}=-Ch^{4},$

where $C=C(A,\rho_{s},\rho_{w},\bm{g})$ is a constant dependent on material
parameters.
The above equation can be solved exactly

$\displaystyle h(t)=(h_{0}^{-3}+3Ct)^{-1/3}.$ (37)

If we require that the initial crevasse depth $d_{0}$ is equal to
$Rh(t_{crit})$ at a critical time $t_{crit}$ when the crevasse propagates:

$\displaystyle Rh(t_{crit})=d_{0}=R_{0}h_{0},$

from (37) we obtain an estimate on timing of the crevasse propagation

$\displaystyle
t_{crit}=\frac{\left(\frac{R}{R_{0}}\right)^{3}-1}{3C}h_{0}^{-3}.$ (38)

It is clear that the simplified estimate (38) cannot accurately capture the
exact timing of calving as it neglects flexure near the ice front and other
processes. However, the assertion that the time to failure $t_{crit}$ scales
with $h_{0}^{-3}$ seems appropriate, as can be seen in Figure 4(a).

Because our numerical model contains a detailed representation of the fracture
process, and because we simulate calving in such a simplified geometrical
situation (an iceberg at flotation, with no lateral variations across flow)
our simulations can serve as a simple reference case for comparison with other
calving laws.

A number of heuristic calving laws are currently under evaluation to test
their suitability as boundary conditions for large scale simulations of ice
sheets. Four different laws are compared in (Choi et al., 2018). These are 1)
height above buoyancy, 2) eigencalving, 3) crevasse-depth, and 4) von Mises
calving. It is worth briefly considering the qualitative similarities and
differences between these four calving laws and the results of our
simulations.

The height above buoyancy criterion qualitatively captures the idea that
calving rate increases with ice thickness. However, because the height above
buoyancy is zero in our simulations, this law is not appropriate for the
floating iceberg that we simulate here: the law would not predict any calving
to occur, which is in contradiction to our simulations.

The eigencalving law considers the calving rate to be proportional to the
product of two principal horizontal strain rates. When applied to the
situation that we simulate, with no lateral strain rate across flow, this
would also predict zero calving rate, so this is also in contradiction to our
results.

The von Mises criterion as implemented by Choi et al. (2018) adds positive
principal strain rates in quadrature, so this law would predict calving to
occur in our simulations. Furthermore, because the calving rate is assumed to
scale with stress it will increase with ice thickness, just as occurs in our
simulations. However, this law also has inconsistencies with our simulations.
First, the calving rate is prescribed to increase with stress, and this
increases linearly with ice thickness, rather than to the fourth power, as
occurs in our results, so this law may seriously underestimate calving rates
at large ice thickness. Secondly, the calving rate in most implementations of
this scheme is assumed to scale with the ice flow speed. This introduces a
problematic frame-dependence to the calving law. Adding a lateral translation
to all horizontal velocities would not alter the calving rate in our iceberg
simulations, so any explicit dependence of calving-rate upon ice velocity
seems problematic. In this law it seems likely that the ice speed is playing a
role as a proxy for other variables such as ice thickness.

Turning to the crevasse-depth criterion developed by Benn et al. (2007) and
used in some large-scale ice sheet models (DeConto & Pollard, 2016), we do
find some similarity in behaviour. However, even here, there is a crucial
difference. In the Benn et al. (2007) calving model a crevasse at the Nye
depth that does not reach the waterline will not produce calving. However, our
results show that a very important control on the timing of calving is the
time that it takes for the bottom of the crevasse to advect downwards towards
the water surface, whereupon sudden hydro-fracture occurs. In fact, in our
coupled viscoelastic fracture model, this is the only control on time to
failure, since we are modelling brittle failure as a rate independent process,
with no inherent timescale, and any elastic adjustments to stress can also
happen instantly. Thus, our results support the idea that the crevasse-depth
criterion is based on a sound principle, which is that hydrofracture occurs
when the deepest part of a crevasse reaches the waterline. However, there is
more to the process than this and a full treatment of the advection of
crevasses in both horizontal and vertical directions is an important component
of the calving problem.

Other calving models have been proposed in addition to those considered by
Choi et al. (2018). One example is the calving parameterisation advanced by
Crawford et al. (2021). This exhibits a power-law dependence of calving rate
on ice thickness, with an exponent in the range 6.0 to 7.3. This is
qualitatively similar behaviour to our simulations, but the exponent is higher
than the fourth power that is recovered from our simulations. Possible reasons
for this difference are that Crawford et al. (2021) consider grounded ice
fronts, while we consider floating ice fronts. Other differences are that the
workflow used by Crawford et al. (2021) transfers broad-scale geometry from a
viscous continuum model to a brittle–elastic discrete-element model, but does
not transfer the location of cracks back to the viscous model. Under this
workflow, preexisting cracks are not advected downwards towards the waterline
by the viscous flow, so the mechanism that we identify as controlling the time
to failure would not be reproduced.

In other models, the time to failure is controlled by a timescale inherent in
the fracture process itself via a damage evolution equation (Krug et al.,
2014; Mercenier et al., 2018; Pralong & Funk, 2005). Such rate dependence may
be important in ice, but it is interesting that our simulations give a
plausible explanation for the delay before calving occurs, and show that it is
not a necessary condition to have a rate dependence attached to the damage
process in order to have a finite time between the calving of successive
icebergs.

Another hypothesised mechanism for calving and cliff collapse is failure under
compressive shear (Bassis & Walker, 2012; Schlemm & Levermann, 2019). In our
simulations, failure can be generated either by a tensile component of
deviatoric elastic strain or by tensile volumetric elastic strain. This means
that failure under compressive shear can occur in our model. However, Figure 3
also shows that the interaction between the stress field, the fracture network
and the water pressure allows tensile forces to develop ahead of the crack
tip, eventually promoting full-depth failure.

The above comparisons illuminate some of the similarities and differences
between the results from the phase-field model and various other calving laws
that have been proposed. These comparisons do not represent a complete
validation of the phase-field model and a much fuller comparison with
observational datasets will be needed to assess the performance of the model
and its ability to capture the rates of calving in more realistic three-
dimensional settings found in Antarctica.

##  6 Conclusion

Importantly, the model presented here simulates the material behaviour of ice
from first principles, but was not constrained to follow any preconceived
style of calving. Despite this, it reproduces the commonly observed phenomenon
of full-depth block calving: separation of icebergs with horizontal dimension
comparable to thickness, leaving a newly-exposed vertical ice front. We have
demonstrated that the system of equations presented here can be used to model
a freely floating iceberg, numerically simulating the energetic transfers
among gravitational potential energy, stored elastic energy, irrecoverable
surface energy and dissipated heat. This is done in a thermodynamically
consistent way.

By representing the lower-dimensional crack network with a regularised phase-
field that can be evaluated on the computational mesh, a standard finite-
element approach can be used, even for complicated crack networks. The
regularisation length provides control on this regularisation, so that
realistic behaviour of cracks is recovered when it is sufficiently small
relative to the size of the simulated iceberg. The implementation of water
pressure within cracks takes full advantage of the phase field representation,
which will aid computations of the influence that pressurised water has within
more complicated crack topology.

Although we have simulated very simple geometries, there is no reason that the
same system of equations cannot be solved in three-dimensions, taking account
of the complicated stress patterns that would be generated by lateral shear,
lateral convergence or interactions between multiple cracks.

Having analysed results from the numerical model and the simple model, we can
conclude that Meier’s statement that “iceberg calving is largely a problem in
fracture mechanics coupled to ice dynamics" is indeed a useful way to address
the problem of calving.

## Funding Statement

This publication was supported by PROTECT. This project has received funding
from the European Union’s Horizon 2020 research and innovation programme under
grant agreement No 869304, PROTECT contribution number XX, and from the NERC
National Capability International grant SURface FluxEs In AnTarctica
(SURFEIT): NE/X009319/1.

## Competing interests

The authors report no conflict of interest.

## References

  * Ambrosio & Tortorelli (1990) Ambrosio, Luigi & Tortorelli, Vincenzo Maria 1990 Approximation of functional depending on jumps by elliptic functional via t-convergence. Comm. Pure Appl. Math. 43 (8), 999–1036.
  * Babuška & Suri (1992) Babuška, Ivo & Suri, Manil 1992 Locking effects in the finite element approximation of elasticity problems. Numer. Math. 62 (1), 439–463.
  * Bassis & Walker (2012) Bassis, J. N. & Walker, C. C. 2012 Upper and lower limits on the stability of calving glaciers from the yield strength envelope of ice. Proc. R. Soc. A. 468 (2140), 913–931.
  * Bathe et al. (1975) Bathe, Klaus-Jürgen, Ramm, Ekkehard & Wilson, Edward L. 1975 Finite element formulations for large deformation dynamic analysis. Int. J. Numer. Meth. Engng. 9 (2), 353–386.
  * Benn et al. (2007) Benn, Douglas I., Warren, Charles R. & Mottram, Ruth H. 2007 Calving processes and the dynamics of calving glaciers. Earth-Science Reviews 82 (3-4), 143–179.
  * Biot (1941) Biot, Maurice A. 1941 General Theory of Three-Dimensional Consolidation. Journal of Applied Physics 12 (2), 155–164, publisher: American Institute of Physics.
  * Bourdin et al. (2012) Bourdin, Blaise, Chukwudozie, Chukwudi & Yoshioka, Keita 2012 A Variational Approach to the Numerical Simulation of Hydraulic Fracturing. In All Days, pp. SPE–159154–MS. San Antonio, Texas, USA: SPE.
  * Bourdin et al. (2000) Bourdin, B., Francfort, G.A. & Marigo, J-J. 2000 Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids 48 (4), 797–826.
  * Braess (2007) Braess, Dietrich 2007 Finite elements: theory, fast solvers, and applications in elasticity theory, 3rd edn. Cambridge ; New York: Cambridge University Press.
  * Choi et al. (2018) Choi, Youngmin, Morlighem, Mathieu, Wood, Michael & Bondzio, Johannes H. 2018 Comparison of four calving laws to model Greenland outlet glaciers. The Cryosphere 12 (12), 3735–3746.
  * Clayton et al. (2022) Clayton, Theo, Duddu, Ravindra, Siegert, Martin & Martínez-Pañeda, Emilio 2022 A stress-based poro-damage phase field model for hydrofracturing of creeping glaciers and ice shelves. Engineering Fracture Mechanics 272, 108693.
  * Clerc et al. (2019) Clerc, Fiona, Minchew, Brent M. & Behn, Mark D. 2019 Marine Ice Cliff Instability Mitigated by Slow Removal of Ice Shelves. Geophysical Research Letters 46 (21), 12108–12116, _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2019GL084183.
  * Crawford et al. (2021) Crawford, Anna J., Benn, Douglas I., Todd, Joe, Åström, Jan A., Bassis, Jeremy N. & Zwinger, Thomas 2021 Marine ice-cliff instability modeling shows mixed-mode ice-cliff failure and yields calving rate parameterization. Nat Commun 12 (1), 2701, number: 1 Publisher: Nature Publishing Group.
  * Cuffey & Paterson (2010) Cuffey, Kurt & Paterson, W. S. B. 2010 The physics of glaciers, 4th edn. Burlington, MA: Butterworth-Heinemann/Elsevier, oCLC: ocn488732494.
  * DeConto & Pollard (2016) DeConto, Robert M. & Pollard, David 2016 Contribution of Antarctica to past and future sea-level rise. Nature 531 (7596), 591–597.
  * Francfort & Marigo (1998) Francfort, G.A. & Marigo, J.-J. 1998 Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids 46 (8), 1319–1342.
  * Glen (1955) Glen, J.W. 1955 The creep of polycrystalline ice. Proc. R. Soc. Lond. A 228 (1175), 519–538.
  * Goodman (1980) Goodman, D. J. 1980 Critical Stress Intensity Factor (KIc) Measurements at High Loading Rates for Polycrystalline Ice. In Physics and Mechanics of Ice (ed. Per Tryde), International Union of Theoretical and Applied Mechanics 1, pp. 129–146. Berlin, Heidelberg: Springer.
  * Greene et al. (2022) Greene, Chad A., Gardner, Alex S., Schlegel, Nicole-Jeanne & Fraser, Alexander D. 2022 Antarctic calving loss rivals ice-shelf thinning. Nature 609 (7929), 948–953.
  * Greve & Blatter (2009) Greve, Ralf & Blatter, Heinz 2009 Dynamics of Ice Sheets and Glaciers. Advances in Geophysical and Environmental Mechanics and Mathematics 1\. Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Griffith (1921) Griffith, A A 1921 VI. The phenomena of rupture and flow in solids.
  * Jansen et al. (2005) Jansen, Daniela, Sandhäger, Henner & Rack, Wolfgang 2005 Model experiments on large tabular iceberg evolution: ablation and strain thinning. J. Glaciol. 51 (174), 363–372.
  * Krug et al. (2014) Krug, J., Weiss, J., Gagliardini, O. & Durand, G. 2014 Combining damage and fracture mechanics to model calving. The Cryosphere 8 (6), 2101–2117.
  * Lhermitte et al. (2020) Lhermitte, Stef, Sun, Sainan, Shuman, Christopher, Wouters, Bert, Pattyn, Frank, Wuite, Jan, Berthier, Etienne & Nagler, Thomas 2020 Damage accelerates ice shelf instability and mass loss in Amundsen Sea Embayment. Proc. Natl. Acad. Sci. U.S.A. 117 (40), 24735–24741.
  * Lipovsky (2020) Lipovsky, Bradley Paul 2020 Ice shelf rift propagation: stability, three-dimensional effects, and the role of marginal weakening. The Cryosphere 14 (5), 1673–1683.
  * Mang et al. (2020) Mang, Katrin, Wick, Thomas & Wollner, Winnifried 2020 A phase-field model for fractures in nearly incompressible solids. Comput Mech 65 (1), 61–78.
  * Maugin (1980) Maugin, G. A. 1980 The method of virtual power in continuum mechanics: Application to coupled fields. Acta Mechanica 35 (1), 1–70.
  * Maugin (1990) Maugin, G. A. 1990 Internal Variables and Dissipative Structures. Journal of Non-Equilibrium Thermodynamics 15 (2).
  * Meier (1958) Meier, M.F. 1958 The Mechanics of Crevasse Formation. IAHS publication .
  * Meier (1997) Meier, M.F. 1997 The iceberg discharge process: observations and inferences drawn from the study of columbia glacier 15, 109–114.
  * Mercenier et al. (2018) Mercenier, Rémy, Lüthi, Martin P. & Vieli, Andreas 2018 Calving relation for tidewater glaciers based on detailed stress field analysis. The Cryosphere 12 (2), 721–739.
  * Miehe (1998) Miehe, C. 1998 Comparison of two algorithms for the computation of fourth-order isotropic tensor functions. Computers & Structures 66 (1), 37–43.
  * Miehe (2011) Miehe, Christian 2011 A multi-field incremental variational framework for gradient-extended standard dissipative solids. Journal of the Mechanics and Physics of Solids 59 (4), 898–923.
  * Miehe (2012) Miehe, Christian 2012 Mixed variational principles for the evolution problem of gradient-extended dissipative solids. GAMM-Mitteilungen 35 (1), 8–25.
  * Miehe et al. (2017) Miehe, Christian, Aldakheel, Fadi & Teichtmeister, Stephan 2017 Phase-field modeling of ductile fracture at finite strains: A robust variational-based numerical implementation of a gradient-extended theory by micromorphic regularization: PHASE FIELD MODELING OF DUCTILE FRACTURE. Int. J. Numer. Meth. Engng 111 (9), 816–863.
  * Miehe et al. (2010a) Miehe, Christian, Hofacker, Martina & Welschinger, Fabian 2010a A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering 199 (45-48), 2765–2778.
  * Miehe & Schänzel (2014) Miehe, Christian & Schänzel, Lisa-Marie 2014 Phase field modeling of fracture in rubbery polymers. Part I: Finite elasticity coupled with brittle failure. Journal of the Mechanics and Physics of Solids 65, 93–113.
  * Miehe et al. (2016) Miehe, C., Teichtmeister, S. & Aldakheel, F. 2016 Phase-field modelling of ductile fracture: a variational gradient-extended plasticity-damage theory and its micromorphic regularization. Phil. Trans. R. Soc. A. 374 (2066), 20150170\.
  * Miehe et al. (2010b) Miehe, C., Welschinger, F. & Hofacker, M. 2010b Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. Int. J. Numer. Meth. Engng. 83 (10), 1273–1311.
  * Mikelić et al. (2015a) Mikelić, Andro, Wheeler, Mary F. & Wick, Thomas 2015a A Phase-Field Method for Propagating Fluid-Filled Fractures Coupled to a Surrounding Porous Medium. Multiscale Model. Simul. 13 (1), 367–398.
  * Mikelić et al. (2015b) Mikelić, A., Wheeler, M. F. & Wick, T. 2015b Phase-field modeling of a fluid-driven fracture in a poroelastic medium. Comput Geosci 19 (6), 1171–1195.
  * Mikelić et al. (2015c) Mikelić, Andro, Wheeler, Mary F & Wick, Thomas 2015c A quasi-static phase-field approach to pressurized fractures. Nonlinearity 28 (5), 1371–1399.
  * Mosbeux et al. (2023) Mosbeux, Cyrille, Jourdain, Nicolas, Gagliardini, Olivier, Råback, Peter & Gilbert, Adrien 2023 On the effect of damage on the recent changes in the Amundsen Sea Sector. Tech. Rep.. Copernicus Meetings.
  * Nye (1955) Nye, J. F. 1955 Comments on Dr. Loewe’s Letter and Notes on Crevasses. Journal of Glaciology 2 (17), 512–514, publisher: Cambridge University Press.
  * Pattyn et al. (2018) Pattyn, Frank, Ritz, Catherine, Hanna, Edward, Asay-Davis, Xylar, DeConto, Rob, Durand, Gaël, Favier, Lionel, Fettweis, Xavier, Goelzer, Heiko, Golledge, Nicholas R., Kuipers Munneke, Peter, Lenaerts, Jan T. M., Nowicki, Sophie, Payne, Antony J., Robinson, Alexander, Seroussi, Hélène, Trusel, Luke D. & van den Broeke, Michiel 2018 The Greenland and Antarctic ice sheets under 1.5 °C global warming. Nature Clim Change 8 (12), 1053–1061, number: 12 Publisher: Nature Publishing Group.
  * Pralong & Funk (2005) Pralong, A. & Funk, M. 2005 Dynamic damage model of crevasse opening and application to glacier calving. Journal of Geophysical Research: Solid Earth 110 (B1), _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2004JB003104.
  * Reeh (1968) Reeh, Niels 1968 On The Calving of Ice From Floating Glaciers and Ice Shelves. J. Glaciol. 7 (50), 215–232.
  * Schlemm & Levermann (2019) Schlemm, Tanja & Levermann, Anders 2019 A simple stress-based cliff-calving law. The Cryosphere 13 (9), 2475–2488.
  * Shen et al. (2019) Shen, Rilin, Waisman, Haim & Guo, Licheng 2019 Fracture of viscoelastic solids modeled with a modified phase field method. Computer Methods in Applied Mechanics and Engineering 346, 862–890.
  * Sun et al. (2021) Sun, Xiangming, Duddu, Ravindra & Hirshikesh 2021 A poro-damage phase field model for hydrofracturing of glacier crevasses. Extreme Mechanics Letters 45, 101277.
  * van der Veen (1998) van der Veen, C.J. 1998 Fracture mechanics approach to penetration of surface crevasses on glaciers. Cold Regions Science and Technology 27 (1), 31–47.
  * Weertman (1973) Weertman, J. 1973 Can a water filled crevasse reach the bottom surface of a glacier? IAHS publication .
  * Zarrinderakht et al. (2022) Zarrinderakht, Maryam, Schoof, Christian & Peirce, Anthony 2022 The effect of hydrology and crevasse wall contact on calving. The Cryosphere 16 (10), 4491–4512.

##  Appendix A Inclusion of pressure inside of cracks

In the following section, we derive an approximation that allows the external
power due to the hydrostatic pressure boundary condition imposed inside of
cracks $\mathcal{C}_{i}$ to be rewritten as

$\displaystyle\int_{\cup_{i}\partial\mathcal{C}_{i}}p_{w}\mathbf{n}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S$
$\displaystyle\approx\int_{\Omega}2p_{w}(1-d)\nabla
d\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V.$ (39)

Note that $\nabla d$ can be understood as an approximation of the unit normal
$\mathbf{n}$ and the term $(1-d)$ ensures that the water pressure $p_{w}$ acts
only inside the softened zone that surrounds the crack. We note that other
approaches involving Biot’s theory of poro-elasticity (Biot, 1941) are
examined in the literature (Bourdin et al., 2012; Clayton et al., 2022;
Mikelić et al., 2015a, c).

We start with the following expression of the divergence theorem,

$\displaystyle\int_{\Omega_{B}}\nabla\cdot(p_{w}\dot{\mathbf{u}})\mathop{}\\!\mathrm{d}V$
$\displaystyle=\int_{\partial_{E}\Omega_{B}}p_{w}\mathbf{n}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S+\int_{\cup_{i}\partial\mathcal{C}_{i}}p_{w}\mathbf{n}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S,$
(40)

Note that the integral is defined only over the material part of the domain
which is time dependent due to potential fracture propagation. To avoid this
issue we approximate the volume integral using the degradation function $g(d)$
as follows:

$\displaystyle\int_{\Omega_{B}}\nabla\cdot(p_{w}\dot{\mathbf{u}})\mathop{}\\!\mathrm{d}V$
$\displaystyle\approx\int_{\Omega}g(d)\nabla\cdot(p_{w}\dot{\mathbf{u}})\mathop{}\\!\mathrm{d}V$
$\displaystyle=\int_{\Omega}g(d)\nabla
p_{w}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V+\int_{\Omega}g(d)p_{w}\nabla\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V.$

We then rewrite the second integral using the Divergence Theorem and obtain

$\displaystyle\int_{\Omega_{B}}\nabla\cdot(p_{w}\dot{\mathbf{u}})\mathop{}\\!\mathrm{d}V$
$\displaystyle\approx\int_{\Omega}g(d)\nabla
p_{w}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V-\int_{\Omega}\nabla(g(d)p_{w})\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V+\int_{\partial_{N}\Omega}g(d)p_{w}\mathbf{n}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S$
$\displaystyle=-\int_{\Omega}p_{w}\nabla
g(d)\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V+\int_{\partial_{N}\Omega}g(d)p_{w}\mathbf{n}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S.$

The last integral in (40) may be degraded in a similar fashion:

$\displaystyle\int_{\partial_{E}\Omega_{B}}p_{w}\mathbf{n}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S\approx\int_{\partial_{N}\Omega}g(d)p_{w}\mathbf{n}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S.$

Combining all of the above terms yields an approximation of the hydrostatic
pressure boundary condition imposed inside of cracks $\mathcal{C}_{i}$

$\begin{split}\int_{\cup_{i}\partial\mathcal{C}_{i}}p_{w}\mathbf{n}\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}S&\approx-\int_{\Omega}p_{w}\nabla
g(d)\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V\\\
&=\int_{\Omega}2p_{w}(1-d)\nabla
d\cdot\dot{\mathbf{u}}\mathop{}\\!\mathrm{d}V.\end{split}$ (41)

Note that $\nabla d$ can be understood as an approximation of the unit normal
$n$ and the term $(1-d)$ ensures that the water pressure $p_{w}$ acts only
inside of cracks and the surrounding softer area.

##  Appendix B Numerical implementation

In this appendix we define matrices and vectors entering into the finite
element iterative solvers from Section 3. We employ the Taylor-Hood elements
to ensure numerical stability of the finite element space. The displacements
$\mathbf{u},\mathbf{w}$ are approximated by vectorial piecewise quadratic
functions and pressure $p,\ q$ by piecewise linear functions. Fracture phase
field $d$ is approximated by piecewise linear functions.
At each time step $t_{k}$ and each iteration we have

$\displaystyle\mathbf{u}$
$\displaystyle\simeq\sum_{m}\mathfrak{u}_{m}\underline{\hat{\varphi}_{m}}(\mathbf{x}),$
$\displaystyle\mathbf{w}$
$\displaystyle\simeq\sum_{m}\mathfrak{w}_{m}\underline{\hat{\varphi}_{m}}(\mathbf{x}),$
$\displaystyle p$
$\displaystyle\simeq\sum_{j}\mathfrak{p}_{j}\varphi_{j}(\mathbf{x}),$
$\displaystyle q$
$\displaystyle\simeq\sum_{j}\mathfrak{q}_{j}\varphi_{j}(\mathbf{x}),$
$\displaystyle d$
$\displaystyle\simeq\sum_{j}\mathfrak{d}_{j}\varphi_{j}(\mathbf{x}).$

The summation is defined over all degrees of freedom stemming from the finite
element discretisation. We define $\varphi_{j}(\mathbf{x})$ to be the
piecewise linear basis functions and
$\underline{\hat{\varphi}_{m}}(\mathbf{x})$ vectorial piecewise quadratic
basis functions. Variables denoted in fraktur font are the coefficients of the
finite element basis functions.

We also need to define projection tensors that allow us to decompose the
strain tensors. Recall that a second order tensor $\mathbf{A}$ can be
decomposed

$\displaystyle\mathbf{A}=\sum_{i=1}^{3}\alpha_{i}\mathbf{m}_{i}\otimes\mathbf{m}_{i},$

where $\alpha_{i},\mathbf{m}_{i}$ are principal strains and principal
directions of $\mathbf{A}$. We define

$\displaystyle\mathbf{A}_{\pm}:=\sum_{i=1}^{3}\langle\alpha_{i}\rangle_{\pm}\mathbf{m}_{i}\otimes\mathbf{m}_{i},$

where $\langle\ \cdot\rangle_{\pm}=\frac{1}{2}(\cdot\pm|\cdot|)$ is the
MacAulay bracket. Note that $\mathbf{A}=\mathbf{A}_{+}+\mathbf{A}_{-}$. The
derivative of such decomposition then defines two fourth order projection
tensors (Miehe, 1998)

$\displaystyle\mathbb{P}_{+}$
$\displaystyle:=\partial_{\mathbf{A}}\mathbf{A}_{+}(\mathbf{A}),$
$\displaystyle\mathbb{P}_{-}$
$\displaystyle:=\partial_{\mathbf{A}}\mathbf{A}_{-}(\mathbf{A})=\mathds{1}-\mathbb{P}_{+}.$

The projection tensors for a given iteration $(k,i)$ are obtained from the
strain tensor from the previous iteration $(k,i-1)$. Therefore, we compute the
following fourth order tensors

$\displaystyle\mathbb{P}_{+}$
$\displaystyle:=\partial_{\mathbf{A}^{(k,i-1)}}\mathbf{A}_{+}(\mathbf{A}^{(k,i-1)}),$

We define the relevant matrices to solve the linear systems of equations (35)
and (B):

$\displaystyle\bm{K}^{j,l}$
$\displaystyle=\int_{\Omega}g(d^{(k,i-1)})\nabla_{s}\underline{\hat{\varphi}_{j}}:\mathbb{P}_{+}:\nabla_{s}\underline{\hat{\varphi}_{i}}+\nabla_{s}\underline{\hat{\varphi}_{j}}:\mathbb{P}_{-}:\nabla_{s}\underline{\hat{\varphi}_{l}}\mathop{}\\!\mathrm{d}V,$
(42) $\displaystyle\bm{\tilde{K}}^{j,l}$
$\displaystyle=\int_{\Omega}g(d^{(k,i-1)})\eta(\Delta
t^{-1}(\nabla_{s}\mathbf{w}^{(k,i-1)}-\nabla_{s}\mathbf{w}^{(k-1,I)}))\nabla_{s}\underline{\hat{\varphi}_{j}}:\nabla_{s}\underline{\hat{\varphi}_{i}}\mathop{}\\!\mathrm{d}V,$
$\displaystyle\bm{V}^{j,l}$
$\displaystyle=\int_{\Omega}g(d^{(k,i-1)})\langle-q^{(k,i-1)}\rangle_{+}\varphi_{j}\nabla\cdot\underline{\hat{\varphi}_{l}}+\langle-q^{(k,i-1)}\rangle_{-}\varphi_{j}\nabla\cdot\underline{\hat{\varphi}_{l}}\mathop{}\\!\mathrm{d}V,$
$\displaystyle\bm{M}^{j,l}$
$\displaystyle=\int_{\Omega}\varphi_{j}\varphi_{i}\mathop{}\\!\mathrm{d}V,$
$\displaystyle\tilde{\bm{M}}^{j,l}$
$\displaystyle=\int_{\Omega}\mathcal{H}^{(k,i)}\varphi_{j}\varphi_{l}\mathop{}\\!\mathrm{d}V,$
$\displaystyle\hat{\bm{K}}^{j,l}$
$\displaystyle=\int_{\Omega}\nabla\varphi_{j}\cdot\nabla\varphi_{l}\mathop{}\\!\mathrm{d}V.$

The external force vectors are similarly given by

$\displaystyle\bm{F}^{j}$
$\displaystyle=\int_{\Omega}g(d^{(k,i-1)})\mathbf{f}\cdot\underline{\hat{\varphi}_{j}}\mathop{}\\!\mathrm{d}V,$
(43) $\displaystyle\bm{P}^{j}$ $\displaystyle=\int_{\Omega}p_{w}\nabla
g(d^{(k,i-1)})\cdot\underline{\hat{\varphi}_{j}}\mathop{}\\!\mathrm{d}V,$
$\displaystyle\bm{T}^{j}$
$\displaystyle=\int_{\partial_{N}\Omega}g(d^{(k,i-1)})\mathbf{t}\cdot\underline{\hat{\varphi}_{j}}\mathop{}\\!\mathrm{d}S,$
$\displaystyle\bm{\mathcal{H}}^{j}$
$\displaystyle=\int_{\Omega}\mathcal{H}^{(k,i)}\varphi_{j}\mathop{}\\!\mathrm{d}V.$

As suggested in Section 3 we adopt an approach described in (Bathe et al.,
1975) referred to as material-nonlinearity-only based on the updated
Lagrangian formulation, which allows us to distinguish to the body
configuration at different times. We assume that the displacement, pressure,
and Lagrange multiplier variables at time $t_{k}$ can be decomposed as

$\displaystyle\mathbf{u}^{k}$
$\displaystyle=\mathbf{u}^{k-1}+\delta\mathbf{u}^{k},$
$\displaystyle\mathbf{w}^{k}$
$\displaystyle=\mathbf{w}^{k-1}+\delta\mathbf{w}^{k},$ $\displaystyle p^{k}$
$\displaystyle=p^{k-1}+\delta p^{k},$ $\displaystyle q^{k}$
$\displaystyle=q^{k-1}+\delta q^{k},$

where the variables $\mathbf{u}^{k-1},\ \mathbf{w}^{k-1},\ p^{k-1},\ q^{k-1}$
are assumed to be known and determine the equilibrium state at time $t_{k-1}$.
We rewrite the system (35) to be solved for variable updates
$\delta\mathbf{u}^{k},\ \delta\mathbf{w}^{k},\ \delta p^{k},\ \delta q^{k}$ so
that the small displacement assumption in each time step remains valid. We
introduce a notation for a vector $\bm{S}$ that accounts for the existing
stress from the previous time step $t_{k-1}$:

$\displaystyle\bm{S}=-\bm{V}^{T}\mathfrak{p}^{(k-1,I)}+2\bm{K}(\mathfrak{u}^{(k-1,I)}-\mathfrak{w}^{(k-1,I)}).$

The equivalent system is then given by

$\displaystyle\begin{bmatrix}2C_{2}^{-1}/\Delta
t\bm{\tilde{K}}+2\bm{K}&-2\bm{K}&-\bm{V}^{T}&\bm{V}^{T}\\\
-2\bm{K}&2\bm{K}&\bm{0}&-\bm{V}^{T}\\\ -\bm{V}&\bm{0}&\bm{0}&\bm{0}\\\
\bm{V}&-\bm{V}&\bm{0}&-\frac{3(1-2\nu)}{2(1+\nu)}\bm{M}\\\
\end{bmatrix}\begin{bmatrix}\delta\mathfrak{w}^{(k,i)}\\\
\delta\mathfrak{u}^{(k,i)}\\\ \delta\mathfrak{q}^{(k,i)}\\\
\delta\mathfrak{p}^{(k,i)}\end{bmatrix}=\begin{bmatrix}\bm{S}+\bm{V}^{T}\mathfrak{q}^{(k-1,I)}\\\
C_{1}\left(\bm{F}+\bm{P}\right)+\bm{T}-\bm{S}\\\ \bm{0}\\\
\bm{0}\end{bmatrix}$

The matrices and vectors entering into the discrete system are computed as in
(42) and (43) with the difference that the integrals are computed over the
domain configuration $\Omega_{k-1}$ at time $t_{k-1}$. As water pressure
varies with depth the integrals (43) depend on the the displacement updates
which we resolve using an inexact Newton iteration.