A nonconforming hydroelastic triangle for ice shelf modal analysis

Hydroelastic oscillations of ice shelves, induced by the action of ocean waves, produce deflection and stresses that could potentially lead to calving events. Due to the large horizontal span of several Antarctic ice shelves, like the Ross, Ronne or Larsen C, hydroelastic models for the ice shelf/ice shelf cavity configuration based on long wave approximations can be very effective. Such a model, based on the linearised Shallow Water Equations and the Kirchhoff–Love bending theory for slender plates is considered. For ice shelf modal analysis, in the framework of the specific model, a nonconforming hydroelastic finite element is developed. The new hydroelastic triangle is based on coupling Specht’s plate element with a linear triangle for the velocity potential approximation. It enables explicit computation of the hydroelastic coupling matrix and optimal convergence rates for the eigenpairs. The element efficiency is verified against a semi-analytical solution and the theoretically predicted convergence rates are validated for solutions with sufficient regularity. The SHEEL element (Specht HydroElastic ELement) can be used for cases of variable bathymetry and mild variations of the ice shelf thickness. The same element can be employed for time domain hydroelastic analysis with very slight modifications. A model of the Larsen C ice shelf is considered as a case study. © 2019 TheAuthors. Published by Elsevier Ltd. This is an open access article under the CCBY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Hydroelastic flexure of ice shelves, caused by the impact of ocean waves, is a challenging wave-solid interaction problem related to the integrity of Antarctic ice shelves, formation of ridges and the potential occurrence of calving events (Massom et al., 2018;Bromirski et al., 2017Bromirski et al., , 2010Brunt et al., 2011;Squire et al., 1994;MacAyeal et al., 2006;Bromirski and Stephen, 2012;Banwell et al., 2017;Lescarmontier et al., 2012). Other phenomena observed in Antarctic regions, like persistent atmospheric waves monitored over the McMurdo station at the vicinity of Ross ice shelf (Chen et al., 2016), might be related to the resonant hydroelastic oscillations of ice shelves as well (Godin and Zabotin, 2016).
Ice shelves are, fundamentally, platforms of ice, suspended on the ocean surface where glaciers or ice sheets flow down to coastal regions. The largest ice shelves can be found in Antarctica (Fig. 1 left) with Ross ice shelf being approximately 800 km across and of average thickness about 500 m. Basic features of the ice shelf geometry are the span and the water filled cavity between the ice shelf and the seabed. The boundary defined by the intersection of the floating ice shelf and the anchor ice that rests on bedrock is termed the grounding line (Fricker and Padman, 2006). Water waves impact the ice shelf at the boundary exposed to the ocean and propagate in the ice shelf/sub ice shelf cavity as hydroelastic waves. At the same time these waves undergo reflection at the incidence boundary, locations of abrupt changes in geometry, and at the grounding line as well, if their wave lengths permit propagation up to this point and attenuation is not strong (Fig. 1  right). Typically, pulses with larger wavelengths propagate further inside the ice shelf/sub ice shelf cavity region and are furthermore less influenced by the buffering effect of the Marginal Ice Zone (MIZ) (Zhang et al., 2015;Montiel and Squire, 2017;Squire and Williams, 2008;Meylan et al., 2018). For these reasons, hydroelastic analysis due to the impact of long waves (ranging from tidal waves to long period swells) is particularly important for ice shelf integrity assessment.
Advanced modelling techniques for the interaction of waves with ice shelves have been presented by several authors Holdsworth (1969), Vinogradov and Holdsworth (1985), Vaughan (1995), Sergienko (2010), Papathanasiou et al. (2014), Papathanasiou et al. (2015), Ilyas et al. (2018) and Mattsson et al. (2018). One of the proposed approaches regarding long wave excitation is the coupling of linearised Shallow Water Equations (SWE) with floating cantilevers, modelled as Euler-Bernoulli beams. Both assumptions, i.e. shallow water conditions and Euler-Bernoulli approximation for bending can be justified, in most cases, due to the large horizontal extent of the spatial domains involved and the wave conditions affecting the response of the ice shelf. Regarding the bending hypothesis, although the thickness of ice shelves can be several hundreds of metres, the horizontal span is typically much larger (several kilometres) and hence the contribution of shear deformation energy is much less than that of bending. Sergienko (2013) and Meylan et al. (2017) applied such a model for the calculation of ice shelf eigenperiods. Papathanasiou et al. (2019) solved the same problem using a finite element method and proposed the use of a different homogeneous boundary condition at the sub ice shelf cavity opening to obtain eigenperiods that match the characteristic values maximising the flexural response of the ice-shelf. However, most of the models and solution techniques presented so far refer to configurations in one horizontal dimension, that is along transects from the grounding line to the ocean exposed edge. Studies in two horizontal dimensions include the analytical results in Godin and Zabotin (2016) and a 3D finite element model in COMSOL, regarding the hydroelastic response of the Ross ice shelf, presented recently by Sergienko (2017).
Due to the irregular shape of ice shelves, it is expected that realistic modelling will depend heavily on the availability of numerical techniques for the governing set of coupled hydroelastic equations, involving two horizontal dimensions. At the same time, due to the high slenderness ratio of large Antarctic ice shelves, 3D finite element models might have increased computational cost, particularly if parametric studies are involved. Reduced dimensionality approaches, combining plate theories and long wave approximation models could therefore offer an optimum balance of accuracy and cost effectiveness. The aim of the present study is to introduce a hydroelastic finite element scheme for the coupled problem of a slender plate interacting with long waves. In particular, the linearised SWE, coupled with the Kirchhoff-Love model for a floating plate that is clamped at a part of the boundary are considered. The envisioned finite element is expected to efficiently model the hydroelastic response of large ice shelves for long wavelength excitations.
The study is organised as follows: The coupled system of governing equations along with appropriate boundary conditions to represent the grounding line and the ice shelf front exposed to the ocean is presented in Section 2. In the same section, the continuous variational form of the ice shelf modal analysis problem is presented. Next, in Section 3, the discrete variational form is considered and the proposed hydroelastic element for the discrete eigenvalue problem is introduced. Several attributes of the finite element are analysed and the exact form of the local hydroelastic coupling matrix is derived and included in Appendix A. The computation of eigenfrequencies and eigenmodes is based on a resulting quadratic eigenvalue problem and is analysed in Section 4. Properties of the global mass, stiffness and hydroelastic coupling matrices are discussed. An added mass formulation for the same problem is presented as well and used to derive a priori error estimates and asymptotic convergence rates for the eigenpairs. In Section 5 verification of the method against a semi analytical-solution detailed in Appendix B is presented, and several numerical experiments are presented and discussed, demonstrating the effectiveness of the proposed Finite Element. Moreover, the theoretically derived convergence rates are verified and the effects of variable bathymetry and thickness are briefly discussed.
Finally, in Section 6 a model for the simulation of the Larsen C ice shelf and comparisons with previous less sophisticated models (Papathanasiou et al., 2019) are included. The paper concludes with some remarks on the effectiveness of the newly derived hydroelastic finite element and suggestions for future developments regarding its applicability in the analysis of ice shelf vibrations.

The hydroelastic ice shelf model
Let Ω ⊂ R 2 be closed and bounded with boundaryΓ . This domain is used to represent the horizontal projection of the ice shelf/ice shelf cavity configuration. The bathymetry function b(x, y): Ω → R + and the positive constants ρ, ρ i , with ρ < ρ i , denoting the density of water and the density of each plate respectively, are introduced. According to Archimedes principle, the draft of the ice shelf is d = ρ i τ /ρ, where τ = τ (x, y) is the thickness distribution of the ice shelf. The domain occupied by the fluid is thus defined as (2.1) The basic assumptions introduced are: (A1) The fluid motion is irrotational and governed by linear shallow water hydrodynamics. This assumption is compatible with the long wave excitation of ice shelves, due to the very large horizontal dimensions involved.
(A2) The ice shelf is a slender, elastic solid, subject to the Kirchhoff-Love assumptions for bending of thin plates. In particular, it is τ /diam(Ω) ≪ 1. The plate material is an isotropic solid satisfying Hooke's law. Its elastic modulus is E > 0 and its Poisson's ratio ν ∈ (0, 0.5]. Similar assumptions apply to the case of Very Large Floating Structures under long wave excitation as well (Squire, 2008;Papathanasiou and Belibassakis, 2018). The scalar function Φ: Ω × (0, T ] → R, representing the fluid velocity potential is considered. The velocity vector, with velocity components along the x and y direction, in each region is then defined as u = ∇Φ, (x, y) ∈ Ω. (2.2) Denoting the ice shelf deflection in Ω as η: Ω × (0, T ] → R, continuity of fluid mass underneath the floating plate, leads to the equation In order to formulate the hydroelastic interaction model, Eq. (2.3) must be supplemented by a plate deflection model. Assuming that the ice shelf is always in contact with water, the fluid upper surface elevation coincides with the floating plate deflection. Thus, the pressure on the upper surface of the fluid acts as a distributed load on the plate. If no abrupt variations in the thickness of the ice shelf occur, the formulae for bending and twisting moments of plates with constant thickness can still be used with sufficient accuracy (see also Timoshenko and Woinowsky-Krieger, 1959, p. 173). The assumption of small variations in the thickness of the ice shelf can be justified in most cases. It is linked to the large horizontal dimensions of the ice shelf, as changes in the thickness occur within several kilometres and the slopes are generally mild. Hence, the expressions for the bending and twisting moments are 12(1−ν 2 ) . Considering mild slopes for the ice shelf profile and assumptions (A1) and (A2), the momentum equilibrium equation governing the response of the ice shelf is Freely floating boundary conditions dictate imposing the normal bending moment and active shear force to be zero at the lateral boundary of the ice shelf exposed to the ocean. These boundary conditions have the form M n (η) = n T M(η)n = 0 and Q(η) + ∇T(η) · t = 0, on Γ f , where T(η) = n T M(η)t is the twisting moment and Q(η) = ∇ · M(η)n is the shear force, with t being the unit tangent vector on Γ and M(η) is the bending moment matrix.
Since different boundary conditions need to be applied at the grounding line and the ice shelf/ocean front, the boundary Γ of Ω is the union Γ = Γ g ∪ Γ f , where Γ g = {(x, y) ∈ Γ | η = 0, ∇η · n = 0 and ∇Φ · n = 0} , and Γ g ∩ Γ f ≡ ∅ (see Fig. 2). Following Papathanasiou et al. (2019), a homogeneous Dirichlet condition for the velocity potential is applied on Γ f . This boundary condition was found to produce an eigenvalue problem with eigenfrequencies that closely match the characteristic frequencies of the more realistic reflection/transmission formulation for ice shelf resonant behaviour models (see also Fig. 1b). Although the latter work was restricted to 1D wave propagation, the extension of the present method to 2D horizontal direction to treat incident wave from various directions and backscattering effects is possible by truncating the water subdomain and implementing appropriate nonreflecting boundary conditions at far distances from the elastic body. This is left for future work.
Introducing a characteristic length L, e.g. L = diam (Ω), the following nondimensional quantities are defined Introducing now standing wave solutions of the form η(x, , where i 2 = −1, and using nondimensional quantities the governing equations, in steady state conditions, become (after dropping tildes) The variational form of this problem can be derived by multiplying Eqs. (2.10), and (2.11) by the conjugate of suitable weight functions χ, −w respectively, integration over Ω and application of the Green-Gauss theorem. Regarding the fourth order term ∆ (K ∆η), two repeated applications of the Green-Gauss theorem and appropriate handling of the resulting boundary integrals, using the surface and normal gradient operators (Timoshenko and Woinowsky-Krieger, 1959;Ciarlet, 2002), the surface divergence theorem and Stokes theorem, are required. The weak form of (2.10) and (2.11) is where an superimposed bar indicates complex conjugate and a(χ , η) = ∫ Ω K ν∆χ∆ηdxdy + ∫ Ω K (1 − ν)(∂ xx χ∂ xx η + ∂ yy χ∂ yy η + 2∂ xy χ∂ xy η)dxdy . (2.14) Introducing now the function spaces associated with the weak solution of the ice shelf/ice shelf cavity eigenstate problem and taking into account the homogeneous natural boundary conditions in (2.6), the variational form of the eigenvalue problem can be stated as (2.17) where the standard notation H k (Ω ; C) is used for the Sobolev space W k,2 (Ω; C) of complex-valued functions f defined on Ω that have Lebesgue-square integrable derivatives ∂ a 1 +a 2 f ∂x a 1 ∂y a 2 for all non-negative integers a 1 +a 2 ≤ k ∈ N 0 . In the following, the notation ∥·∥ k is adopted for the inner product induced norm in H k (Ω ; C), k ≥ 1 while the symbol |·| k denotes the corresponding semi norm. For k = 0, the notation (·, ·) and ∥·∥ will be adopted for the L 2 (Ω ; C) inner product and the induced norm respectively.

The nonconforming hydroelastic triangle
Let T h denote a quasi-uniform triangulation of Ω, where h . = max h K , h K being the characteristic size of triangle K. A nonconforming approximation is considered for η. In particular, the Specht (1988) triangle shape functions are selected, and the solution η h of the discretised problem is such that The restriction of η h inside each element K is a polynomial of degree 3 and the bilinear functional associated with plate bending is (3.1) Conforming approximation using linear shape functions in each triangle K is adopted for the velocity potential, and in The structure of the nonconforming hydroelastic triangle and the form of the finite element matrices in element level will be analysed in the following part of this section.
The K triangle vertices (x i , y i ) are numbered locally and counter-clockwise, as shown in Fig. 3, and ℓ i , A K denote the length of the edge opposite to node i and the triangle area respectively. Furthermore the parameters µ i , b i , c i are introduced, defined as and i, j, k are cyclic permutations of 1, 2, 3.
Linear shape functions, represented in area coordinates as L i (x, y) = A i (x, y)/A K , i = 1, 2, 3 are employed for the interpolation of the velocity potential Φ. The restriction of the finite element solution Φ h in element K is therefore 3 are the nodal velocity potential values at the triangle vertexes. In the following, the acronym SHEEL (Specht HydroElastic ELement) will be used when referring to this element.
For the approximation of the ice shelf deflection, a nonconforming set of shape functions, namely the shape functions introduced by Specht (1988) for triangular thin plate elements, will be used. The restriction of the finite element solution where the nodal degrees of freedom are y denote the rotations with respect to x and y axis respectively, and the nine Specht shape functions are with i, j, k cyclic permutations of 1, 2, 3. Arranging the nodal unknowns inside element K as the matrix of the hydroelastic shape functions takes the form (3.12c) Table 1 Gauss points and weights for the bending stiffness integrals in b (as suggested in Ref. Specht, 1988).

Point
Area coordinates Weight where the b term corresponds to the flexural energy, p corresponds to potential energy and v term corresponds to water kinetic energy. Obviously it is where (ic K ) * is the conjugate transpose of ic K . Although several of the integrals above can be calculated in closed form, assuming constant depth and uniform material properties, Gauss quadratures for triangles need to be employed in the general case. Exact integration is always possible for c K , the hydroelastic coupling term though. The closed form of c K is given in Appendix A.
Exact formulae for the integrals in m K for cases where M (or equivalently the ice shelf thickness) is constant can be found in the appendix of reference (Dmitrochenko and Mikkola, 2008). The same formulae apply to the calculation of For the calculation of b, the stiffness matrix term corresponding to flexural energy, in the original paper of Specht, a Gauss quadrature with three points was suggested as the optimal choice. In particular, using area coordinates and a reference triangle, the integration rule proposed in Specht (1988), and used in the present study as well, is summarised in Table 1. Integral v involves the standard linear shape functions for triangles and in particular their derivatives. For cases of variable bathymetry Gauss integration can be applied.

Computation of eigenpairs and error estimates
Assembly of the element contributions in global matrices leads to a quadratic eigenvalue problem of the form where the vector of unknowns is U = [ η θ Φ ] T and the global mass, hydroelastic coupling and stiffness matrices have the form: The element definition in these three matrices, according to Eqs. (3.12a)-(3.12c), implies that: (i) Matrix M is symmetric, and positive semi-definite. In particular this matrix is singular as the velocity potential does not appear under the second derivative sign in the governing equations. Regarding the number of zero eigenvalues in M, it is equal to the total number of Φ h degrees of freedom.
(ii) Matrix C is skew symmetric and consequently iC is Hermitian. They are both singular. (iii) Matrix K is symmetric and negative definite. This result is straightforward to show. First, the contributions associated with the plate deflection satisfy (see also Ciarlet, 2002) Given the boundary conditions on the grounding line Γ g this quantity is itself a norm over V h η . These terms appear now under the minus sign in the definition of the stiffness matrix. Furthermore, matrix K ΦΦ is produced from the discrete bilinear functional Assuming that B is essentially bounded and since meas(Γ f ) > 0, the Poincare-Friedrichs inequality implies that there exists a positive constant C f such that (4.5) The Quadratic Eigenvalue Problem (QEP) above has real eigenvalues (Tisseu and Meerbergen, 2001) with some of them being infinite due to the singular mass matrix. The eigensolutions can be computed, setting λ h = 1/ω h , from the reverse characteristic polynomial (see Lancaster, 1991), (4.6) using direct QEP solvers or linearisation and Generalised Eigenvalue Problem (GEP) solvers; see Tisseu and Meerbergen (2001).

Error estimates and the added mass method
Due to the nonconformity of the proposed finite element scheme, convergence characteristics are bound to depend on two sources of error, namely the approximation error of the interpolation functions adopted and the consistency error due to the inexact nature of the sesquilinear functionals used in the discrete form. A priori error estimates for the eigenfrequencies and eigenfunctions can be derived using an equivalent added mass formulation and standard error estimates for nonconforming thin plate elements (Ciarlet, 2002;Rannacher, 1979). In particular, the following result can be shown Theorem 1. Assuming sufficient regularity of the eigenfunctions, i.e. η k ∈ H 5 (Ω) and Φ k ∈ H 2 (Ω) and h > 0 sufficiently small, there exist constants C ω , C η > 0 such that (4.7) Proof.
Since meas(Γ f ) > 0, there exists a bijection L: (4.8) Then iω , and the following Rayleigh quotient can be defined for a plate with mass distribution MI + L −1 where I denotes the identity operator in V η . This quotient corresponds to an Added Mass Method for the resonating ice shelf problem. For the fully discrete problem it is iω h the Rayleigh quotient reads where, without any loss of generality, it is   η k,h   = ∥η k ∥ = 1. The denominator of the quotient is a weighted inner product and there exists a positive constant C > 0 such that with the same result being true for the fully discrete case as well. Then, the continuity and strong ellipticity of a, a h imply that the continuous problem has an infinite sequence of real, isolated eigenfrequencies and the spectrum of the discrete problem has a total of dim V h η isolated real eigenvalues. The formulation now falls into the category of nonconforming Finite Elements for the thin plate problem analysed by Rannacher (1979). Regarding the principal error terms, the consistency error E a,h of a h is O(h) (Specht, 1988) and the same is true for the error E L,h in the approximation of L −1 by L −1 h using linear triangles. Then, following Rannacher, it is ) .
with Π h the interpolation operator from V η to V h η and c, C , C 1 , C 2 , C 3 positive constants. Inequalities (4.7) now follow if the polynomial degree in the approximation space is considered along with the basic identity ω 2 Other nonconforming element choices for the ice shelf bending problem are also possible. Of particular interest is the Morley element (Morley, 1971) that features quadratic interpolation and only six degrees of freedom per element. Implementation and analysis of the Morley element for ice shelf vibration problems will be the subject of a future study.
Remark (Added Mass Method). Since K ΦΦ is invertible, the last equation of system (4.2) can be written as ) . (4.12) Substituting into the first two equations ] . (4.13) The added mass matrix M a = of the fluid structure interaction system is symmetric and positive definite. To show that it suffices to show that K −1 ΦΦ is negative definite, since for any vector with rows as many as the rows in M a and one column it is (4.14) But it has already been shown that K ΦΦ is negative definite and therefore so is its inverse. This solution approach corresponds to the added mass method and although it might be convenient for small scale problems, the inversion of the sparse matrix K ΦΦ is typically a formidable approach when the dimension of the approximation spaces is large.

Approximation for small eigenfrequencies
For the fundamental and lower order modes, and in particular when |ω| ≪ (1/M) 1/2 , the governing equation (2.10), that contains the term (1 − ω 2 M)η, can be approximated as (4.15) Since typically M = d/L ≪ 1 for large ice shelves, the ice shelf inertia term is neglected in this approximation and Eq. (4.15) is expected to yield very accurate results for relatively small values of the eigenfrequency. In this case, the finite element discretisation using SHEEL produces the Generalised Eigenvalue Problem (GEP) (iω h C + K) U = 0, (4.16) and significant computational savings can be expected. One key issue in this case is that matrix C is singular and skew symmetric. The number of zero eigenvalues for this matrix is even if the dimension is even and odd if the dimension is odd. The following result relates the number of the zero eigenvalues to the mesh attributes.
Let N T denote the total number of nodes in T h with N g the number of nodes in Γ g and N f the number of nodes in Γ f .

Then
Theorem 2. Assume that N g + N f = 2n, n ∈ N. Then the dimension of C is even.
Proof. The total number of degrees of freedom for the problem, is equal to the dimension of C. Each node has a total of 4 degrees of freedom. Three of them are associated with the plate deflection and one with the velocity potential. Nodes on Γ g have therefore only one unknown, while nodes on Γ f have three unknowns. The total number of unknowns is (4.17) This result describes the proper selection of meshes, so that the extra zero eigenvalue is eliminated. The solution of GEP (4.16) can be realised through the use of the reverse polynomial (Lancaster, 1991), similarly to the case of the QEP (4.6). □ In the following section, the effectiveness of the SHEEL element will be verified against a semi-analytical solution and the convergence rates of Theorem 1 will be validated using specific examples.

Numerical results
The SHEEL element has been implemented in a MATLAB R ⃝ code, compatible with the meshing capabilities offered by MATLAB for quasiuniform meshes of 2D areas. This specific code is verified in the following example and in Section 5.2. Subsequently, it is used for the analysis of variable bathymetry and ice shelf thickness cases, as well as a model of the Larsen C ice shelf.

Verification against a semi analytical-solution
In this section the developed hydroelastic element will be verified against a semi-analytical solution for a model problem involving a circular ice-shelf. Although the specific geometry is not realistic, the selected configuration enables the study of convergence characteristics and computational efficiency.
A circular ice-shelf with diameter L = 100 km is considered. The ice-shelf thickness is τ = 300 m and the depth b = 500 m. The ice and water density are set to ρ i = 900 kg/m 3 and ρ = 1000 kg/m 3 respectively. The mechanical properties of ice used are E = 11 GPa and v = 0.3 and the acceleration of gravity is g = 9.81 m/s 2 . Half of the ice-shelf perimeter is assumed to represent the grounding line, while the remaining part is exposed to the ocean. The configuration geometry (in the nondimensional setting) is depicted in Fig. 4. The definition of Γ g and Γ f is A hierarchy of meshes has been applied and the numerical method has been found to converge with increasing number of elements. A mesh of 4064 elements is shown in Fig. 4 along with the resulting sparse pattern of the global mass, stiffness and hydroelastic coupling matrix.
The first 6 eigenmodes for the deflection η are shown in Fig. 5 and the corresponding modes for the velocity potential Φ along with quiver plots for the velocity field in Fig. 6. For η, the real part is shown normalised using the value max Ω |Re(η)|. For the velocity potential the imaginary part is shown and the values are normalised using the value max Ω |Im(Φ)|. The abrupt change of boundary conditions that occurs at the Γ g − Γ f interface creates a weak singularity in the velocity field that manifests as a sudden change in the velocity potential slope along the perimeter of the circle and in particular at locations (x, y) = (0, −0.5) and (x, y) = (0, 0.5). This effect can be clearly identified in the 1st, 2nd and 4th eigenmodes of Φ. In all meshes employed, a node has been allocated to these two locations. The effect of the weak singularity is more clearly demonstrated in Fig. 7, where the measure of the velocity field |∇Φ| is plotted for the first six modes. Although standard finite element analysis implies that in such cases the convergence rate might be compromised and governed by the order of the singularity, still the error reduces with increasing number of elements and the results are indicative of the method's capability to simulate complicated geometries (e.g. re-entrant corners) or complicated boundary conditions that could lead to solutions of reduced smoothness. This attribute is very significant for modelling realistic geometries of ice shelves due to their irregular perimeter shapes.
For the benchmark example selected, a semi-analytical solution in terms of a Fourier-Bessel series expansion is derived.
Setting S = K ∆η, the governing equations (2.10)-(2.11) become  (5.6) The above system of second order Partial Differential Equations can be diagonalised to produce three second order differential equations of the form ∆F n + λ n F n = 0, n = 0, 1, 2. where only the linearly independent solutions that produce bounded fields in the axes origin have been retained in the series expansion.
The solution (5.8) is transformed back to the original fields as follows where a mn , b mn denote the unknown modal amplitudes. This series solution can be truncated and employed, in conjunction with the boundary conditions of the problem and Galerkin's method, to form an algebraic system for the unknown modal amplitudes a nm , b nm . Details of this methodology for the problem considered are provided in Appendix B.
In Fig. 8, the convergence of the SHEEL element as the number of elements N in the employed mesh increases is compared to the convergence of the semi-analytical solution with increasing number of terms N m in the Fourier-Bessel series expansion. In the horizontal axis, for reasons of calibration, the quantity 0.1N 1/2 is selected to indicate the increase in the number of finite elements. Due to the symmetries in the geometry and boundary conditions considered, the eigenmodes appear as symmetric and anti-symmetric. The left hand side graph in Fig. 8 shows the convergence of eigenfrequencies 1, 3, 4 and 6 that correspond to symmetric deformation modes. The right hand side graph shows the convergence of eigenfrequencies 2 and 5 that correspond to anti-symmetric modes. The exact value is considered to be the SHEEL solution with a mesh of 260 096 elements. Both the SHEEL solution and the semi-analytical are found to converge to the same value. The FEM solution though has the advantage that can predict the weak singularity whereas the semi-analytical solution cannot, unless other approaches are considered (e.g. analytical treatment of the singularity in the modal expansions). Finally, in Fig. 9 (left hand side plot), the approximation of eigenfrequencies using SHEEL element and the GEP (4.16) for small ω values is compared to the solution of the full QEP. Up to the first 500 modes, both solutions are in good agreement and the relative error is in all cases less that 3%. As expected the error increases for increasing number of mode. In Fig. 9, the eigenfrequency solution is also compared to the eigenfrequencies of the basin without the ice shelf (blue line) and to that of an ice shelf in vacuo that is represented by the eigenfrequencies of a plate (red line). Similarly to the findings in Papathanasiou et al. (2019), the basin approximation appears to be better for the fundamental and lower order modes, while the hydroelastic solution cannot be approximated adequately as the number of mode increases after a specific threshold.

Convergence characteristics
In order to verify the convergence rate of the SHEEL element, a configuration that allows for solutions with the desired smoothness is selected. In particular, a square ice shelf with edge length of L = 150 km is considered. The ice-shelf thickness is τ = 300 m and the depth b = 500 m. As was the case for the circular ice shelf, the ice and water density are set to ρ i = 900 kg/m 3 and ρ = 1000 kg/m 3 respectively. Young's modulus and Poisson's ratio are E = 11 GPa and v = 0.3 and the acceleration of gravity is g = 9.81 m/s 2 .
Two different geometries for the grounding line are considered. In the first configuration, one side of the square is fixed and hence defines the grounding line. The remaining three sides of the square are exposed to the action of ocean      waves. In the second case, three sides of the square ice shelf are fixed and their fourth side defines the ice shelf ocean interface boundary. The configuration geometry (in the nondimensional setting) is depicted in Fig. 10. The convergence of the 1st, 3rd and 6th eigenfrequency is shown in Fig. 11 for both configurations in double logarithmic scale. The observed rates vary between 1.85 and 2.15 in all cases and are in very good agreement with the theoretically predicted value 2. It is also evident from the results in Fig. 11 that the approximation deteriorates as the mode number increases.

Effects of variable bathymetry and thickness
The SHEEL formulation can be used for cases of variable bathymetry in the sub ice shelf cavity basin and for cases of mild thickness variations as well. In this section a numerical example that demonstrates these capabilities will be presented. A square ice shelf of length L = 150 km is considered. The grounding line extends on three sides, while the remaining one is exposed to the ocean (Fig. 12). The ice and water density are ρ i = 900 kg/m 3 and ρ = 1000 kg/m 3 respectively. Young's modulus and Poisson's ratio are E = 11 GPa and v = 0.3 and the acceleration of gravity is g = 9.81 m/s 2 .
Finally, in Fig. 14, the first 6 eigenmodes for the ice shelf deflection are plotted for cases (a), (b), (c). The results in the first column correspond to profile (a), the results of the second column to (b) and those of the third column to profile (c).
Three different bathymetry and thickness profiles are considered, designated as cases (a), (b) and (c). For case (a) both the seabed profile and thickness have no variation. For case (b) the thickness is constant while the depth varies with respect to the x coordinate. Finally, in case (c) both the bathymetry and thickness profile vary with respect to x. The specific profiles as a function of the nondimensional coordinate x (−0.5 ≤ x ≤ 0.5) are: For the above given data, a mesh with 82 944 elements was selected after a convergence study with coarser meshes. The first 30 nondimensional hydroelastic eigenfrequencies are calculated and depicted in Fig. 13. As expected, the eigenfrequencies increase with increasing basin depth. The eigenfrequencies are found to increase further with the thickness reduction but this increase is very mild. This is the result of two different mechanisms. First, the second moment of area of the cross-section decreases faster than the area itself and this would produce a decrease in the eigenfrequencies due to a decreasing stiffness-mass ratio. However, at the same time, the depth of the basin increases further due to the thickness reduction. This results to an increase of the eigenfrequencies. For the fundamental and lower order modes the effect of the basin depth is found to be dominant. It is worth noting, using the results of this example, that the shift in the spectrum can be significant for variations of the seabed profile and the thickness that could be characterised as mild.
It is evident that the mode shapes are similar for the three bathymetry/thickness profiles with some small discrepancies that manifest as the mode number increases. In particular, differences in the deflection of the 6th mode, near the ocean exposed side, are observed between cases (a), (b) and (c). These changes appear to be significant, to the extent that there is a relocation of the nodal lines, i.e. the zero deflection contours along the ice shelf span. Differences in the modal shapes are expected to manifest more clearly as the number of mode increases.

A model for the Larsen C ice shelf
In this last section, an approximate model for the Larsen C ice shelf will be considered and the SHEEL element will be applied for the calculation of normal modes and eigenperiods. Larsen C is part of the Larsen ice shelf and the fourth largest ice shelf in Antarctica, with area of approximately 44 200 km 2 . It is worth mentioning that in July 2017 a large part, approximately 5800 km 2 , broke from the main body of the Larsen C ice shelf, forming an iceberg (designated A − 68) more than 200 m thick (Hogg and Gudmundsson, 2017). In the following analysis, the fundamental and lower harmonics of Larsen C, using a geometry that represents the ice shelf before the calving event of July 2017, will be studied. An image of the Larsen C ice shelf (before the calving event of July 2017) is presented in Fig. 15. The horizontal geometry of the ice shelf is approximated using the half-disc area with radius L = 200 km defined by the magenta line in Fig. 15 (right hand side). The ice shelf cavity depth and ice shelf thickness are chosen to approximately represent mean values for Larsen C, based on data presented and discussed in Griggs and Bamber (2009) and are set to b = 500 m and τ = 300 m respectively. The ice shelf density is ρ i = 917 kg/m 3 , the water density ρ w = 1027 kg/m 3 and the value for Young's modulus E = 11 GPa.
The problem is nondimensionalised in terms of the spatial domain using the diameter of the ice shelf 2L = 400 km.
Regarding the grounding line boundary and the boundary exposed to the ocean, it is A quasi uniform mesh with approximately 100000 elements is used. Convergence of the solution has been established using a hierarchy of meshes until the solution did not change at four decimal places. Contour plots of the upper surface elevation η for the first eight modes are shown in Fig. 16. The same result for the velocity potential Φ is shown in Fig. 17. Fig. 18 depicts a quiver plot for the velocity field ∇Φ.     Finally, the first 20 eigenperiods of the ice shelf, calculated in hours are depicted in Fig. 19. The fundamental period is approximately 3.9 h. This result is compared with the one derived by the simplified 1D hydroelastic model in Papathanasiou et al. (2019). In that case a value of approximately 4.6 h was obtained. The relative difference between these values is 18%. The 1D model performs reasonably well, while it is expected that the 2D case with the half circular area would give a lower period. This is due to the fact that the 1D model assumes a rectangular area of very large extent along the y direction. Consequently the actual half circular area is in comparison reduced, a fact that implies smaller periods. The same observation applies equally to all the modes that are characterised by oscillations along the y direction, like modes 4 and 10 in Fig. 19. The 1D model yields slightly larger periods and the relative error is 13.2% for mode 4 and 8.2% for mode 10. The proposed Finite Element Methodology can of course be applied for the prediction of normal modes if the reduction in the ice shelf span, due to the 2017 calving event, is taken into account as well and if more detailed descriptions for the seabed topography and ice shelf thickness are considered (Griggs and Bamber, 2009;Fretwell et al., 2013). This analysis will be the subject of future studies.
Conclusions and future research A nonconforming hydroelastic triangle has been developed and employed for modal analysis of ice shelves. The governing equations are coupled and expressed in terms of the ice shelf deflection and the velocity potential of the water basin defined by the sub ice shelf cavity. The new element is based on Specht's interpolation for the bending of Kirchhoff-Love plates and linear interpolation for the velocity potential. Details on the formulation and the explicit form of the hydroelastic coupling matrix are included. The discrete variational form produces a Quadratic Eigenvalue Problem, while a Generalised Eigenvalue Problem, involving only the stiffness and hydroelastic coupling matrix, can be used as a small frequency approximation with very good accuracy. The SHEEL (Specht HydroElastic ELement) is verified against a semi-analytical solution and a priori estimates for the convergence rates have been derived with the use of an equivalent added mass formulation. Several numerical examples indicate that the element is robust and straightforward to apply. For very large ice shelves, like the major Antarctic ones, SHEEL can produce reliable results for eigenperiods in the range of approximately five hours to a few minutes. Future research is aimed towards two different directions. First, for larger eigenperiods the effects of Coriolis acceleration should be included in the formulation (Godin and Zabotin, 2016;Papathanasiou et al., 2019) and future extension of Finite Element solvers for the normal modes of ice-shelves will be focusing on this direction. On the opposite side of the spectrum, for eigenperiods smaller than 1 min i.e. in the range of swells, shorter wavelengths are expected. Hence, the shallow water model (Kalyanaraman et al.) or the thin plate assumption could be inaccurate. In such cases, the present formulation should be enhanced in terms of the plate model (e.g. Reissner-Mindlin theory based elements), as well as the water wave model (e.g. Coupled Mode Systems Belibassakis and Athanassoulis, 2005). Given the SHEEL element and the two aforementioned extensions for larger and lower periods respectively, an integrated simulation tool for the hydroelastic response of ice shelves is envisioned. Finally, the incorporation of viscoelastic effects for ice flexure (MacAyeal et al., 2015) is a straightforward extension given the current SHEEL element formulation.

Acknowledgement
T K Papathanasiou gratefully acknowledges support from the Brunel Research Initiative & Enterprise Fund BRIEF 2018/2019-award number 1069.

Appendix A. The local hydroelastic coupling matrix c K
The exact form of the hydroelastic coupling matrix c K for SHEEL is formulated. The elements of this skew symmetric matrix are calculated form the terms in the bilinear form − ∫ K (χ Φ − wη)dxdy inside element K, using the standard integration formula for linear shape functions in a triangle K. The sparse pattern in matrix c K is The nonzero elements of the upper triangular part can be calculated using (3.7)-(3.10), (3.12b) and (A.1) as row 1 where R is the radius of the circular plate and J Finally, using the requirement that the normal bending moment and the active shear force are zero at the free part of the boundary (see, e.g., Timoshenko and Woinowsky-Krieger, 1959), M n = M r = 0 and Q + ∇T · t = Q − r −1 ∂M rt ∂θ = 0, at r = R and − π 2 < θ < π . One term of the coupling is defined on the angular interval containing the free side of the plate (−π/2 < θ < π/2) and the other on the fixed side of the plate ( π/2 < θ < 3π /2), respectively. Subsequently, the system of equations is projected on the Fourier basis {cos(mθ), m = 0, 1, . . . , M} and {sin(mθ), m = 1, . . . , M} exploiting the orthogonality properties of the latter basis. The following algebraic system is The determinant of the above homogeneous system (B.10)-(B.12) is a function of the frequency. Consequently, a direct search method is applied for the determination of its roots in the vicinity of zero, corresponding to the first eigenfrequencies.
Furthermore, an asymptotic approximation, that has been found to provide reasonable predictions for the lowerorder eigenfrequencies is obtained from the solution of the Helmholtz equation with parameter k 2 = ω 2 /B, where B = (b − d)/diam(Ω). Following an approach similar to the one presented above, in this case the general solution can be represented as follows Φ (r, θ) = ∞ ∑ m=0 J m (kr) [â m cos(mθ ) +b m sin (mθ ) ] . (B.20) The Dirichlet boundary condition at the part of the ice shelf interacting with the ocean and the Neumann boundary condition at the grounding line are now written as The solution can be approximated using the Galerkin method and projecting the above set of equations to the Fourier basis. This procedure will lead to a system similar to that in Eq. (B.10).