Hydroelastic analysis of ice shelves under long wave excitation

The transient hydroelastic response of an ice shelf under long wave excitation is analysed by means of the finite element method. The simple model, presented in this work, is used for the simulation of the generated kinematic and stress fields in an ice shelf, when the latter interacts with a tsunami wave. The ice shelf, being of large length compared to its thickness, is modelled as an elastic Euler-Bernoulli beam, constrained at the grounding line. The hydrodynamic field is represented by the linearised shallow water equations. The numerical solution is based on the development of a special hydroelastic finite element for the system of governing of equations. Motivated by the 2011 Sulzberger Ice Shelf (SIS) calving event and its correlation with the Honshu Tsunami, the SIS stable configuration is studied. The extreme values of the bending moment distribution in both space and time are examined. Finally, the location of these extrema is investigated for different values of ice shelf thickness and tsunami wave length.


Introduction
The catastrophic impact of climate change on the Antarctic Peninsula is examined in the works of Scambos et al. (2003) and Skvarca et al. (1999), where attempts to identify the mechanisms of climate-induced, ice shelf disintegration are made.Ice shelf stability is being re-evaluated as wave trains are becoming rougher and elevated temperatures lead to the further weakening of ice formations (Young et al., 2011).In fact, the question of whether ocean wave forcing acts as a collapse triggering mechanism is thoroughly explored in the literature.In particular, gravity wave forcing is depicted as a major cause of rift propagation within an ice shelf, preceding breakup events (Bromirski and Stephen, 2012).Additionally, the effects of infra-gravity waves and intense storm activity are also considered crucial for ice shelf stability (Bromirski et al., 2010).
The present contribution is motivated by the calving event triggered by the Honshu earthquake-generated tsunami, in March 2011.Observational data showed that the Tsunami generated by the aforementioned earthquake in Japan reached the Sulzberger Ice Shelf in Antarctica and caused the formation of two icebergs, the largest being the size of Manhattan island (Brunt et al., 2011).It is evident that the oscillatory flexural bending, induced by wave excitation, is a primary mechanism for ice shelf and ice tongue calving.Ice-tsunami wave interaction is also manifested in the run-up stage, when drifting ice formations are swept by the incoming long wave.The Tohoku Tsunami exhibited the rare feature of transporting large ice masses, causing significant disruptions on the Kuril Islands shoreline as documented in Kaistrenko et al. (2013).
Due to their structural characteristics, namely their negligible bending rigidity and large horizontal dimension, the dynamic response of ice shelves when interacting with the ocean wave field can be effectively modelled as an initial-boundary value problem of hydroelasticity.Hydroelastic analysis is also applied for the study of ice floes sub-T.K. Papathanasiou et al.: Hydroelastic analysis of ice shelves under long wave excitation jected to ocean forcing (see Squire, 2007).Under the above considerations, ice shelves can be modelled as constrained semi-infinite plates floating over a water region with either zero or non-zero draft (see Sergienko, 2010).Related to ice shelf modelling, a recent work of Bhattacharjee and Guedes Soares (2012) focuses on the frequency domain problem of a floating semi-infinite plate in the vicinity of a vertical wall.A variety of plate edge conditions are examined, including a free, a fixed and a pinned condition at the vertical wall interface.Brocklehurst et al. (2010) present an analytical solution to the problem of a clamped semi-infinite, homogeneous, elastic plate over a constant bathymetry region.Tkacheva (2013) employs an eigenmode expansion for the solution of a fixed plate on a vertical wall under regular wave loading.
The majority of studies consider the case of harmonic wave excitation, which enables the calculation of the floating body response in the frequency domain.In this case, a common line of work is the modal expansion technique, where the elastic deformation is deduced by the superposition of distinct modes of motion (Belibassakis and Athanassoulis, 2005).The hydrodynamic forces are treated primarily through the employment of the Green function method or the eigenfunction expansion matching method.A number of studies have focused on transient analysis of elastic floating bodies, allowing for non-harmonic wave forcing and timedependant loads on the body.These attempts incorporate direct time integration schemes, Fourier transforms, modal expansion techniques and other methods (Meylan and Sturova, 2009;Sturova, 2009;Watanabe et al., 1998).For a nonuniform elastic plate floating on shallow waters of variable depth, Papathanasiou et al. ( 2015) developed a higher-order finite element for the time domain solution of the hydroelastic problem composed of a freely floating or semi-fixed body, while the non-linear transient response is examined in Sturova et al. (2010) by means of a spectral-finite difference method.
In the present contribution, the previous work of the authors on higher order FE schemes (i.e.Papathanasiou et al., 2015) will be applied in the hydroelastic analysis of iceshelves under long-wave excitation.In Sect. 2 the physical domain and the governing equations are presented.The variational formulation of the previously defined initial boundary value problem is discussed in Sect.3. In Sect. 4 a case study, with parameters resembling that of the Sulzberger Ice shelf, is analyzed by means of the proposed methodology.The temporal distributions of the maximum and minimum bending moment values, along with their corresponding location along the semi-fixed floating body are given.Finally, a parametric analysis regarding the location of the occurred extreme bending moment values is performed for different ice shelf thickness and initial disturbance wavelength values.

Physical domain geometry and governing equations
The ice shelf is represented by an elastic, heterogeneous, thin plate with a fixed edge, extending infinitely at the y direction (vertical to the page).The plate of horizontal length L, rests on a layer of inviscid, incompressible fluid over an impermeable bottom.Assuming shallow water conditions, the long wave approximation (i.e.wavelength much greater than water depth) can be employed.The last assumption allows for dimensionality reduction, resulting in a 1-D system of equations, since now the z component of the fluid velocity is considered negligible.The domain is divided into regions S 0 ≡ (0, L) and S 1 ≡ (L, ∞), with the hydroelastic coupling taking place at the former (Papathanasiou et al., 2015).In S 0 , the plate deflection coincides with the water upper surface elevation η(x, t).The fluid velocity potential in the two regions is denoted as ϕ 0 and ϕ 1 respectively.In order to account for the draft of the plate, the variable bathymetry function b(x) = H (x) − d(x), where B(x) is the water depth and d(x) = τ (x) ρ i /ρ w the draft of the plate, τ (x) being the plate thickness, is defined.The ice and water density are ρ i and ρ w , respectively.The flexural rigidity of the plate is given by D(x) = E τ 3 /12(1 − v 2 ), with E being the Young's Modulus of ice and v the Poisson's ratio.The mass per unit length of the plate is denoted by m and dropping tildes, the governing system of differential equations is reduced to (see also Sturova, 2009), where a superimposed dot denotes differentiation with respect to time while an index following the a function denotes differentiation with respect to the spatial variable.In addition, the coefficients appearing in the above equations are defined as The bending moment and shear force along the shelf are given by the expressions M b = K η xx and V = (K η xx ) x , respectively.

Stress distribution within the ice beam
In agreement with the Euler-Bernoulli beam model, the normal stress varies linearly along the z direction.The maximum normal stress value at any given cross section is The shear stress distribution, as derived from equilibrium relations, varies quadratically along the vertical direction.Maximum shear stress, located at the neutral fibre is, (5) The above system of equations is supplemented with boundary, interface and initial conditions.At the fixed end, simulating the ice shelf grounding line, the deflection and slope are set to zero.At the free edge of the plate, representing the ice shelf tip facing the ocean, zero bending moment and shear force is imposed.These conditions read The water velocity is assumed zero underneath the grounding line and thus the velocity potential gradient vanishes, At the interface between S 0 and S 1 , assuming energy and mass conservation, the following matching conditions are derived (Stoker, 1957;Sturova, 2009): The ice shelf is assumed to be initially at rest, while an incoming long wave transverses region S 1 and reaches the free edge of the shelf.The initial boundary value problem formulation is thus completed with the following conditions, η(x, 0) = η(x, 0) = ϕ 0x (x, 0) = 0, x ∈ S 0 and (9a) In the last of Eq. (9b), F (x) denotes the free surface elevation caused by the Tsunami wave at an initial time, at an area distant to the ice shelf edge.

Finite elements -variational formulation of the governing equations
In order to derive the variational formulation of the above problem, Eqs. ( 1)-( 3) are multiplied by the weight functions ν, −w 0 and w 1 , respectively.Integration by parts yields Using the conditions described in Eqs. ( 6)-( 8), and adding Eqs. ( 10)-( 12) the equivalent semi-discrete variational problem is formulated as Papathanasiou et al. (2015).Find η, ϕ 0 and ϕ 1 , such that for every ν, w 0 and w 1 at any given moment in time it holds where and while superscript h denotes spatially discrete quantities.
A special hydroelastic element is developed and employed for the solution of Eq. ( 13).The reader is directed to the previous work by Papathanasiou et al. (2015) for more details concerning the proposed finite element scheme.The interpolation degree selected features 5th order Hermite polynomials for the beam deflection/upper surface elevation in the hydroelastic region and 4th order Lagrange polynomials for the approximation of the water velocity potential.Hence, within each element, η(xt) and ϕ 0 (x, t) are approximated by where the positive constant R 1 is selected large enough so that any disturbance propagating inside S 1 does not reach point R within the time interval of interest.Fourth order Lagrange shape functions are used for the interpolation of the velocity potential ϕ 1 .By substituting the above approximate solutions into Eq.( 13) and letting the weight functions ν assume the form of the Hermite C 1 shape functions while w 0 and w 1 are substituted by the Lagrange C 0 shape functions, the resulting system is derived in the form, M ü + C u + K u = 0, with the vector of unknowns being u = [η η x ϕ 0 ϕ 1 ] T .
In the present section, the simplistic, mechanical model described above will be employed for the calculation of the hydroelastic response of the Sulzberger Ice shelf under long wave forcing.The SIS is simulated by a semi-fixed plate of 100 km in length.For the employed bathymetric profile, mean depth values were used.In Brunt et al. (2011) it is mentioned that the water column depth in front of the ice-shelf is 150 m, while it increases to 800 m within 100 km from the ice shelf front.Thus, the ocean depth under the ice shelf is assumed to be 150 m, while a mildly sloping bottom is considered over a distance of 100 km from the edge of the ice shelf (see Fig. 1).The water depth increases from 150 m, at x = 1 to 800 m at x = 2 (Brunt et al., 2011).The initial, bellshaped free-surface elevation considered in the following examples is where A represents the amplitude, x 0 the origin location, w the half-wavelength and µ a smoothness parameter controlling the steepness of the initial pulse.Finally, the material constants selected are as follows: ice shelf density ρ i = 922.5 kg m −3 , water density ρ w = 1025 kg m −3 , Young's modulus E = 5 × 10 9 Pa and Poisson's ratio v = 0.3 (see also Sturova, 2009).The acceleration of gravity is g = 10 m s −2 .In the following analysis 400 hydroelastic elements have been used, while the number of elements in region S 1 is selected such that the element size is the same for both regions S 0 and S 1 .Numerical experiments have shown that this discretization ensures convergence, as any further refinement induces virtually no variation of the results.Finally, the Newmark method has been employed for time integration.The non-dimensional time interval T = 70 (corresponding to 7000 s) is considered.At first, the effects of an initial pulse with A = 0.5 m, µ = 50 × 10 5 m −2 and w = 8000 m are considered.In Fig. 1, a visual representation of the given pulse propagation is plotted.The bell-shaped disturbance is split into two waves travelling in opposite directions.The pulse propagating to the left (towards the negative x axis) is partially reflected when reaching the bathymetric variation at x = 2.As the wave propagates over decreasing depth, its amplitude increases, while the velocity decreases.The velocity reduction is evident in the curved trajectory path for 1 ≤ x ≤ 2, as shown in Fig. 1.At x = 1, the wave impacts the ice shelf free edge, initiating the propagation of the hydroelastic wave while it is partially reflected.The hydroelastic wave, featuring dispersive characteristics, is fully reflected at the grounding line (fixed end, x = 0), at t ≈ 67.The dispersive nature of the hydroelastic pulse can also be seen in Figs. 2 and 3, manifested as the formation of smaller disturbances preceding the main elevation wave.These disturbances reach the grounding line at earlier times than the main pulse and lead to an increase of the bending moment locally.This phenomenon is displayed in the maximum and minimum bending moment time profiles (Figs. 2 and 3) as spikes, located at x = 0 and appearing in the time interval from t ≈ 55 to t ≈ 65.
The present analysis aims to provide some first and simple means for the estimation of long wave impact on floating, slender formations and their response, as a first step towards the hydroelastic modelling of ice shelves.As illustrated in Fig. 1, phenomena, such as wave reflection, hydroelastic dispersion, bending moment variation are well reproduced.In order to investigate the generated stress field within the floating body, the bending moment distribution is examined.Bending moment distributions are directly linked to maximum normal stress values.In particular, for notched or pre-cracked specimens it is usually those normal stresses that mostly influence crack initiation and propagation.The latter phenomena are crucial when a pre-existing crack happens to be inside a tensile zone of large magnitude.
Typically, for the bending of thin beams, the normal stresses due to bending are dominant and as a first approach, shear stresses may be neglected.In Figs. 2 and 3 the maximum and minimum bending moment temporal and spatial distributions are shown.For the maximum bending moment, the temporal distribution is shown in a thick red line, while the location of the corresponding values along the ice shelf is given by the thin black line (Fig. 2).When at rest, the maximum bending moment is zero in absence of flexural effects.Immediately after impact, t ≈ 30 the maximum bending moment is seen to increase.The location of the maximum bending moment value is found to follow the main pulse towards the fixed end.At t = 34, the entire pulse has passed underneath the floating cantilever, causing an increase in maximum bending moment.At the same time, the location of the maximum value for the bending moment is shifted back near the free edge.This is due to the fact that the entire wavelength of the initial pulse has passed underneath the floating cantilever, causing the tip to bend again as it recovers  to the initial undeformed state.At t = 44 the location of the maximum bending moment value is shifted towards the ice shelf tip once again.The above can be attributed to flexural effects taking place at the right side of the propagating disturbance.As the hydroelastic wave propagates away from the free edge, the tip is restored to its original position causing additional flexing in the interior of the cantilever.Due to the fact that, in the present work, the grounding line is simplistically modelled as a fixed boundary, the global bending moment extrema are found at the fixed edge, at the time of reflection t = 67.Prior to full reflection, a series of spikes in the maximum bending moment distribution are caused by the dispersed hydroelastic waves reaching the fixed edge before the main pulse.
As shown in Fig. 3, the minimum bending moment intensifies until the entire pulse wavelength has passed under the floating cantilever, at which point the minimum bending moment value remains virtually constant up to the arrival of the dispersed wave train at the fixed edge.
Notably, the notion that the pulse will reach the fixed end is rather unrealistic.The induced flexural effects will cause the bending failure of the semi-fixed floating body long before the hydroelastic pulse arrives at the grounding line.As seen in Figs. 2 and 3, the maximum and minimum bending moment values reach a plateau approximately after the full disturbance passes underneath the ice shelf.Considering the effects before the hydroelastic wave train reaches the grounding line, namely for short times after the long wave impact, the corresponding location of the given extreme bending moment value along the ice shelf may be linked to both ice shelf thickness and initial disturbance form.Figure 4 displays a parametric study of the extreme bending moment value location for different ice shelf thickness and tsunami wavelength values.Variable ξ denotes the distance from the free edge up to the location of the extreme value along the semifixed floating body x = 1.In all cases, the extreme bending moment values have been considered in a time interval excluding the effects of the fixed end forcing (t ≥ 50).In that manner, Fig. 4 demonstrates the location of extreme bending moments for the phase during which the main pulse enters the region of hydroelastic interaction.As can be seen in the aforementioned figure, the location of the extreme bending moment is relatively insensitive to variations of the wavelength.For thickness values of 80 and 100 m, this location is found to be at about 2 % of the ice shelf length (2 km into the 100 km long ice shelf), calculated from the free edge.The above results are found in agreement with the work of Squire (1993), where the breakup of shore fast ice, modelled as a semi-infinite, thin floating plate, is investigated in the frequency domain.Furthermore, as thickness increases, the location of extreme values seems to shift towards the interior of the ice shelf.For a thickness of 120 m, location ξ is placed at approximately 10 % of the ice shelf length and features a slight variation with increasing initial disturbance wavelength.However, this variation is very small when compared to the total length of the beam.Another interesting feature is that in this last case (120 m thickness) the maximum absolute value found corresponds to negative values of the bending moment (see Fig. 5), whereas for thickness val- ues of 80 and 100 m the maximum absolute bending moment values are found to be positive (sagging moments).This feature explains the different shape of the 120 m curve in Fig. 4, when compared to the curves corresponding to 80 and 100 m, which closely resemble one another.The fact that in the case of 120 m the extreme bending moment values are negative might be attributed to the beam thickness being very large compared to the water depth under the ice shelf.Finally, these results are strongly dependent on the form of the incoming wave, in the sense that if another wave profile instead of an elevation pulse is chosen, the bending moment fields will be of a different nature.

Conclusions
In the present work, the transient hydroelastic response of a semi-fixed floating cantilever, resembling an ice shelf, is studied by means of a higher order finite element scheme.The simple model derived above is able to provide valuable information regarding the kinematic and stress fields induced by long wave forcing on an ice shelf.An illustrative case study is presented with parameters selected so as to approximately simulate the Sulzberger Ice Shelf topology and the relevant calving event conditions, initiated by the 2011 Honshu Tsunami.Bending moment profiles, as generated by a long wavelength elevation pulse, are studied and critical points of the induced stress field are located.During the wave entry in the hydroelasticity dominated region, the locations of extreme bending moments is found to be relatively insensitive to the excitation wavelength for given ice shelf thickness values.Important extensions of the present study include 3-D hydroelastic interaction, as well as the investigation nonlinearity effects of both in the hydrodynamic model and in the elastic subregion.Finally, the study of tsunami-ice interaction in the run-up stage constitutes another possible future research direction in the context of the present model applications.

Figure 1 .
Figure 1.Space-Time plot of the bell-shaped pulse propagation.The bathymetric profile is shown in a schematic below.All dimensions are normalized with respect to the plate length L = 100 km.

Figure 2 .
Figure 2. Maximum bending moment temporal profile (red thick line) and location of corresponding values along the floating cantilever (black line).A detailed figure of the profile after the wave impact is presented, along with representative snapshots of the deformed ice shelf.

Figure 3 .
Figure 3. Minimum bending moment temporal profile (blue thick line) and location of corresponding values along the floating cantilever (black line).A detailed figure of the profile after the wave impact is presented, along with representative snapshots of the deformed ice shelf.

Figure 4 .
Figure 4. Location of extreme bending moment along the semifixed floating body for various values of thickness and initial disturbance wavelength.Variable ξ measures the distance of the point of occurrence of the extreme value from the free edge.

Figure 5 .
Figure 5. Plot of maximum and minimum bending moment value distributions for w = 8000 m and τ/L = 0.0012.Extreme bending moment value is negative, during the entry phase, for an ice shelf thickness value of 120 m.