Coupling an XFEM code for linear elastic fracture propagation to a viscous flow model. GeoComputation 2019

Introduction: Ice shelves, floating extensions of grounded ice sheets and glaciers, border much of the Antarctic coastline and are responsible for most of the continent's ice mass loss. This is accomplished by melting and by iceberg calving, each process accounting for about half the total. Calving occurs as numerous smaller icebergs and by a small number of large tabular icebergs. Tabular icebergs separate from ice shelves at long, through-cutting rifts that form transverse to ice flow near the seaward fronts of the shelves. The formation and propagation of rifts transpires over time scales from days to decades, making observational studies rare. Thus, while ice shelf rifts represent a critical component in understanding and making projections of future Antarctic mass-balance, they are not well understood or represented in models. In this work we seek to adapt current day best practices in fracture propagation modelling to the Ice Sheet System Model (ISSM), therefore provide a diagnostic and predictive numerical tool for the study of ice shelf rifts.

Fracture mechanics of ice shelves: Fractures are difficult to model because they introduce complex stress fields at their tips and because they may evolve and modify the geometry of the material in which they exist. The extended finite element method (XFEM) minimizes both these difficulties by the introduction of additional sets of degrees of freedom to a classical elastic finite element method. In the XFEM, fractures can be implicitly incorporated into a continuum with a heavy-side function added to all the nodes of the elements intercut by the fracture, avoiding re-meshing at every evolution of the fracture. Asymptotic displacement functions are also used to enrich the nodes in the vicinity of the crack tip to describe the strains that LEFM studies have shown to decay in this manner from a crack tip.

In 1921, A. A. Griffith developed a theory of brittle fracture based on a balance between the strain energy at a fracture tip and the free surface energy necessary for the creation of new surfaces through propagation. This criteria is the foundation of modern linear elastic fracture mechanics (LEFM), which is now widely applied numerically in conjuncture with different continuum mechanics modelling approaches. In LEFM theory, stresses at fracture tips are described by means of scaleless stress intensity factors (SIFs). For a given material, there is a threshold SIF past which a fracture will propagate called fracture toughness. Three different SIFS, each associated to a mode of deformation relate to the stresses in the vicinity of a fracture tip and their general formulation is as follows,

\sigma = \frac{K}{\sqrt{\pi r}} f(\theta)

in which \sigma , K, r, f(\theta), are, respectively, the stress tensor, stress intensity factor, distance to the crack tip and a geometry dependent function of the angle from the crack plane.

Ice shelves are very thin relative to their horizontal dimensions (100s of m to 100s of km). It is thus typical to simplify the problem by assuming plane stress and bimodal (modes 1 and 2) crack propagation. Although there are several criteria that can be used to infer the orientation of propagation in the case of bimodal crack propagation, most criteria will promote an orientation in which the mode 2 SIF tends to zero (). This is analogous to crack tip propagation being in the direction of least extension.

Coupling method - Equivalent nodal forces:

A one-way coupling between the models is accomplished by interpolating stresses from the ice flow model to the XFEM domain. The stresses are imposed on the XFEM domain as body forces. In the standard finite element solution, nodal force vectors can be expressed as,

f_e = \int_{V_e} B^T \sigma dV

in which f_e , V_e, B, are, respectively, nodal forces, the volume occupied by said element (or area for 2D problems) and the element shape function gradients.

Preparation of the stress field: An established, open access finite element code, the Ice Sheet System Model (ISSM) is used to simulate ice shelf stresses. The coordinate system used for our ISSM model is the same as the GIS projection supporting the whole modelling process. The model geometry is established by digitizing natural features and boundaries from satellite images of the study area. Ice thickness and velocity inferred from remote sensing are used to initialise the ice flow model. Although transient simulations are possible, we only solve the instantaneous balance of stresses for this work.

preparation of the XFEM code: The XFEM implementation is an adaptation of an open access XFEM Matlab code developed by Nguyen and Bordas. This implementation is easy to use and compatible with ISSM.

A digitized rift placed within a square XFEM domain. The rift tip is located with its linear trace mid-way between the transverse extents of the domain and the tip approximately third of the way across, in the along-rift direction. This ensures that the tip remains away from the edges of the domain even after a few propagation steps. Domain size is determined by a scaling factor.

Nodal forces are evaluated for each element as described in Equation (2). Stresses are evaluated at the center of each element and as the stiffness matrix is being assembled the nodal forces are computed from these stresses and added to the force vector. Displacements are then computed using conventional matrix inversion methods.

Application: We chose the North-East corner of the Ross Ice Shelf as an application example because it is host to a number of rifts that have evolved significantly since the advent of satellite imagery (Figure 1 and 3). This region, bounded by Roosevelt Island and the Shirase Coast, has a relatively simple flow field, and thus, glaciological stress field (Figure 2). In this simple case, a single rift tip is propagated through the domain in order to demonstrate that the coupling works. More complicated geometries are also possible.

The XFEM model reproduces observed propagation geometry of the chosen ice shelf rift with high fidelity. This demonstrates that the coupling works and the models can be used to investigate a variety of issues, including propagation of complicated geometries and rift interactions. The coupled models can also be used to examine more complicated ice flow situations. It is worth noting that propagation direction of the large rift considered here was determined by the same glaciological stresses that govern the shelf motion, that is, near-field effects were not important.

Conclusion: The one-way coupling of the LEFM code to ISSM demonstrated here is new. It works well and can be applied to problems in which the ice flow field is not expected to change in response to the presence of a rift (or in which the observed field is already in equilibrium with the rift being studied). Two-way coupling would allow the ice flow (glaciological stresses) to respond to a growing rift, which would in turn feed back to the rift geometry. This is a future goal. The method described here is not limited to this medium. It could be applied to any medium where stress distribution is relatively well known and where fractures of some kind are present.