iceFEM: A FreeFem package for wave induced ice-shelf vibrations

Vibrations of ice-shelves in response to ocean waves were first investigated by Holdsworth & Glynn (1978) who proposed that resonant vibrations lead to icebergs calving from the shelf front. Since then seismometric measurements on the Ross ice-shelf, the largest Antarctic ice-shelf, confirmed the presence of this ocean wave-induced ice-shelf vibration (Bromirski et al., 2015; Massom et al., 2018). The period of vibration ranged from the long infragravity/tsunami waves to shorter, swell waves. More recently, Brunt et al. (2011) presented the first observational evidence that a Northern Hemisphere tsunami triggered calving on the Sulzberger ice-shelf. Mathematical models based on linear wave theory have been proposed to study these ocean-wave induced ice-shelf vibrations.


Summary
Vibrations of ice-shelves in response to ocean waves were first investigated by Holdsworth & Glynn (1978) who proposed that resonant vibrations lead to icebergs calving from the shelf front. Since then seismometric measurements on the Ross ice-shelf, the largest Antarctic ice-shelf, confirmed the presence of this ocean wave-induced ice-shelf vibration (Bromirski et al., 2015;Massom et al., 2018). The period of vibration ranged from the long infragravity/tsunami waves to shorter, swell waves. More recently, Brunt et al. (2011) presented the first observational evidence that a Northern Hemisphere tsunami triggered calving on the Sulzberger ice-shelf. Mathematical models based on linear wave theory have been proposed to study these ocean-wave induced ice-shelf vibrations. Consider the geometry and the schematic shown in Figure 1. The motion of the ice-shelf (solid) is governed by the elastodynamic equations which are coupled with the fluid motion, governed using the linear potential flow theory. The fluid is modelled as a semi-infinite rectangular region of uniform depth, whereas the sub-shelf cavity region is assumed to be non-uniform. The time-domain problem is converted to the frequency domain by applying a transformation at a prescribed frequency ω. The semi-infinite region is truncated by constructing analytic expressions and deriving a non-local boundary condition of the form The computational domain is then restricted to the sub-shelf cavity region where the finite element method is used to solve the resulting governing equations. The finite dimensional weak formulation of the coupled problem is to find where V h and W h are appropriate finite element spaces. For more details on the construction of the non-local boundary condition, refer to Ilyas et al. (2018).

Solution Method
The solution method is based on the modal expansion technique (Ilyas et al., 2018;B. Kalyanaraman et al., 2020). The displacement and the velocity potential can be written as where λ j 's are the unknown dofs. Substituting equation (5) into the weak formulation of the linear elasticity equations (3), we obtain f which corresponds to the (reduced) linear system where the entries of the matrices are given by The functions η j ∈ W h are the in-vacuo vibration modes of the ice-shelf which corresponds to solving the eigenvalue problem The diffraction potential ϕ 0 ∈ V h and the radiation potential ϕ j ∈ V h corresponding to the vibration mode η j can be obtained by solving: The diffraction and radiation potentials can be computed in parallel once the in-vacuo modes are obtained. The entries in the linear system (Equation 6) are analytic functions of the incident frequency ω. The dimension of the linear system is much smaller than the finite element degrees of freedom and can be solved efficiently. Using the analyticity of the resulting system, a large number of frequency domain solutions can be obtained by interpolating the linear system without having to solve the much larger finite element problem on a finer frequency grid.

Example
We consider the problem of simulating the vibrations of iceberg which is a modification of the ice-shelf vibration problem. This example is given in the program iceberg.edp. The following bash script solves a small set of the iceberg problem using the finite element method for different values of incident wave frequencies:  This computes the solution on a coarse ω-space, (ω c ) containing 27 points. The program then writes a set of files containing the real and imaginary part of the LHS and RHS of the linear system (6) inside the working directory along with the diffraction and radiation reflection coefficients. Once the reduced system is obtained, the entries inside the files can be interpolated on a finer ω-space (ω f ) to obtain the solution. Once the solutions for λ on the finer grid is obtained, quantities like the reflection coefficients can be computed using the new solution and the diffraction and radiation reflection coefficients (Figure 2). Several MATLAB routines are available in the modules folder to compute the solution and the reflection coefficients in real and complex ω space. For example, to perform interpolation on the real ω-space and obtain the solution λ j , the function interpolateFreq() could be used. An example script using interpolateFreq() to compute the reflection and transmission coefficients is given in RealRefIce.m. The PDF manual contains more details on the MATLAB interface along with tutorials to use the macros available within the package. More examples can be found in the README.md file located in the repository.

Statement of need
FreeFem (Hecht, 2012) is an open-source domain specific language to implement finite element methods and is based on C++. FreeFem is an excellent choice for studying ice-shelf vibrations problems due to its flexibility and ease of implementation. Similar packages include deal.II (Arndt et al., 2020), FEniCS project (Logg et al., 2012), which are examples of automated differential equation solvers and more modern packages like Gridap (Badia & Verdugo, 2020), based on Julia. iceFEM is a FreeFem package for simulating ice-shelf vibrations, heavily inspired by the ffddm module available in FreeFem for implementing the domain decomposition methods. Numerical methods have been proposed to study the vibrations of these ice-shelves, predominantly based on the thickness averaged thin-plate model for the iceshelf and the depth averaged shallow-water models for the fluid flow in the sub-shelf cavity region (Meylan et al., 2017;Sergienko, 2013). More recently, finite element methods have been proposed to model non-uniform shelf/cavity regions (Ilyas et al., 2018;B. Kalyanaraman et al., 2020). The method also involves implementing a general non-local boundary condition based on analytic expressions to handle semi-infinite regions. The numerical methods can be extended to complex valued inputs for the incident frequency and hence require flexible and versatile finite element solvers while being easy to implement. iceFEM implements the method to solve the vibration problem for complex valued incident frequencies which are useful to study resonances in the complex plane (B. Kalyanaraman et al., 2020), shown in Figure 3. The poles correspond to the resonance frequencies of the shelf/cavity system. The color denotes the phase angle and the brightness, the magnitude of the complex number (Wegert, 2012).
While a priori knowledge of finite elements are useful to write scripts using iceFEM, some macros which yield the essential stiffness matrix and load vector are available. More work is being done to improve the user-friendliness of iceFEM in future releases. Real-life examples using the BEDMAP2 dataset can also be imported and solved using iceFEM. The finite element algorithms in the iceFEM package was validated using the thin-plate solutions obtained using the eigenfunction matching methods in B. Kalyanaraman et al. (2019).