SWIIFT v0.10: a numerical model of wave-induced sea ice breakup with an energy criterion

The wave-induced breakup of sea ice contributes to the formation of the marginal ice zone in the polar oceans. Understanding how waves fragment the ice cover into individual ice floes is thus instrumental for accurate numerical simulations of the sea ice extent and its evolution, both for operational and climate research purposes. Yet, there is currently no consensus on the appropriate fracturing criterion, which should constitute the starting point of a physically sound wave–ice model. While fracture by waves is commonly treated within a hydroelastic framework and parametrised with a maximum strain-based criterion, in this study we explore a different, energy-based, approach to fracturing. We introduce SWIIFT (Surface Wave Impact on sea Ice – Fracture Toolkit), a one-dimensional model based on linear plate theory, that can produce time-domain simulations of wave-induced fracture, into which we incorporate this energy fracture criterion. We demonstrate SWIIFT with simple simulations that reproduce existing laboratory experiments of the fracture by waves of an analogue material, allowing qualitative comparisons and validations of the energy fracture criterion. We find that under some wave conditions, identified by a dimensionless wavenumber, corresponding to in situ or laboratory wave-induced fracture, the model does not predict fracture at constant curvature; thereby calling into question the appropriateness of parametrising sea ice fracture with a maximum strain criterion.

1 Introduction

In the Arctic, the newly available open ocean areas (Raphael et al., 2025) have exposed sea ice to the effects of stronger and more frequent wave events (Thomson and Rogers, 2014). The remaining sea ice is also overall younger, thinner, more fragile and therefore more likely to be fragmented by winds, ocean currents and waves (Stroeve et al., 2012; Stroeve and Notz, 2018). In unconsolidated or fragmented ice, waves are less attenuated (Collins et al., 2015; Ardhuin et al., 2020) and can therefore propagate further into the consolidated part of the ice cover – the ice pack – and break it to a greater extent, thus enabling a wave–ice positive feedback loop (Thomson, 2022; Horvat, 2022).

This wave-induced breakup results in an assembly of floes with sizes ranging from a few metres to hundreds of metres, defining what is generally referred to as the marginal ice zone (MIZ; see Dumont, 2022, and many others), a region whose dynamics is affected by wave propagation. Fragmented sea ice behaves very differently from the consolidated ice pack. It is more mobile, potentially reaching a free drift state, with ice internal stress no longer resisting motions imparted by winds or currents (Alberello et al., 2020), tides (Watkins et al., 2023), or waves (Auclair et al., 2022; Womack et al., 2022), even at high ice concentration. It is also more sensitive to melt (Horvat et al., 2016; Thomson, 2022), as the ratio of lateral surface (proportional to the perimeter and exposed to the ocean) to top surface (exposed to air) increases when the horizontal extent of a floe diminishes, eventually accelerating the disintegration of smaller floes (Toyota et al., 2025). As a result, the MIZ response to storms can result in quick and large sea ice losses (Smith et al., 2018; Blanchard-Wrigglesworth et al., 2022; Cavallo et al., 2025), which could amplify the observed sea ice decline (Asplin et al., 2012; Thomson and Rogers, 2014). Concomitantly, high-frequency sea ice extent variability is missing in state-of-the-art climate models (Blanchard-Wrigglesworth et al., 2021), in which sea ice fragmentation by waves is not accounted for. This suggests an improved representation of the MIZ in sea ice models is essential to deliver accurate predictions of the sea ice evolution, over both short-term and climate time scales. However, it remains challenging as the physical processes controlling sea ice breakup are still largely unascertained, or rest on hypotheses that are not fully backed up by observations.

The first step in this model development should be the identification of a fracture or disintegration criterion, allowing to determine under which wave forcing (amplitude, wavelength, spectral distribution) and ice conditions (thickness, mechanical stiffness) the ice will fragment into floes. To our knowledge, there is actually neither complete physical evidence nor clear consensus within the sea ice community on this criterion. To our knowledge again, all the current wave breakup modelling approaches are based on local flexural stress or strain reaching a prescribed critical value, or threshold. Behind this viewpoint is the consideration that maximum deformation will either occur at the crests and troughs of waves (for example,  Dumont et al., 2011) or at the wave front (Tkacheva, 2001). Voermans et al. (2020) extended this formalism by combining a strain threshold with wave characteristics into a dimensionless quantity, the value of which separates breakup from non-breakup. The universality of this approach was later called into question (Passerotti et al., 2022). When modelling individual floes, any local threshold is however susceptible to be exceeded over large spans of the floes, which makes super-parametrisations necessary. The location of maximum strain or stress is often considered for the fracture location (Dumont et al., 2011; Williams et al., 2017; Montiel and Squire, 2017; Mokus and Montiel, 2022), but other methods, such as computing the strain between successive wave crests and troughs only have been used (Horvat and Tziperman, 2015).

These local threshold-based criteria are consistent with (and usually come hand in hand with) an hydroelastic representation of the wave–ice system, on which a large fraction of the modelling research on wave–ice interaction lies (Squire, 2020). Wave-induced sea ice fracture has thus naturally been considered through this lens (for example,  Fox and Squire, 1991; Montiel and Squire, 2017; Zhang and Zhao, 2021; Mokus and Montiel, 2022); even though more novel and computationally involved approaches exist (Herman, 2017; Ren et al., 2021; He et al., 2022). In the hydroelastic framework, the ice is assimilated to an elastic plate that is thin enough for the variations in the buoyancy forces acting on it to be negligible, and that therefore conforms exactly to the shape of the ice–ocean interface. When associated with a critical strain fracturing criterion, this framework has shown agreement with observations (Kohout et al., 2016; Voermans et al., 2020). It has also allowed wave and floe-resolving numerical simulations to generate steady-state floe size distributions (Kohout and Meylan, 2008; Horvat and Tziperman, 2015; Mokus and Montiel, 2022), and has therefore percolated into coupled global sea ice models (Roach et al., 2019; Bateson et al., 2020; Yang et al., 2024).

The current contribution digs into the question of the criterion for the fracturing of consolidated sea ice by waves, that is, flexural brittle failure. In this, we are motivated by recent laboratory results investigating the response of an ice analogue material to wave forcing (Auvity et al., 2025). In particular, these authors highlighted that the curvature at which their material broke is not constant, but depends monotonically on the applied wavelength. The spread in reported sea ice critical strains (Kohout and Meylan, 2008; Voermans et al., 2020) could thus be an artefact hiding such a relationship. In this context, we investigate an approach based on a model of fracture propagation in elastic solids (Griffith, 1921) which is common in the field of fracture mechanics. It opposes the energetic cost of creating new surfaces to the elastic energy stored in a material. The resulting energy-based fracture criterion includes the effect of bending deformation as it depends on the associated elastic deformation energy, but is non-local as it is integrated over the length of the deformed ice floe. It leads to a unique solution. Since the original work of Griffith (1921), this model has been updated (Francfort and Marigo, 1998; Francfort, 2021) and built upon specifically for application to sea ice (Mulmule and Dempsey, 1997; Balasoiu, 2020; Ren et al., 2021). Measurements of sea ice fracture toughness, which can be linked to the energetic cost of fracturing, have been compiled (Dempsey, 1991; Schulson and Duval, 2009) and an extensive body of work also exists on freshwater ice (for example,  Gharamti et al., 2021a, b).

With the intent on focusing on the wave-induced deformation and resulting fracturing of brittle ice, we have implemented this energy-based criterion in a framework that differs from the hydroelastic representation in that the ice is not assumed to conform to the water surface, but freely deforms within the wave field as a result of the local buoyancy and gravity forces. We neglect other processes affecting the seasonal ice zone (SIZ; see Roach et al., 2025), such as thermodynamics; in particular, we do not handle ice formation within a wave field, and we restrict our study to the case of brittle fracture, excluding the disintegration of a more granular material (dislocation or melting of forming ice). The resulting simple, yet versatile, 1D model also accommodates a strain-based fracture criterion. It thereby allows investigating the effect of using either criterion on the occurrence of the fracture and, eventually, on the extent of a simulated MIZ, and shape of the associated floe size distribution. Importantly, our model can be stepped forward in time, so that we can use it to follow the propagation of a fracture front as a function of time; in contrast to being able to solely recover the final, fractured state (Horvat and Tziperman, 2015; Mokus and Montiel, 2022). We present an illustration of this capability in Sect. 2.6. In the present paper, we exploit this model in another use case, the comparison to laboratory experiments on fracture, conducted on an analogue material to sea ice. We pursue this comparison with the particular aim of validating-invalidating the applicability of the energy-based fracturing criterion.

The paper is divided as follows: in Sect. 2, we give a general mathematical description of the model, including the treatment of sea ice as an elastic plate, the formulation of the breakup criteria, and the representation of waves. In Sect. 3, we give specific information pertaining to the numerical results we present in Sect. 4, that is, the particular setup of the model in this study. We discuss these results in Sect. 5.

2 Floe and fracture model

In light of the objectives motivating our approach, stated in the introduction, we present in Sect. 2.1 the main physical hypotheses made to achieve a simple, versatile, numerically efficient, yet physically sound model. In Sect. 2.2, we detail our approach to deriving floe deformation, used as an input for the fracture parametrisations presented in Sect. 2.3, and forced by waves discussed in Sect. 2.4. We present numerical aspects in Sect. 2.5 and conclude this Section with an example of time simulation in Sect. 2.6.

2.1 Main hypotheses

A common assumption behind wave–elastic plate interaction models (for example,  Fox and Squire, 1991; Tkacheva, 2001; Mokus and Montiel, 2022) is to consider the plate thin enough for variations of the buoyancy force acting on it to be negligible. The plate is, however, subjected to the fluid pressure acting on its bottom side, and it is assumed that the plate conforms to the fluid motion at all times. Fluid pressure is determined by solving for the fluid flow, typically by assuming a potential flow and harmonic solutions, with the plate exerting a boundary condition on the fluid domain. To develop the model presented herein, we adopt a different approach, motivated by our interest in the ice deformation and fracture, whereas the focus of fluid-centred models has historically been that of wave scattering and attenuation by the plate. The interested reader can found a comparison between the two approaches in Appendix B.

In this study, the ice cover is considered thick enough for the local changes in buoyancy force not to be negligible. We do not explicitly resolve the fluid flow underneath the plate, and impose no condition on the ice–ocean interface. Instead, we solve for the vertical deflection of the ice stemming from the local balance of gravity and buoyancy forces driven by the sea surface displacement. The fluid surface thus acts as a forcing term, which is made aware of the presence of the ice floes only through parametrised attenuation; a consequence is that floes can locally be immersed. While we limit ourselves to the case of linear elasticity and linear wave forcing, our mechanical formulation interpolates between the limits of an elastic floe that conforms perfectly to the wave surface, and of a rigid floe only capable of solid motion, which can therefore be submerged. We quantify this behaviour with the dimensionless wavenumber kL_D, that relates the wavenumber of the forcing wave k (formally introduced in Sect. 2.4.1) to the flexural length of the floe L_D (formally introduced in Sect. 3.1). The small kL_D limit (long wave, compliant floe) corresponds to the strain formulation of Dumont et al. (2011), while the large kL_D limit (short wave, rigid floe) corresponds to their stress formulation.

Ice formation and melt are assumed to happen at timescales beyond that of wave-induced fracture, so that they can be neglected. We do not consider the reflection of waves at the ice–water interface (for example, Mokus and Montiel, 2022), nor viscous deformation of the plate, the compression of an array of floes due to wave radiation (for example, Herman, 2018), or any surge motion, and take note that the model would be more complete if the pressure forcing associated with the contact between the ice and the water was explicitly taken into account.

2.2 Governing mechanical equations

Our one-dimensional model considers a fluid volume of finite or infinite depth, equipped with a Cartesian coordinate system (x,z) where z is the vertical coordinate oriented upward, as shown in Fig. 1. We assume translational invariance in the second horizontal direction. The domain is populated with floating ice floes of prescribed positions and lengths, which may not overlap. Any part of the domain not covered with ice is deemed to be open water.

As in the work of Meylan et al. (2015), we model floes as elastic plates, and we derive the deformation of the ice cover using the Kirchhoff–Love thin-plate theory. The ice is thus considered homogeneous, isotropic, and transversally loaded by body forces. We assume a constant thickness h along a given floe, although different floes can have different thicknesses. The two forces acting on the ice are buoyancy and gravity. In the case of a fluid at rest (no waves), equating the gravity force per unit area and the buoyancy force per unit area (thus applying Archimedes' principle) allows expressing the draught of a floe, d, as

$$\begin{array}{}\text{(1)}& {\mathit{\rho }}_{\mathrm{i}}h\mathit{g}-{\mathit{\rho }}_{\mathrm{w}}d\mathit{g}=\mathrm{0}\iff d={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}}h,\end{array}$$

where ρ_i is the density of the ice, ρ_w the density of the ocean water and g the gravitational field.

We now move away from this rest state, and impose a perturbation of the fluid surface η(x). Because the propagation speed of elastic waves in sea ice is several orders of magnitude greater than that of surface gravity waves (Moreau et al., 2020a), we consider this perturbation and the resulting floe deformation to be quasi-static. Thus, we only consider one independent variable, the space coordinate x, and no explicit time dependency. The vertical displacement w(x) corresponding to this perturbed state is determined by the local balance between gravity and buoyancy. The weight per unit area of the ice is still ρ_ihg. However, the height of displaced fluid now corresponds to the difference between the fluid surface η(x) and the displaced bottom of the ice floe, w(x)−d. Locally, the buoyancy force per unit area is thus $-{\mathit{\rho }}_{\mathrm{i}}\left[\mathit{\eta }\right(x)-(w\left(x\right)-d\left)\right]\mathit{g}$; this is illustrated in Fig. 1. The floe is then subjected to the resulting body force

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{q}\left(x\right)& {\displaystyle }={\mathit{\rho }}_{\mathrm{i}}h\mathit{g}-{\mathit{\rho }}_{\mathrm{w}}\left[\mathit{\eta }\right(x)-(w\left(x\right)-d\left)\right]\mathit{g}\text{(3)}& {\displaystyle }& {\displaystyle }=-{\mathit{\rho }}_{\mathrm{w}}\left[\mathit{\eta }\right(x)-w(x\left)\right]\mathit{g},\end{array}$$

and projecting q onto the vertical axis gives

$$\begin{array}{}\text{(4)}& q\left(x\right)={\mathit{\rho }}_{\mathrm{w}}g\left[\mathit{\eta }\right(x)-w(x\left)\right].\end{array}$$

Figure 1 Schematic of a deformed ice floe in a wave field. The horizontal dashed line represents the sea surface at rest in the absence of ice, which we use as the reference level. The dashed rectangle represents a floe at rest. A perturbation of the fluid surface (solid line) in the free-surface regions imposes a deflection of the ice floe (solid-lined shape). This corresponds to a model output, with ice thickness 50 cm, floe length 120 m, forcing wave with amplitude 20 cm and wavelengths 76.4 m (open water) and 84.3 m (ice-covered water).

Using the bending equation of a loaded plate, we then obtain a differential equation on the deflection of the floe,

$$\begin{array}{}\text{(5)}& D{\displaystyle \frac{{\mathrm{d}}^{\mathrm{4}}w}{\mathrm{d}{x}^{\mathrm{4}}}}={\mathit{\rho }}_{\mathrm{w}}g\left[\mathit{\eta }\right(x)-w(x\left)\right],\end{array}$$

where we take advantage of the simplifying hypotheses made on our geometry. In this equation, we introduced the flexural rigidity $D=\frac{Y{h}^{\mathrm{3}}}{\mathrm{12}(\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}})}$, characterising the ability of the plate to resist bending, with Y and ν the Young's modulus and Poisson's ratio of the plate. We complete Eq. (5) with free-edge boundary conditions, that is, vanishing moment and force at both ends of the ice floe of length L. We choose a reference frame local to the floe, where x=0 corresponds to its left edge, so that the complete boundary problem can be written

$$\left\{\begin{array}{llr}\frac{{\mathrm{d}}^{\mathrm{4}}w}{\mathrm{d}{x}^{\mathrm{4}}}=\left(\frac{D}{{\mathit{\rho }}_{\mathrm{w}}g}{)}^{-\mathrm{1}}\right[\mathit{\eta }\left(x\right)-w\left(x\right)]& x\in [\mathrm{0},L]& \left(\mathrm{6}\mathrm{a}\right)\\ \frac{{\mathrm{d}}^{\mathrm{2}}w}{\mathrm{d}{x}^{\mathrm{2}}}=\mathrm{0}& x\in \mathit{\{}\mathrm{0},L\mathit{\}}& \left(\mathrm{6}\mathrm{b}\right)\\ \frac{{\mathrm{d}}^{\mathrm{3}}w}{\mathrm{d}{x}^{\mathrm{3}}}=\mathrm{0}& x\in \mathit{\{}\mathrm{0},L\mathit{\}}.& \left(\mathrm{6}\mathrm{c}\right)\end{array}\right.$$

For prescribed wave conditions and material properties, solving Eq. (6) thus provides the deflection w of the floe.

We focus here on the bending undergone by an elastic plate because of a perturbation of its fluid foundation. We recall that we do not explicitly resolve the fluid motion itself, nor the translational motions imparted to the plate by the fluid. In particular, we thus neglect surge motion, that is, ice drift in the direction of wave propagation.

2.3 Fracture

2.3.1 Energy criterion

Unlike the prevalent maximum strain formalism commonly used by the sea ice community when modelling wave–ice interactions (for example, Kohout and Meylan, 2008; Dumont et al., 2011; Horvat and Tziperman, 2015; Mokus and Montiel, 2022), we develop a breakup criterion from the framework of fracture mechanics, based on the consideration that in solid, brittle materials, fracture happens to minimise the internal energy associated with deformation (Griffith, 1921). In this framework, the total energy to be minimised is the sum of the elastic energy E_el associated with the deformation (in our case, bending) of the material, and of the fracture energy E_fr associated with the creation of new surfaces around a crack. This energy decomposition is consistent with mode I fracturing, which in the case of our model translates to vertical fractures due to in-plane tensile stress.

The elastic energy density (per unit length in the transverse horizontal direction, and per unit thickness) stored in a material that is elastic, isotropic, and homogeneous, and stretched only in the longitudinal direction, is

$$\begin{array}{}\text{(7)}& W_{\mathrm{el}}={\displaystyle \frac{D}{\mathrm{2}h}}{\mathit{\kappa }}^{\mathrm{2}}\left(x\right)\end{array}$$

with

$$\begin{array}{}\text{(8)}& \mathit{\kappa }\left(x\right)={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}w}{\mathrm{d}{x}^{\mathrm{2}}}}\end{array}$$

the local linear floe curvature due to the deformation. Equation (7) stems from integrating the density of elastic energy (per unit volume) along the axis normal to the neutral plane of the plate (in our case, the z-direction), and takes into account stretching or compression in the directions of the plane (in our case, simply the x-direction). By integrating Eq. (7) along the floe, we obtain the surface energy density,

$$\begin{array}{}\text{(9)}& E_{\mathrm{el}}={\displaystyle \frac{D}{\mathrm{2}h}}\underset{\mathrm{0}}{\overset{L}{\int }}{\mathit{\kappa }}^{\mathrm{2}}\left(x\right)\mathrm{d}x.\end{array}$$

The fracture energy density E_fr relates to the energy required to create a new surface. In the case where the only admissible fractures vertically break the ice through its entire thickness, we simplify the formulation from Francfort and Marigo (1998) as

$$\begin{array}{}\text{(10)}& E_{\mathrm{fr}}=N_{\mathrm{fr}}G\end{array}$$

with N_fr the number of fractures, and G the energy release rate. Again, note that this energy is expressed per unit surface normal to the x-direction.

To determine whether a floe breaks, we compare two energy states: that of the unbroken, deformed floe, and a hypothetical state in which this floe has fractured into several fragments. In the former state, the elastic energy, noted $E_{\mathrm{el}}^{\mathrm{0}}$, is that of the deformed floe, as defined in Eq. (9). In the latter state, the total elastic energy, noted $E_{\mathrm{el}}^{\mathrm{s}}$, is the sum of the elastic energies of the individual newly broken floes. If the broken state is – from an energy standpoint – favourable, it should replace the unbroken state. Formally, we look for the finite set of the fracture locations, ${\mathit{x}}^{\mathrm{fr}}=\mathit{\{}{x}_j\mid x_j\in (\mathrm{0},L\left)\mathit{\right\}}$. This set has size N_fr, the number of fractures, dividing the original floe into N_fr+1 fragments. It should minimise the free energy F defined as

$$\begin{array}{}\text{(11)}& F\left({\mathit{x}}^{\mathrm{fr}}\right)=E_{\mathrm{fr}}\left(N_{\mathrm{fr}}\right)+E_{\mathrm{el}}^{\mathrm{s}}\left({\mathit{x}}^{\mathrm{fr}}\right)-E_{\mathrm{el}}^{\mathrm{0}}\end{array}$$

with the additional constraint that for breakup to occur, we must have F<0. In other words, a floe breaks if the elastic energy released by the breakup exceeds the energetic cost of the breakup. If no such set can be found, we conclude that the current deformation of the floe is not sufficient to fracture it. The post-fracture elastic energy expands to

$$\begin{array}{}\text{(12)}& E_{\mathrm{el}}^{\mathrm{s}}\left({\mathit{x}}^{\mathrm{fr}}\right)={\displaystyle \frac{D}{\mathrm{2}h}}\sum _{j=\mathrm{0}}^{N_{\mathrm{fr}}}\underset{x_j}{\overset{x_{j+\mathrm{1}}}{\int }}{{\mathit{\kappa }}_j}^{\mathrm{2}}\left(x\right)\mathrm{d}x\end{array}$$

with x₀=0 and $x_{N_{\mathrm{fr}}+\mathrm{1}}=L$. The curvatures κ_j are obtained from solving Eq. (6) individually for every (at this stage, still hypothetical) fragments.

Equation (11) has an explicit dependency to the number of fractures allowed to happen for a given quasi-static state, that is, at a given time. It suggests that the size of x^fr should be a dimension to the minimisation problem. In practice, when considering travelling waves, floes of reasonable size, and the succession of such quasi-static states, at most one single fracture is admissible at a given timestep, which greatly diminishes the numerical cost of the procedure. In what follows, we will thus use N_fr=1, x^fr={x₁}, and we will have

$$\begin{array}{}\text{(13)}& E_{\mathrm{el}}^{\mathrm{s}}\left(x_{\mathrm{1}}\right)=E_{\mathrm{el}}^{<}\left(x_{\mathrm{1}}\right)+E_{\mathrm{el}}^{>}\left(x_{\mathrm{1}}\right)\end{array}$$

with $E_{\mathrm{el}}^{<}\left(x_{\mathrm{1}}\right),E_{\mathrm{el}}^{>}\left(x_{\mathrm{1}}\right)$ the elastic energies of the left (x<x₁) and right (x>x₁) fragments obtained from that single fracture, while the fracture energy reduces to E_fr=G. Hence, we look for

$$\left\{\begin{array}{lr}{x}_{\mathrm{fr}}=\arg {\min }_{x_{\mathrm{1}}\in (\mathrm{0},L)}(E_{\mathrm{el}}^{<}+E_{\mathrm{el}}^{>}),& \left(\mathrm{14}\mathrm{a}\right)\\ F\left(x_{\mathrm{fr}}\right)<\mathrm{0}.& \left(\mathrm{14}\mathrm{a}\right)\end{array}\right.$$

For given ice and wave conditions, fracture search can thus be conducted in a completely deterministic manner. In practice, we proceed by sampling $E_{\mathrm{el}}^{\mathrm{s}}\left(x_{\mathrm{1}}\right)$ regularly on (0,L), ensuring the sampling rate is sufficient to capture its oscillations. We find the set of arguments of the peaks (the local maxima) of this discretised $E_{\mathrm{el}}^{\mathrm{s}}$, which we augment with the bounds {0,L} of the domain, to obtain an ordered sequence of at least two coordinates bounding, two by two, local minima of $E_{\mathrm{el}}^{\mathrm{s}}$. We conduct local minimisation between the bounds using Brent's method (Virtanen et al., 2020). Finally, the smallest of these minima is validated against Eq. (14b). If this condition is verified, its argument is the fracture location x_fr. These steps are summarised in Fig. 2, and a fracture search is illustrated in Fig. 3. Note that in this case, the asymmetry of the total energy profile comes from differences in wave phase at the edges of the floe, and wave attenuation by the ice cover, discussed in Sect. 2.4.

Figure 2 Algorithmic steps behind a fracture search.

Figure 3 Illustration of a fracture search, from a situation corresponding to Fig. 1. In (a), contributions to the system's elastic energy change according to the coordinate of the fracture, and the total energy (that includes the fracture energy) is compared to the energy of the initial, unfractured floe in order to determine whether fracture should occur. The vertical line locates the global minimum, x_fr, of total energy which is, according to our model, where the floe should break. In (b), representation of the deformed floe, before and after fracture at x_fr. The shaded rectangle represents the energy relaxation length, as defined in Eq. (27), centred on x_fr.

2.3.2 Strain criterion

To allow future comparisons, we additionally implement a conventional strain criterion for fracture. Under that formulation, the floe is allowed to fracture if the bending strain ε locally exceeds a prescribed critical strain ε_cr, that is if

$$\begin{array}{}\text{(15)}& \exists x\mid \left|\mathit{\epsilon }\right(x\left)\right|>{\mathit{\epsilon }}_{\mathrm{cr}}\end{array}$$

with

$$\begin{array}{}\text{(16)}& \mathit{\epsilon }\left(x\right)=-{\displaystyle \frac{h}{\mathrm{2}}}{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}w}{\mathrm{d}{x}^{\mathrm{2}}}}\end{array}$$

the maximum (when taking the absolute value) bending strain, here defined as evaluated at the top of the floe.

Typically, if Eq. (15) holds anywhere, it holds on continuous intervals along the floe. We illustrate this in Fig. 4. A second criterion must then be chosen to constrain the fracture. Herein, we arbitrarily choose to consider the global strain extremum (single fracture). We thus have

$$\left\{\begin{array}{lr}{x}_{\mathrm{fr}}=\arg {\max }_{x\in (\mathrm{0},L)}\left|\mathit{\epsilon }\right(x\left)\right|,& \left(\mathrm{17}\mathrm{a}\right)\\ \left|\mathit{\epsilon }\right(x_{\mathrm{fr}}\left)\right|>{\mathit{\epsilon }}_{\mathrm{cr}}.& \left(\mathrm{17}\mathrm{b}\right)\end{array}\right.$$

This criterion can be straightforwardly extended to the case of multiple fracture by considering all the local extrema exceeding the critical value.

2.3.3 Values of fracture parameters

Our two fracture parametrisations rely on two different parameters: the energy release rate G (energy criterion) or the critical strain ε_cr (strain criterion). These parameters have to be measured or estimated from sea ice samples. Sea ice properties can vary greatly based on its history and environmental conditions, such as temperature, brine fraction, or past loading rate. Nevertheless, G and ε_cr are physical quantities that can be measured, and are not completely free parameters. For a detailed compilation of sea ice property measurements, we refer the reader to Timco and Weeks (2010).

Additionally, ice strength depends on the direction of the applied stress. We are here solely interested in failure from bending (mode I, or opening mode; see Saddier et al., 2024), which is compatible with Griffith's model of fracture as well as with wave action. In particular, in the plane strain approximation, ε_cr can be related to the flexural strength σ_f so that ${\mathit{\epsilon }}_{\mathrm{cr}}=\frac{\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}}}{Y}{\mathit{\sigma }}_{\mathrm{f}}$, and G to the fracture toughness K_1c so that $G=\frac{\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}}}{Y}{K}_{\mathrm{1}\mathrm{c}}^{\mathrm{2}}$. Schulson and Duval (2009) compile values of K_1c in the range 75 to 150 kPa m^½. Wei and Dai (2021) measured values down to 26 kPa m^½ for floating samples in the lab, a reduction that could be attributed to temperature or size effects (Dempsey et al., 1999). Reported values of ε_cr are typically in the range of 10⁻⁵ to 10⁻⁴ (Kohout and Meylan, 2008), even though larger value (on the order of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$) have been reported for lab-grown, saline ice (Herman et al., 2018). Sea ice is subject to fatigue, and repeated cyclic loading was shown to lower its apparent flexural strength (Langhorne et al., 1998).

Figure 4 Strain-based fracture parametrisation, for a situation corresponding to Fig. 1. The line represents the maximum bending strain along the floe. The shaded vertical strips indicate where the strain exceeds, in absolute value, a typical critical strain of ${\mathit{\epsilon }}_{\mathrm{cr}}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$. The crosses indicate local extrema, and the vertical line the global extrema.

2.4 Forcing waves

The main focus of this study being ice deformation and fracture, the wave component of the model is kept relatively simple. To align with the linearity hypothesis made on elastic plates, we only consider linear plane waves. This is also in line with previous studies (for instance, Kohout and Meylan, 2008; Dumont et al., 2011; Horvat and Tziperman, 2015; Mokus and Montiel, 2022).

2.4.1 Dispersion relations

For a prescribe angular frequency ω, we derive wavenumbers k from the dispersion relations

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{{\mathit{\omega }}^{\mathrm{2}}}{g}}=k\tanh \left(kH\right)\end{array}$$

in the open-water parts of the domain, and

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{{\mathit{\omega }}^{\mathrm{2}}}{g}}=({\displaystyle \frac{D}{{\mathit{\rho }}_{\mathrm{w}}g}}{k}^{\mathrm{4}}+\mathrm{1}-{\displaystyle \frac{{\mathit{\omega }}^{\mathrm{2}}}{g}}d)k\tanh \left(k\right(H-d\left)\right)\end{array}$$

in the ice-covered parts. We use the single symbol k for brevity, and the appropriate dispersion relation should be understood from context. In the right-hand side of Eq. (19), the term $\frac{D}{{\mathit{\rho }}_{\mathrm{w}}g}{k}^{\mathrm{4}}$ corresponds to the elastic response of the ice cover, while the term $\frac{{\mathit{\omega }}^{\mathrm{2}}}{g}d$ corresponds to its mass-loading response. Whether the former has a significant contribution to the dispersion relation when the ice is heavily fragmented is debated (Sutherland and Dumont, 2018; Dumas-Lefebvre and Dumont, 2023). As it can easily be turned off, we include it here for completeness.

We note that this relation dispersion can be derived by considering the bending of a plate conforming to a fluid surface excited by harmonic waves. This can therefore be seen as a soft coupling of the fluid to the plate. The dispersion relation of plate excited by harmonic waves, without a fluid foundation, would otherwise be ${\mathit{\omega }}^{\mathrm{2}}=\frac{D}{{\mathit{\rho }}_{\mathrm{i}}h}{k}^{\mathrm{4}}$.

2.4.2 Sea state

For any given floe in the domain, a linear monochromatic wave can be parametrised with two complex variables, amplitude $\widehat{a}$ and wavenumber $\widehat{k}$. The modulus $a=\left|\widehat{a}\right|$ denotes the amplitude of the wave at the left edge of the floe, while the argument $\mathit{\varphi }=\mathrm{Ang}\widehat{a}$ denotes the phase of that wave mode at the left edge of the floe. The real part of the wave number, k=Re k, describes wave propagation while its imaginary part $\mathit{\alpha }=\mathrm{Im}\widehat{k}$ describes the spatial rate of attenuation in the direction of propagation. Following Sutherland et al. (2019), we implement a parametrisation with attenuation linear in ice thickness and quadratic in wavenumber, so that

$$\begin{array}{}\text{(20)}& \mathit{\alpha }={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}h{k}^{\mathrm{2}}.\end{array}$$

Other parametric attenuations can easily be added to the current framework; either directly to SWIIFT's codebase¹ or by a user at run time. Attenuation can also be turned off altogether.

A linear polychromatic plane wave can be defined by superposition. The wave state is then

$$\begin{array}{}\text{(21)}& \mathit{\eta }\left(x\right)=\sum _j\mathrm{Im}\left[{\widehat{a}}_j\exp \right({\widehat{k}}_{j}x\left)\right]\end{array}$$

where the subscript j denotes spectral modes. The modal amplitudes can be derived from any spectral density S(ω), using the relationship (Horvat and Tziperman, 2015)

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{a}_j^{\mathrm{2}}=S\left({\mathit{\omega }}_j\right)\mathrm{\Delta }{\mathit{\omega }}_j\end{array}$$

where Δω_j is the width of the angular frequency bin corresponding to the amplitude a_j.

2.4.3 Wave propagation over a finite distance

To allow for the advection of a developing sea into the ice-covered domain, we apply a semi-Gaussian kernel to the sea state. We implement this modification to avoid the non-realistic situation of a fully developed sea appearing under a potentially kilometre-wide MIZ. Therefore, Eq. (21) is modified into

$$\begin{array}{}\text{(23)}& \mathit{\eta }\left(x\right)=\sum _j\mathrm{Im}\left[{\widehat{a}}_j\exp \right({\widehat{k}}_{j}x\left)\right]K_j\left(x\right)\end{array}$$

with

$$\begin{array}{}\text{(24)}& K_j\left(x\right)=\left\{\begin{array}{ll}\exp (-\frac{(x-{\mathit{\mu }}_j{)}^{\mathrm{2}}}{\mathrm{2}{{\mathit{\sigma }}_j}^{\mathrm{2}}})& x\ge {\mathit{\mu }}_j,\\ \mathrm{1}& x\le {\mathit{\mu }}_j.\end{array}\right.\end{array}$$

To each wave mode, we associate a coordinate μ_j. The wave is considered fully developed in the half-plane left of that coordinate. The parameters σ_j control the width of the transition between a fully developed wave, and a near-rest state (as the wave envelop of a given mode is reduced to about 1 % of its maximum at $x={\mathit{\mu }}_j+\mathrm{3}{\mathit{\sigma }}_j$).

2.5 Numerical scheme

We have so far presented our framework for modelling fracture in a quasi-static state. Here, we give more details on how we iterate from a quasi-static state to the next, and summarise the steps leading to evaluating ice floe fractures.

2.5.1 Wave propagation–attenuation

Let τ be a model timestep. Each wave mode in Eq. (23) propagates at a phase speed

$$\begin{array}{}\text{(25)}& c_j={\displaystyle \frac{{\mathit{\omega }}_j}{k_j}}.\end{array}$$

Between time t_n and time t_n+1, the phase ϕ_j of mode j increases by −ω_jτ. The limit μ_j of the fully developed wave advances of a distance c_jτ. Therefore, we iterate in time by updating the values of our ${\widehat{a}}_j$ and μ_j by these quantities. A new quasi-static wave profile η(x) can then be computed across the domain.

2.5.2 Sea ice deformation and fracture

Once the sea surface has been computed, the resulting sea ice deflection w(x) is computed for all individual floes, which are scanned for possible fractures. For a given floe, if ${\mathit{\mu }}_j>L\forall j$ (that is, the wave acting on the floe is fully grown), the deflection can be determined analytically, as developed in Appendix A. Otherwise, we obtain a solution to Eq. (6) with a numerical solver (Virtanen et al., 2020). In turn, the deflection is used to compute the curvature. Depending on the chosen fracture mechanism (energy-based or strain-based), the floes are considered for breakup, as described in Sect. 2.3.1 and 2.3.2.

If using the energy criterion, we evaluate the post-fracture total energy along the discretised floe. We use a peak detection algorithm (Virtanen et al., 2020) to separate intervals of convex free energy (as can be identified in Fig. 3a), onto which Eq. (13) is minimised. If the global minimum among these local minima satisfies Eq. (14b), fracture occurs; these steps are summarised in Fig. 2. If using the strain criterion, we evaluate the bending strain along the discretised floe. Again, a peak detection algorithm is run on −ε²(x) to detect convex intervals, and we conduct local minimisation on these, which is equivalent to maximising $\left|\mathit{\epsilon }\right(x\left)\right|$. If the global minimum satisfies Eq. (17b), fracture occurs.

The input necessary for both parametrisations is thus the floe curvature. In Sect. 2.2 and 2.4, we merely suggest a simple mechanical model to infer this curvature from wave forcing. Other 1D models outputting the curvature field, or actual curvature measurements, can be substituted without having to alter the fracture formalism presented in Sect. 2.3. However, for the fracture parametrisation to be sensible, it is necessary that the mechanical model can be stepped forward in small time increments. Thus, we choose here not to rely on harmonic solutions to the bending problem, such as in Mokus and Montiel (2022), as these rely on the hypothesis that a steady state has been reached in the whole fluid domain. We do so at the cost of neglecting floe bending inertia and relaxing constraints on the fluid itself. In particular, wave scattering induced by different boundary conditions imposed on the fluid when transitioning between open water and ice-covered water regions is not accounted for. A comparison between the two types of solution can be found in Appendix B. We find minor differences in terms of floe curvature (impacting the strain parametrisation) and resulting elastic energy (impacting the energy parametrisation). We compute the ratio of elastic energy (E_el, defined in Eq. (9)), derived from the scattering model, to that same energy derived from SWIIFT, for the case of a polychromatic forcing. The elastic energy is, generally (77.5 % of the ensemble considered), overestimated by SWIIFT; but the distribution of these ratios being skewed, the two models yield, on average, a similar value (mean of the ensemble considered: 1.02, geometric mean: 0.69). Additionally, we do not find these ratios to depend on model parameters such as ice thickness or floe length. We thus conclude that even though differences exist between the two solutions, they are less meaningful than random fluctuations of the wave state.

2.5.3 Timestep selection

Care must be taken when selecting a model timestep, τ. The theoretical upper limit for crack propagation in an elastic, isotropic, and homogeneous material is set by the speed of Rayleigh waves, $c_{\mathrm{R}}=\sqrt{\frac{Y}{\mathrm{2}{\mathit{\rho }}_{\mathrm{i}}(\mathrm{1}+\mathit{\nu })}}$. Using Y=3.8 GPa and ν=0.33, values of the Young's modulus and Poisson's ratio estimated in situ for sea ice (Moreau et al., 2020b), this speed is on the order of $c_{\mathrm{R}}\approx \mathrm{1250}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. As we consider cracks to instantaneously fracture floes through their thickness, we must have $\mathit{\tau }>\frac{h}{c_{\mathrm{R}}}$. For 1 m thick ice, it corresponds to τ>0.8 ms. Fractures propagate faster in stiffer ice, and increasing the Young's modulus (or reducing the thickness) would lower this bound, which must be considered on a case-by-case basis.

However, we also want to keep the timestep small enough that we can detect fractures as soon as it is possible for them to occur, as delaying the onset of a fracture may affect the length of the resulting floes. Therefore, we aim to keep the ratio of the progression of the wave front to the wave amplitude small. In the monochromatic case, with phase speed c, it translates to having $\frac{\mathit{\tau }c}{a}<r\iff \mathit{\tau }<r\frac{ak}{\mathit{\omega }}$, with r<1. Setting $r=\frac{\mathrm{1}}{\mathrm{5}}$ ensures sufficient convergence. An analogous relationship can be derived for polychromatic cases, by substituting the amplitude by the significant wave height H_S of the spectrum, and the phase speed c by the maximum phase speed of within the sampled spectrum, so that $\mathit{\tau }<r\frac{H_{\mathrm{S}}}{{max}_j{c}_j}$.

2.6 Example of time simulation

The study of floe size distributions in relation to ice or wave parameters is out of the scope of this study. However, as an illustration of the capabilities of SWIIFT, we present in this Section the result of a single simulation.

We initialise the domain with a single floe of length L=600 m, thickness h=50 cm, Young's modulus Y=4 GPa, and Poisson's ratio ν=0.3. We choose to parametrise fracture with the energy criterion, and set the fracture toughness to K_1c=100 kPa m^½, which together with the other mechanical parameters corresponds to an energy release rate G=2.275 J m⁻².

This floe is forced with waves issued from a (one-parameter) Pierson–Moskowitz spectrum, with significant wave height H_S=0.5 m (corresponding to a peak period of 3.84 s), truncated to the period interval T∈1 to 15 s. We discretise this spectrum onto 50 linearly spaced frequency bins, and thus obtain 50 tuples of amplitudes and wavelengths, which we complete with 50 initial phases sampled from a uniform distribution in 0 to 2π rad. Spatial attenuation is parametrised as defined in Eq. (20). We initialise the growth kernels K_j with standard deviations σ_j equal to the wavelengths, and means μ_j equal to three respective standard deviations upstream from the floe. At time t=0 s, the magnitude of the surface perturbation at the left edge of the floe, issued from wave superposition as defined in Eq. (23), is about 0.2 mm.

We set the timestep $\mathit{\tau }=\frac{\mathrm{1}}{\mathrm{5}}\frac{H_{\mathrm{S}}}{{max}_j{c}_j}=\mathrm{8.58}\mathrm{ms}$. We run the simulation for 120 s; the first fracture occurs at t=9.097 s, the last one at t=105.497 s. We present results of this fracture experiment in Fig. 5, showing a snapshot of the simulated fluid and ice displacement along with the evolving number and lengths of the simulated fragments. A video of the simulation is available as supplementary material (Mokus, 2025a).

Figure 5 Snapshots of a fracture experiment. In (a), view of the domain at t=60 s. The continuous, dark line represents the fluid surface (η(x)), and the discontinuous, lighter lines the vertical displacements (w(x)) of individual floes. The marks along the bottom spine indicate the boundaries between fragments; the last 80 m, at the right of the domain, have not yet been affected by the waves. Note that the vertical scale is greatly exaggerated: the aspect ratio of the graph, in physical units, is $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$. Because of the thickness of the lines, some floes appear to overlap, they actually do not. In (b), horizontal bars show the extent of individual floes. The height of the bars indicates the order of the floe in the array, and each group of bars, or “stair”, corresponds to a snapshot. The time of the snapshots are indicated on the y-axis, and darker colours correspond to later times. In (c), we show size distributions as swarmplots, omitting the rightmost fragment. Each dot corresponds to a length as indicated by the x-axis, and within a group, the y-axis only serves to separate dots. Vertical clusters thus indicate a concentration of observations around the corresponding length. From t=0 to 120 s, there are respectively 1, 6, 27, 52, 58, 59, and 60 fragments.

3 Numerical experiment

To evaluate the capabilities of our model, in particular, validate the energy-based fracturing approach and highlight the difference between energy and strain criteria, we choose to replicate breakup experiments conducted at the laboratory scale on a material that served as an analogue for solid, cohesive ice (Auvity et al., 2025). These focused on quantifying the onset of breakup, by progressively increasing the amplitude of a forcing stationary wave, at different frequencies. The experiment setup was as follows: a water tank of length 80 cm and depth 11 cm was covered with a brittle layer of varnish, with thickness on the order of 100 µm. The layer was detached from the walls of the tank prior to the experiment. Stationary surface waves were generated with a wave maker. A one-dimensional profilometry system and image-processing method were used to extract the wave properties (frequency, amplitude, wavenumber) and determine when fracture occurred. This work is similar to that of Saddier et al. (2024), who also conducted wave-induced fracture experiments on an analogue material at the laboratory scale, under stationary but also progressive forcing. However, in their experiment, the material is a granular raft hold together by capillary forces more than a continuous solid, and breaks because of viscous stress rather than because of bending stress. The former is directly relevant for representing the disintegration of an already fragmented and granular sea ice, that has already been broken up or that is in a consolidation phase (transition from frazil to grease ice). As our work focuses on the fracturing of a solid and cohesive ice cover that we treat as a continuous elastic medium, rather than on the disintegration of a granular ice of low cohesion, we favoured the work of Auvity et al. (2025) for our comparisons.

As we aim to replicate this experiment, in what follows, we will be using a standing wave forcing, and turn off any attenuation. We use our model to determine, for prescribed wavenumbers and material properties inherited from these laboratory experiments (listed in Table 1), the critical amplitude a_cr at which the material starts to fracture. We do so using our energy formulation. Our model being linear, the amplitude directly controls the deflection of the plate, and thus, its curvature and resulting elastic energy. It is therefore an intuitive quantity to control the outputs of the model, as well as a quantity that was measured experimentally.

The critical amplitude can then be related to a critical curvature κ_cr by evaluating Eq. (8) at the coordinate of the fracture. In turn, κ_cr can be used to derive a critical strain, using Eq. (16). We do not run separate experiments based on a strain criterion, as per the results of these authors, a critical strain independent of the wave forcing does not seem to exist for this material and therefore cannot be prescribed in numerical experiments. Even so, the results of energy-based simulations allow us to draw conclusions on the relevance of this type of criterion. These are discussed in Sect. 5.

3.1 Length scales

To replicate the experimental protocol of Auvity et al. (2025), we consider only monochromatic stationary forcings, so that η(x)=asin (k_nx) with $k_n=\frac{n\mathit{\pi }}{L}$. The symbol L represents both the length of the plate, and of the domain. The positive integer n is the harmonic number. The wave tank we simulate is short enough for attenuation to be considered insignificant.

We define two additional lengths: the flexural length

$$\begin{array}{}\text{(26)}& L_{\mathrm{D}}=({\displaystyle \frac{D}{{\mathit{\rho }}_{\mathrm{w}}g}}{)}^{\mathrm{1}/\mathrm{4}},\end{array}$$

and the relaxation length

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{L}_{\mathit{\kappa }}& ={\displaystyle \frac{{\int }_{\mathrm{0}}^{x_{\mathrm{fr}}}(x_{\mathrm{fr}}-x)\left[\mathit{\kappa }\right(x)-{\mathit{\kappa }}^{<}(x){]}^{\mathrm{2}}\mathrm{d}x}{{\int }_{\mathrm{0}}^{x_{\mathrm{fr}}}\left[\mathit{\kappa }\right(x)-{\mathit{\kappa }}^{<}(x){]}^{\mathrm{2}}\mathrm{d}x}}\\ & +{\displaystyle \frac{{\int }_{x_{\mathrm{fr}}}^L(x-x_{\mathrm{fr}})\left[\mathit{\kappa }\right(x)-{\mathit{\kappa }}^{>}(x){]}^{\mathrm{2}}\mathrm{d}x}{{\int }_{x_{\mathrm{fr}}}^L\left[\mathit{\kappa }\right(x)-{\mathit{\kappa }}^{>}(x){]}^{\mathrm{2}}\mathrm{d}x}}.\end{array}\end{array}$$

The former is a natural length scale of our system, appearing in the bending equation (Eq. 6a), and relates the flexural rigidity of the plate (that resists bending) to the reaction of the fluid it rests upon (that sustains bending). The latter gives a measure of the distance over which the curvature of post-fracture fragments is different from the curvature of the original floe that gave rise to these fragments. We show an example of this in Fig. 3b.

We introduced the symbols κ^<(x) and κ^>(x) to denote the curvature of the left and right post-breakup fragments, with $x\in [\mathrm{0},L]$. By definition, κ^< (respectively κ^>) exists only for $x\in [\mathrm{0},x_{\mathrm{fr}}]$ (respectively $x\in [x_{\mathrm{fr}},L]$). We choose this integral definition of L_κ because the differences in pre- and post-breakup curvatures is well-represented (when moving away from the fracture location) by a damped sine with oscillation period and attenuation rate $\sqrt{\mathrm{2}}{L}_{\mathrm{D}}$, that is,

$$\begin{array}{}\text{(28)}& \mathit{\kappa }\left(x\right)-{\mathit{\kappa }}^{>}\left(x\right)\sim \sin \left({\displaystyle \frac{x}{\sqrt{\mathrm{2}}{L}_{\mathrm{D}}}}\right)\exp (-{\displaystyle \frac{x}{\sqrt{\mathrm{2}}{L}_{\mathrm{D}}}}),\end{array}$$

which ensues from the shape of the solution presented in Appendix A1. Therefore, except at the floe boundaries where curvature is 0 m⁻¹ (as imposed by the boundary condition, Eq. 6b), the whole length of the initial floe may participate in releasing energy. For long enough waves, the relaxation length tends to $\sqrt{\mathrm{2}}{L}_{\mathrm{D}}$, that is ${lim}_{k\to \mathrm{0}}{L}_{\mathit{\kappa }}=\sqrt{\mathrm{2}}{L}_{\mathrm{D}}$. This can be shown analytically by assuming Eq. (28).

We thus have three typical horizontal length scales:

• The wavelength $\mathit{\lambda }=\frac{\mathrm{2}\mathit{\pi }}{k}$, imposed by the wave forcing, and linearly tied to the domain length L through the harmonic number n, so that $L=\frac{n\mathit{\lambda }}{\mathrm{2}}$.

• The flexural length L_D, that depends on the properties of the material, the density of the fluid, and gravity. Only the former are varied in this study, with stiffer, thicker materials having a longer L_D.

• The relaxation length L_κ, that quantifies the distance over which fracture modifies the system.

As we consider short wavelengths and a very thin plate, capillarity effects could in principle be important. However, because the flexural length of the material exceeds its capillary length, these are negligible, which Auvity et al. (2025) verified experimentally. Therefore, we will not consider them either, and the dispersion relation we will be using is given in Eq. (19), dominated by the term in ${L_{\mathrm{D}}}^{\mathrm{4}}$.

3.2 Linearity limitation

The analogue material used in the laboratory experiments of Auvity et al. (2025) requires nonlinear waves (ak≈0.14) for fracture to occur. As neither nonlinear plate nor non-linear waves are represented by our numerical model, we have to relax this condition, typically quantified by the wave slope ak, to observe fracture at all. Thus, we set the upper bound of our dichotomic searches so that ak≤0.5, which places us out of the linear framework our model relies on. As here, we are qualitatively showcasing the behaviour of our model rather than quantitatively exploiting the results, we deem this limitation to be inconsequential. For thickness and Young's modulus typical of sea ice, fracture in our model does happen in a linear regime, as illustrated in Fig. 3, where ak=0.015.

Note that we define wave slope with respect to the wave propagating underneath the elastic plate and not with respect to the free surface waves. For a given time period, hydroelastic waves with dispersion relation Eq. (19) are typically longer than free surface gravity waves with dispersion relation Eq. (18), making the former slightly less steep.

4 Results

The results presented in this section focus on detecting a fracture threshold using our energy formalism and comparing this threshold to that obtained in the laboratory experiment of Auvity et al. (2025). To do so, we use in our simulations the material parameters issued from Auvity et al. (2025). Those are given in Table 1. We do not tune model parameters. As we are interested in detecting the fracture threshold, and our model does not have a fatigue term, we work with strictly unrelated quasi-static states, and we find the critical amplitude a_cr systematically by dichotomic search. In Sect. 4.1, we detail the influence of varying the wavenumber exclusively. Then, in Sect. 4.2, we replicate the same protocol, while also varying the mechanical parameters of the plate.

4.1 Reference case

We start by illustrating the response of our model, for a range of prescribed wavenumbers, with four quantities: the normalised position of the fracture $\frac{x_{\mathrm{fr}}}{L}$, the critical amplitude and curvature, and the relaxation length. These are presented in Fig. 6. To exemplify the deviation between our computed deflection field and the forcing fluid surface, we choose here the case of the third harmonic, so that $L=\frac{\mathrm{3}\mathit{\pi }}{k}$. We thus obtain curvature profiles that are symmetric² with respect to $x=\frac{L}{\mathrm{2}}$, with three antinodes, which can be seen in Fig. 7. From left to right, the first and third antinodes (close to the left and right edges of the plate, respectively) are more influenced by the boundary conditions than the second one (located at the middle of the plate).

Table 1 List of model parameters and their values. Parameters followed by an asterisk are inferred from other fixed parameters.

Figure 6 Relationships between the nondimensionalised wavenumber and normalised (with respect to plate length) fracture location (a), critical amplitude (b), critical curvature (c), and relaxation length (d). Model parameters are provided in Table 1. In (a), three horizontal dashed lines represent the asymptotes $\frac{x_{\mathrm{fr}}}{L}=\frac{\mathrm{1}}{\mathrm{6}}$, $\frac{\mathrm{1}}{\mathrm{3}}$, and $\frac{\mathrm{1}}{\mathrm{2}}$. In (b) to (d), vertical lines show delimitations between regions corresponding to different behaviours of fracture location, observed on (a). The regions are numbered in (c). In the second region of (b), the triangle of height twice its horizontal base (in loglog space and data units) gives an indication of the slope. In (d), an horizontal dashed line represents the asymptote $\sqrt{\mathrm{2}}{L}_{\mathrm{D}}$.

In what follows, we multiply k by the flexural length, so that to obtain the (dimensionless) wavenumber kL_D. This allows exploring the model behaviour between two limits: small kL_D values thus correspond to longer waves and lower flexural rigidity (the plate conforms to the fluid), while high values correspond to shorter waves and higher flexural rigidity (the plate is non-deformable). We divide Fig. 6b–d different behaviours of the fracture location as predicted by the energy criterion, identified in Fig. 6a. These regions are separated by kL_D=0.1638, 0.3275, and 0.7578. They correspond, for increasing kL_D, to fracture happening in the middle of the floe (the second curvature antinode); fracture happening close to the first or third curvature antinode; fracture again happening in the middle of the floe; and fracture uncorrelated from any curvature extremum. Because of the particular wave forcing imposed in our model configuration, the free energy profile is symmetrical with respect to the middle of the plate. Therefore, it is not numerically possible to discriminate between an energy minimum happening at x_fr or L−x_fr, and in Fig. 6a, we only show the branch corresponding to $x_{\mathrm{fr}}\le \frac{L}{\mathrm{2}}$.

We note that, in regions 1 to 3, fracture predicted by the energy criterion does not systematically happen at the global curvature extrema. In this third harmonic case, the global extremum is in the middle of the floe, except in the band $k{L}_{\mathrm{D}}\in [\mathrm{0.178},\mathrm{0.357}]$. The overlap with the region 2, defined by $k{L}_{\mathrm{D}}\in [\mathrm{0.1638},\mathrm{0.3275}]$, is thus not one-to-one. Additionally, in region 2, fracture does not happen at the antinode, but in its vicinity. For kL_D<0.242, fracture is on the left of the first antinode, and for kL_D>0.243, to its right (the situation is reversed for the third antinode). In region 4, fracture happens far from either antinode. It seems that for increasing kL_D, $\frac{x_{\mathrm{fr}}}{L}\to \frac{\mathrm{1}}{\mathrm{3}}$, which corresponds to a node of the forcing.

In Fig. 7, we show examples of the behaviour of the free energy and of the along-plate curvature for these different regions. We compare the latter to the “conforming” curvature, ${\mathit{\kappa }}_{\mathrm{conf}}\left(x\right)=\frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{x}^{\mathrm{2}}}{w}_{\mathrm{conf}}$, with $w_{\mathrm{conf}}\left(x\right)=\frac{a}{\mathrm{1}+(k{L}_{\mathrm{D}}{)}^{\mathrm{4}}}\sin \left(kx\right)$ the associated displacement stemming from the fluid surface. The term at the denominator ensures that it satisfies Eq. (6a), and for long waves, ${lim}_{k\to \mathrm{0}}{w}_{\mathrm{conf}}=\mathit{\eta }$. In Fig. 8, we show for the same examples the floe deflection, compared to the forcing amplitude, and we indicate the relaxation length. To aid comparison, we normalise the along-floe coordinate by the floe length.

In the first two regions, both corresponding to small wavenumbers, the difference between curvature and conforming curvature (Fig. 7a, b) is noticeable only in the immediate vicinity of the edges of the domain. Minima of free energy correspond roughly with extrema of curvature, and the floe deflection (Fig. 8a, b) follows the fluid surface. In the first region, and to a lesser degree in the second region, the critical curvature (Fig. 6c) varies little, although a slight positive trend exists. Our strain-based and energy-based criteria in these two regions would therefore predict virtually similar fractures. The critical amplitude (Fig. 6b) varies with the inverse of the squared wavenumber, that is, with the square of the wavelength. The relaxation length (Fig. 6d) is almost constant and tends to $\sqrt{\mathrm{2}}{L}_{\mathrm{D}}$ from below for decreasing wavenumbers. As the floe length is, in the case of stationary wave forcing, inversely proportional to the wavenumber, the relaxation length normalised by the floe length increases with the wavenumber, and is not constant across the different panels of Fig. 8.

In the third region, curvature and conforming curvature (Fig. 7c) are now dissimilar between the left (respectively right) edge and the left (respectively right) antinode. The free energy still shows three troughs, but the trough at $x=\frac{L}{\mathrm{2}}$ is now clearly more pronounced than the other two. There are still three distinct deflection extrema (Fig. 8c), synchronised with the forcing wave, and the amplitude of deflection of the floe is slightly smaller than the forcing amplitude. From kL_D≈0.3057, we locally (around the two positive deflection antinodes) have $\mathit{\eta }-w>h-d$. As h−d corresponds to the freeboard of the floe at rest, this suggests parts of the deformed floe are immersed. This takes place close to the transition from the second region, which happens at kL_D=0.3275. The non-zero deflection near the edges shows that deflection is now significantly different from the sine forcing. In terms of the occurrence of fracture, this region corresponds to a sharp increase in critical curvature, incompatible with a strain-based (that is, constant critical curvature) criteria. The relaxation length, however, is still practically constant with kL_D, and a good indicator of the zone over which pre- and post-fracture modelled deflection differ. An inflexion of the critical amplitude decrease rate is also visible. In regions 1 to 3, fracture locations near antinodes and the shape of post-fracture deflections are consistent with mode I fracturing.

As kL_D increases, the length of the plate diminishes relatively to its flexural length. The impact of the boundary conditions on the deflection profile is therefore amplified. This effect is sizeable in the fourth region (Fig. 8d): curvature and conforming curvature (Fig. 7d) no longer match anywhere along the plate. This is despite staying in a regime where L_D≪L, as $k{L}_{\mathrm{D}}=n\mathit{\pi }\frac{L_{\mathrm{D}}}{L}\approx \mathrm{10}\frac{L_{\mathrm{D}}}{L}$. The central free energy trough has separated into two distinct troughs corresponding to global minima, no longer in phase with curvature extrema. This separation corresponds to the transition from the third region. The fourth region also shows a drop in critical curvature (Fig. 6c), which has been monotonically increasing with kL_D thus far. However, the maximum curvature $\mathit{\kappa }(L/\mathrm{2})$ keeps increasing irregularly with kL_D, as can be seen by comparing Fig. 8c and d. The critical amplitude, which seems to plateau on the right of the third region (Fig. 6b), is singular at the transition, then diminishes again before increasing irregularly. Additionally, the edges of the fragments no longer mirror each other. There is a significant post-fracture discontinuity in deflection, which is a characteristic of region 4, and not consistent with bending (mode I) fracture, but reminiscent of a sliding (mode II) or tearing (mode III) fracture. A minimum of critical curvature is reached for kL_D=1, which also corresponds to a maximum of relaxation length. For n≥2, there is a single positive kL_D, quickly converging to 1, for which the slope at the edges of the floe vanishes, that is, $(\frac{\mathrm{d}w}{\mathrm{d}x}{)}_{x=\mathrm{0},L}=\mathrm{0}$. It corresponds to the two outside-most deflection antinodes vanishing, leaving only the internal ones. It can be seen from the deflection profile in Fig. 8d, that compared to Fig. 8a–c, the slope at the edges has changed sign, and that a single antinode (at the centre of the floe) remains.

Figure 7 Examples of free energy (left vertical axes, green lines) and curvature (right vertical axes, grey lines) profiles for the four regions identified in the text. A sine curvature is shown in dashed lines as a comparison to the curvature derived from Eq. (8). The fracture locations, determined from the energy criterion, are shown with vertical lines: the energy profiles being symmetrical with respect to the middle of the plates, $F\left(x_{\mathrm{fr}}\right)=F(L-x_{\mathrm{fr}})$. The dashed horizontal lines show the zero-curvature reference.

Figure 8 Examples of pre- and post-fracture floe deflection profiles for the four regions identified in the text, normalised by the critical amplitude. The relaxation lengths, centred on the fracture locations, are shown with shaded rectangles. We only consider the left fracture location, where applicable. The dashed horizontal lines show the zero-deflection reference.

4.2 Influence of mechanical parameters

We further investigate the response of our model to varying mechanical parameters, by reproducing the analysis presented in Sect. 4.1 for an ensemble of (h,Y) pairs with 128 members. Doing so, we aim to reproduce the internal variability that stems from laboratory conditions in the experiments of Auvity et al. (2025). We generate this ensemble through Latin hypercube sampling, and enforce that the two variables are independent with normal marginal densities of prescribed means 100 µm and 70 MPa, and prescribed standard deviations 20 µm and 14 MPa, respectively; we show the joint density of our sample in Fig. 9. The resulting distribution of flexural rigidities is positively skewed, with mean $\mathrm{7.88}\times {\mathrm{10}}^{-\mathrm{6}}$ Pa m³ and median $\mathrm{6.93}\times {\mathrm{10}}^{-\mathrm{6}}$ Pa m³.

Figure 9 Joint density and marginal densities of our (h,Y) ensemble.

We further impose the relation

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\mathrm{2}G{h}^{\mathrm{2}}}{D}}=C\end{array}$$

between the energy release rate, the Young's modulus, and the thickness, as derived by Auvity et al. (2025), setting the dimensionless material constant parameter $C=\mathrm{2.8}\times {\mathrm{10}}^{-\mathrm{4}}$. This expression serves as a proxy establishing a value for G, which is poorly constrained experimentally. The other parameters are kept fixed at the values presented in Table 1. We also keep the same constraint on the wave slope, requiring a_crk≤0.5: depending on the precise values assumed by the thickness and Young's modulus, the interval of wavenumbers that leads to fracture may vary.

4.2.1 Comparison to experimental data

We show numerical results in Fig. 10 (colour-coded lines), which we compare to experimental data from Auvity et al. (2025) (circle and square markers). We obtain results similar to those presented in Sect. 4.1. The critical amplitude profiles do not depend on the mechanical properties of the simulated material, up to a multiplicative constant, that increases with flexural rigidity (or, equivalently, flexural length).

This fact extends to the other variables shown in Fig. 6. Within the ranges of mechanical parameters explored, which are in agreement with the values and internal variability estimated by Auvity et al. (2025), the order of magnitude of the critical amplitude and its decreasing tendency with increasing kL_D agree between the simulations and the laboratory experiments. However, and as expected, this agreement is only qualitative. Indeed, we recall that in the experimental setting, fracture was only obtained with nonlinear waves. It is likely the reason why the critical amplitudes measured by Auvity et al. (2025) varies as k⁻¹, not as k⁻² as we find numerically. The differences between experimental and numerical critical amplitudes are therefore deepened at small wavenumbers.

Figure 10 Relation between the dimensionless wavenumber and the critical amplitude, for varying thickness and Young's modulus. Lines are numerical results, with the colour scale indicating the flexural rigidity, that combines the two varying parameters. Grey dots are experimental results from Auvity et al. (2025). Pink squares are the same experimental results, with a different horizontal scaling, as described in the text.

We also note that the way we define the relaxation length L_κ in our linear waves simulations differs from the definition of Auvity et al. (2025), who used the full width at half-maximum (hereinafter, FWHM) of pre-fracture curvature in their nonlinear experiments. On Fig. 10, we thus also represent the experimental critical amplitudes as a function of $\frac{L_{\mathrm{D}}}{\mathrm{FWHM}}$ (pink squares), which horizontally shifts the experimental points, in an attempt to correct the discrepancy between their nonlinear waves forcing and our linear model.

We do not adopt their definition, as it would be incompatible with our region 4 results, where fracture does not happen around curvature peaks or troughs, while experimental fractures always happened in the vicinity of a deflection antinode. If we were to apply it to regions 1 to 3, we would obtain something very similar to the FWHM of sin (kx), that is $\frac{\mathrm{2}\mathit{\pi }}{\mathrm{3}k}=\frac{\mathit{\lambda }}{\mathrm{3}}$. Because of their nonlinear wave forcing, Auvity et al. (2025) measured FWHMs that varied like $\frac{\mathit{\lambda }}{\mathrm{12}}$. Bending was thus concentrated in a smaller fraction of their forcing wavelengths, and $\frac{L_{\mathrm{D}}}{\mathrm{FWHM}}$ can be seen as an alternative, rescaled dimensionless wavenumber. Using this definition, we can improve the overlap between experimental and numerical results, with experimental points falling in regions 1 to 3, where both model end experiments show the critical curvature depends on the forcing wavelength, precluding a constant strain threshold.

Figure 11 Relationship between the nondimensionalised wavenumber and two dimensionless quantities: the energy dissipation length scaled by the flexural length (a), and squared critical curvature scaled by thickness and energy dissipation length (b). We show ensemble averages, for four harmonic numbers. For a given harmonic number, ensemble members are virtually identical, with coefficients of variation well below $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{3}}$ where the means are non-zero.

4.2.2 Impact of harmonic number and dimensionless quantities

So far, we focused on the harmonic number n=3. However, for large enough kL_D, the response of the model depends on n, and therefore on the geometry of the domain. In particular, for even n, fracture in region 4 happens systematically in the middle of the floe. Due to the symmetry property of the forcing, this means our model predicts fracture with ${\mathit{\kappa }}_{\mathrm{cr}}=\mathrm{0}{\mathrm{m}}^{-\mathrm{1}}$, inconsistent with bending fracture but reminiscent of shearing or tearing fracture. The fundamental configuration, with n=1, is a particular case. The maxima of deflection, curvature, and free energy happen at $x=\frac{L}{\mathrm{2}}$ independently of kL_D. We illustrate this difference in Fig. 11 by presenting the relationships between the nondimensionalised wavenumber and two dimensionless quantities, for different harmonic wavenumbers. We choose these two quantities because they exhibit the remarkable property of depending on the dimensionless wavenumber, but not on the individual variations of the mechanical parameters.

The first of these quantities, shown in Fig. 11a, is the relaxation length L_κ (as shown in Fig. 6d for the fixed parameters experiment), normalised by the augmented flexural length $\sqrt{\mathrm{2}}{L}_{\mathrm{D}}$. We retrieve, independently of the harmonic number, the limit ${lim}_{k\to \mathrm{0}}{L}_{\mathit{\kappa }}=\sqrt{\mathrm{2}}{L}_{\mathrm{D}}$ already stated in Sect. 3.1. Behaviours depending on the harmonic number emerge from kL_D≈0.4. If all curves show a downward trend in what corresponds to kL_D in regions 2 and 3, this trend is more pronounced for smaller numbers, in particular $n\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$. The second striking difference, is that between even and odd harmonic numbers. There is a discontinuity at the transition from region 3 to region 4 for even numbers, with an upward jump preceding a sustained downward trend. This transition is continuous for n=3. There is no region 4 behaviour for n=1, as in this configuration, fracture always happens at $\frac{x_{\mathrm{fr}}}{L}=\frac{\mathrm{1}}{\mathrm{2}}$.

The second dimensionless number, shown in Fig. 11b, can be built as the product of two distinct dimensionless quantities: critical curvature multiplied by thickness (that is, twice the critical strain), and critical curvature multiplied by relaxation length. Notably, both these quantities do depend on thickness and Young's modulus, without showing the ordering critical amplitude does in Fig. 10. However, their product, ${{\mathit{\kappa }}_{\mathrm{cr}}}^{\mathrm{2}}h{L}_{\mathit{\kappa }}$, only depends on kL_D. This quantity was also derived by Auvity et al. (2025), who interpreted it as a constant independent of the wave forcing. We do not replicate this result outside of regions 1 and 2, that is, the wavenumber band where neither critical curvature nor relaxation length vary. As in Fig. 11a, the different curves are indistinguishable within these two regions (that is, at small kL_D), and a discontinuity exists for even harmonic numbers between regions 3 and 4, as for these, the critical curvature drops to 0 m⁻¹ in region 4.

Finally, it can be seen that the upper bound of the range of kL_D that sees fracture happens depends on n. It first increases with n but peaks for n=3, and then decreases. This is despite keeping the same a_crk≤0.5 criterion on the wave slope. The lower bound, however, does not change and keeps the value kL_D=0.1167. The differences between the different harmonics in region 4 are explained by the loss of deflection extrema near the boundaries as kL_D increases, which have dissimilar effects on the deflection profile in this region for different n. For n>2, there exists a single kL_D for which the slope of the deflection at the edges of the plate, $(\frac{\mathrm{d}w}{\mathrm{d}x}{)}_{x=\mathit{\{}\mathrm{0},L\mathit{\}}}$, cancels. It converges exponentially towards 1, so that noting it (kL_D)₀, we have $\log |\mathrm{1}-(k{L}_{\mathrm{D}}{)}_{\mathrm{0}}|\sim -n$. The cancellation of the slope conveys the transition from a deflection profile with n antinodes to one with n−2 antinodes. When kL_D keeps increasing, the higher n, the higher the variety of behaviours shown by w. However, for large enough kL_D, $w\left(x\right)\to \frac{\mathrm{2}a}{n\mathit{\pi }}$ for odd n, and $w\left(x\right)\to \frac{\mathrm{2}a}{n}(-\frac{\mathrm{2}x}{L}+\mathrm{1})$ for even n. These are rigid motions, independent of the forcing, corresponding respectively to heave (translation) or pitch (rotation).

5 Discussion and conclusion

We have developed a versatile, lightweight one-dimensional model that simulates the time-dependent fracture of sea ice by waves. This model has the particularity of solving directly for the deflection of a floe caused by the competition between buoyancy and gravity; instead of solving for waves scattered by the presence of the floe at the ice–fluid interface, and assuming that the deflection follows that interface. We can thus use the model to continuously explore the response of a floe to bending between two limits: an elastic plate conforming exactly to a fluid foundation, and an undeformable plate. We have implemented in this model two fracture criteria. One, compatible with continuum fracture mechanics, is based on looking for a global post-fracture energy minimum and comparing it to the pre-fracture energy state to determine whether fracture should occur, and is a novelty of this model. The other is compatible with the more common hydroelastic approach applied to sea ice, based on locally comparing strain to a prescribed, constant threshold.

The present study is centred on presenting the theoretical and numerical aspects of the model itself and validate/invalidate the energy-based fracture criterion. We apply the model to an analogue material used in the laboratory to study wave-induced ice fracture, in the specific setting of monochromatic stationary waves. Because no constant critical strain threshold was observed during these experiments, but a relationship between energy release rate and other parameters of our model exists, we focus on investigating the energy criteria. Even in this particularly simplified configuration, the response of the model in terms of critical amplitude or curvature is not straightforward.

Our results indicate that the critical curvature derived from an energy-based fracture depends on the forcing wave, contradicting the existence of a universal critical strain. This was also observed in the laboratory (Auvity et al., 2025). In the (dimensionless) wavenumber band where our results overlap with the experimental data from Auvity et al. (2025), we obtain comparable critical amplitudes. However, we are not able to replicate their scaling for low kL_D. Additionally, we obtain that for large enough kL_D (region 4), the two criteria, energy-based and critical strain-based, diverge on the predicted fracture location, in that energy-predicted fracture is uncorrelated from curvature, or strain, extrema. The fracturing behaviour in region 4 is inconsistent with bending fracture, and suggests out-of-plane shear or in-plane shear fracturing. The latter would describe fracture propagating perpendicularly to the direction of wave propagation (that is, as someone tearing up a sheet of paper), in contradiction with the invariance hypothesis made on the modelled plate and therefore cannot be represented in the current 1D model.

A possible explanation for the different relationship between critical amplitude and wavenumber observed experimentally ($a_{\mathrm{cr}}\sim k^{-\mathrm{1}}$) and numerically ($a_{\mathrm{cr}}\sim k^{-\mathrm{2}}$, for kL_D≲0.6) may be that experimentally, for fracture to happen, the material considered (varnish) required nonlinear wave forcing, which our model does not represent. In the case of ice, the linearity assumption is, however, valid. For values corresponding to a similar experiment conducted on fresh water ice by the same team (Auvity Baptiste and Zanchi Vasco, personal communication, 2025), that is thicker and stiffer than their varnish, our results in terms of critical amplitude as a function of dimensionless wavenumber are, up to multiplicative constants, identical to those presented here, and the wave slopes required to obtain fracture are typically of the order of 0.02; well within the linear regime. Therefore, increasing the numerical complexity of the model to accommodate nonlinear plate behaviour seems unnecessary at this stage. Another explanation is that, while our model considers homogeneous plates with constant Young's modulus, the material engineered by Auvity et al. (2025) is obtained by layering. Because of introduced vertical inhomogeneities, this process is likely to introduce a dependency of the Young's modulus to the obtained thickness. Further analysis of this new dataset, and in particular whether the dependency of the curvature at failure on the wavenumber exists, is ongoing.

The key features of our model are that bending is driven exclusively by the along-plate variation of buoyancy, which cannot be resolved by hydroelastic models, and that fracture can be controlled by an energy criterion integrated over the entire floe. For high kL_D values, that is either large wavenumber or a stiff elastic plate, this leads to physically questionable behaviours, such as submergence, and post-breakup deflection discontinuity across the fracture (region 4) or even fracture at zero-curvature (for even harmonics). We note that the possibility of submergence is the direct consequence of the weak, one-way coupling between fluid and plate, as we only represent the response of the plate to the fluid, while ignoring the feedback response of the fluid. This one-way coupling is a trade-off allowing us to maintain the theoretical and numerical complexity low. As a rationality check, we verify that the bending energy, transmitted to the plate by the fluid, is orders of magnitude less than the gravitational potential energy of the fluid: there is thus no unaccounted energy leaks into the plate. As none of this energy is returned to the fluid, it is likely we overestimate the likeliness of fracture.

Here, we have described the simulated behaviour as a function of kL_D by distinguishing the results in 4 different regions, based on where the fracture happens. The kL_D thresholds being regions should however be regarded carefully, as they depends slightly on the harmonic number of the forcing, and might be a feature of stationary forcing. At field scale, with waves propagating within the ice cover, the fracture front follows the wave front, in such a way that fragments are typically smaller than the dominant wavelength (Dumas-Lefebvre and Dumont, 2023). Therefore, the increased complexity in the model behaviour at high harmonic numbers may not be representative of natural conditions.

More experimental data is needed to confirm or infirm the behaviour of the model in the large kL_D band, and whether a forcing-dependent trend for critical curvature exists in the case of ice. Previous wave tank fracture experiments (Dolatshah et al., 2018) showed bending failure typically happening at kL_D≈1, though the uncertainties on thickness and Young's modulus are quite large. Other values of Young's modulus reported for such experiments, in the low MPa range (Herman, 2018; Passerotti et al., 2022), seem inconsistent with a cohesive, solid sheet of ice, as represented in our model. Nevertheless, these authors did observe fracture in the range $k{L}_{\mathrm{D}}\in [\mathrm{0.29},\mathrm{0.54}]$ and for kL_D=0.29, respectively. Voermans et al. (2020) compiled a list of studies of wave-induced breakup observations. Mechanical parameters were, for the most part, not measured, but they suggest estimations based on known empirical relations. Following their methods, we can generate ensembles of kL_D wavenumbers, that we find lying in the range $k{L}_{\mathrm{D}}\in [\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{2}},\mathrm{3.1}]$ for both breaking and non-breaking cases. In the case of realistic wave forcing, with material parameters representative of first-year sea ice, the peak of wave energy occurs in kL_D bands corresponding to what we identified as region 3 and 4. However, higher-frequency waves are also the ones most effectively attenuated by the ice cover, so that they contribute less to fracture.

We acknowledge our results are a first step towards the validation of the fracture formalism we propose. Planned future work will involve using our model to study whether the choice of fracturing criterion impacts the floe size distribution resulting from propagating wave-induced breakup, and applying it in configurations corresponding to recent and exciting observations of transient wave-induced breakup of instrumented ice in a natural setting (Kuchly et al., 2025).

Appendix A: Equation for the moment-deformation

We consider the boundary problem

$$\left\{\begin{array}{llr}\frac{{\mathrm{d}}^{\mathrm{4}}w}{\mathrm{d}{x}^{\mathrm{4}}}={\mathit{\sigma }}^{\mathrm{4}}\left[\mathit{\eta }\right(x)-w(x\left)\right]& x\in [\mathrm{0},L]& \left(\mathrm{A}\mathrm{1}\mathrm{a}\right)\\ \frac{{\mathrm{d}}^{\mathrm{2}}w}{\mathrm{d}{x}^{\mathrm{2}}}=\mathrm{0}& x\in \mathit{\{}\mathrm{0},L\mathit{\}}& \left(\mathrm{A}\mathrm{1}\mathrm{b}\right)\\ \frac{{\mathrm{d}}^{\mathrm{3}}w}{\mathrm{d}{x}^{\mathrm{3}}}=\mathrm{0}& x\in \mathit{\{}\mathrm{0},L\mathit{\}}& \left(\mathrm{A}\mathrm{1}\mathrm{c}\right).\end{array}\right.$$

where η denotes the surface undergoing wave forcing, w the floe deflection, and $\mathit{\sigma }=\frac{\mathrm{1}}{L_{\mathrm{D}}}$ is the reciprocate of the flexural length.

The surface η is typically the superposition of propagating, attenuated wave modes, so that

$$\begin{array}{}\text{(A2)}& \mathit{\eta }\left(x\right)=\sum _j{\mathit{\eta }}_j\left(x\right)\end{array}$$

and

$$\begin{array}{}\text{(A3)}& {\mathit{\eta }}_j\left(x\right)=a_j\exp (-{\mathit{\alpha }}_{j}x)\sin (k_{j}x+{\mathit{\varphi }}_j),\end{array}$$

with a wave amplitude, α wave attenuation per unit distance, k wave number, and ϕ wave phase at the left floe edge. The index j is used with respect to a discretised wave spectrum.

Finally, we define the elastic energy per unit cross-sectional area

$$\begin{array}{}\text{(A4)}& E={\displaystyle \frac{D}{\mathrm{2}h}}\underset{\mathrm{0}}{\overset{L}{\int }}({\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}w}{\mathrm{d}{x}^{\mathrm{2}}}}{)}^{\mathrm{2}}\mathrm{d}x.\end{array}$$

In the rest of this document, we will note

$$\begin{array}{}\text{(A5)}& \mathit{\kappa }\left(x\right)={\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}w}{\mathrm{d}{x}^{\mathrm{2}}}}\end{array}$$

the curvature of the floe.

The ODE in Eq. (A1a) is linear, and so are the boundary conditions fourth order ODE, but it is linear, and so are the boundary conditions in Eqs. (A1b), (A1c). Here, we consider the simplified case where D and h are constant, making Eq. (A1a) a constant-coefficients, linear ODE.

A1 Solution to the BVP on the floe deflection

The ODE in Eq. (A1a) is linear and non-homogeneous. Its general solution w is the superposition of an homogeneous solution w_h and a particular solution w_p.

A1.1 Homogeneous solution

The characteristic polynomial of the homogeneous ODE

$$\begin{array}{}\text{(A6)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{4}}w}{\mathrm{d}{x}^{\mathrm{4}}}}+{\mathit{\sigma }}^{\mathrm{4}}w\left(x\right)=\mathrm{0}\end{array}$$

associated with Eq. (A1a) is

$$\begin{array}{}\text{(A7)}& P\left(s\right)=s^{\mathrm{4}}+{\mathit{\sigma }}^{\mathrm{4}}.\end{array}$$

It has solutions

$$\begin{array}{}\text{(A8)}& {\displaystyle }(\mathrm{1}+i)\stackrel{\mathrm{̃}}{\mathit{\sigma }},(\mathrm{1}-i)\stackrel{\mathrm{̃}}{\mathit{\sigma }},(-\mathrm{1}+i)\stackrel{\mathrm{̃}}{\mathit{\sigma }},(-\mathrm{1}-i)\stackrel{\mathrm{̃}}{\mathit{\sigma }}\end{array}$$

with $\stackrel{\mathrm{̃}}{\mathit{\sigma }}=\frac{\sqrt{\mathrm{2}}}{\mathrm{2}}\mathit{\sigma }$. Independent solutions to Eq. (A1a) are thus

$$\begin{array}{}\text{(A9)}& \stackrel{\mathrm{ˇ}}{\mathit{f}}=\left[\begin{array}{c}\exp \left(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\right)\left(\cos \right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x)+i\sin (\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\left)\right)\\ \exp \left(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\right)\left(\cos \right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x)-i\sin (\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\left)\right)\\ \exp (-\stackrel{\mathrm{̃}}{\mathit{\sigma }}x)\left(\cos \right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x)+i\sin (\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\left)\right)\\ \exp (-\stackrel{\mathrm{̃}}{\mathit{\sigma }}x)\left(\cos \right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x)-i\sin (\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\left)\right)\end{array}\right].\end{array}$$

Applying the full-rank linear transformation

$$\begin{array}{}\text{(A10)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\begin{array}{cccc}\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ -i& i& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& -i& i\end{array}\right]\end{array}$$

to $\stackrel{\mathrm{ˇ}}{\mathit{f}}$ yields the real-valued independent solutions

$$\begin{array}{}\text{(A11)}& \mathit{f}=\left[\begin{array}{c}\exp \left(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\right)\cos \left(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\right)\\ \exp \left(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\right)\sin \left(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\right)\\ \exp (-\stackrel{\mathrm{̃}}{\mathit{\sigma }}x)\cos \left(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\right)\\ \exp (-\stackrel{\mathrm{̃}}{\mathit{\sigma }}x)\sin \left(\stackrel{\mathrm{̃}}{\mathit{\sigma }}x\right)\end{array}\right].\end{array}$$

Finally, any linear combination f^Tc, with the real-valued vector

$$\begin{array}{}\text{(A12)}& \mathit{c}=\left[\begin{array}{c}{c}_{\mathrm{1}},c_{\mathrm{2}},c_{\mathrm{3}},c_{\mathrm{4}}\end{array}\right]\end{array}$$

is a solution to Eq. (A6).

A1.2 Particular solutions

The non-homogeneous term in Eq. (A1a) can be written

$$\begin{array}{}\text{(A13)}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{4}}\mathit{\eta }\left(x\right)& {\displaystyle }={\mathit{\sigma }}^{\mathrm{4}}\sum _j{\mathit{\eta }}_j\left(x\right)\text{(A14)}& {\displaystyle }& {\displaystyle }={\mathit{\sigma }}^{\mathrm{4}}\sum _j\mathrm{Im}\left[{\widehat{a}}_j{e}^{\mathrm{i}{\widehat{k}}_{j}x}\right]\end{array}$$

with the complex amplitudes $\widehat{a}=a{e}^{\mathrm{i}\mathit{\varphi }}$ and the complex wavenumbers $\widehat{k}=k+\mathrm{i}\mathit{\alpha }$. Using the exponential response formula, and the characteristic polynomial Eq. (A7), we obtain particular solutions of the form

$$\begin{array}{}\text{(A15)}& {\displaystyle }{w}_{p,j}& {\displaystyle }={\mathit{\sigma }}^{\mathrm{4}}\mathrm{Im}\left[{\displaystyle \frac{{\widehat{a}}_j{e}^{\mathrm{i}{\widehat{k}}_{j}x}}{{{\widehat{k}}_j}^{\mathrm{4}}+{\mathit{\sigma }}^{\mathrm{4}}}}\right]\text{(A16)}& {\displaystyle }& {\displaystyle }=\mathrm{Im}\left[{\displaystyle \frac{{\widehat{a}}_j{e}^{\mathrm{i}{\widehat{k}}_{j}x}}{({\widehat{k}}_j/\mathit{\sigma }{)}^{\mathrm{4}}+\mathrm{1}}}\right]\end{array}$$

so that the particular solution to Eq. (A1a) can be written

$$\begin{array}{}\text{(A17)}& w_{\mathrm{p}}=\sum _j{w}_{p,j}.\end{array}$$

In what follows, we will write

$$\begin{array}{}\text{(A18)}& {\stackrel{\mathrm{̃}}{a}}_j={\displaystyle \frac{{\widehat{a}}_j}{({\widehat{k}}_j/\mathit{\sigma }{)}^{\mathrm{4}}+\mathrm{1}}}\end{array}$$

the complex amplitude of these solutions. As $\mathit{\sigma },k_j,{\mathit{\alpha }}_j>\mathrm{0}$, this amplitude exists only if

$$\begin{array}{}\text{(A19)}& k_j\ne \stackrel{\mathrm{̃}}{\mathit{\sigma }}\wedge {\mathit{\alpha }}_j\ne \stackrel{\mathrm{̃}}{\mathit{\sigma }}.\end{array}$$

A1.3 Coefficients of the homogeneous solution

The coefficients of c have to be determined to enforce the boundary conditions. Enforcing these four conditions leads to the system

$$\begin{array}{}\text{(A20)}& {\mathbf{D}}_{\mathrm{BC}}{\mathbf{M}}_{\mathrm{BC}}\mathit{c}=\mathit{r}\end{array}$$

with

$$\begin{array}{}\text{(A21)}& {\mathbf{D}}_{\mathrm{BC}}=\left[\begin{array}{cccc}\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{3}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{2}\sqrt{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{3}}\end{array}\right],\end{array}$$

$$\begin{array}{}\text{(A22)}& \begin{array}{rl}& {\mathbf{M}}_{\mathrm{BC}}=\\ & \left[\begin{array}{cccc}\mathrm{0}& \mathrm{1}& \mathrm{0}& -\mathrm{1}\\ -e^{\mathit{\beta }}\sin \left(\mathit{\beta }\right)& e^{\mathit{\beta }}\cos \left(\mathit{\beta }\right)& e^{-\mathit{\beta }}\sin \left(\mathit{\beta }\right)& -e^{-\mathit{\beta }}\cos \left(\mathit{\beta }\right)\\ -\mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{1}\\ -e^{\mathit{\beta }}\sin (\mathit{\beta }+\frac{\mathit{\pi }}{\mathrm{4}})& e^{\mathit{\beta }}\cos (\mathit{\beta }+\frac{\mathit{\pi }}{\mathrm{4}})& e^{-\mathit{\beta }}\cos (\mathit{\beta }+\frac{\mathit{\pi }}{\mathrm{4}})& e^{-\mathit{\beta }}\sin (\mathit{\beta }+\frac{\mathit{\pi }}{\mathrm{4}})\end{array}\right],\end{array}\end{array}$$

where $\mathit{\beta }=\stackrel{\mathrm{̃}}{\mathit{\sigma }}L$, and

$$\begin{array}{}\text{(A23)}& \mathit{r}=\mathrm{Im}\sum _j\left[\begin{array}{c}{\widehat{k}}_j^{\mathrm{2}}{\stackrel{\mathrm{̃}}{a}}_j\\ {{\widehat{k}}_j}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{a}}_j{e}^{\mathrm{i}{\widehat{k}}_{j}L}\\ \mathrm{i}{\widehat{k}}_j^{\mathrm{3}}{\stackrel{\mathrm{̃}}{a}}_j\\ \mathrm{i}{\widehat{k}}_j^{\mathrm{3}}{\stackrel{\mathrm{̃}}{a}}_j{e}^{\mathrm{i}{\widehat{k}}_{j}L}\end{array}\right].\end{array}$$

The third and first lines of the system give

$$\left\{\begin{array}{lcr}{c}_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}}(r_{\mathrm{1}}-\frac{r_{\mathrm{3}}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}})+c_{\mathrm{3}}+\mathrm{2}{c}_{\mathrm{4}}& & \left(\mathrm{A}\mathrm{24}\mathrm{b}\right)\\ c_{\mathrm{2}}=\frac{r_{\mathrm{1}}}{\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}}+c_{\mathrm{4}}& & \left(\mathrm{A}\mathrm{24}\mathrm{b}\right)\end{array}\right.$$

and the system in Eq. (A20) can be simplified to the more tractable

$$\begin{array}{}\text{(A25)}& {\mathbf{M}}_{\mathrm{II}}\left[\begin{array}{c}{c}_{\mathrm{3}}\\ c_{\mathrm{4}}\end{array}\right]=\left[\begin{array}{c}-r_{\mathrm{3}}\\ r_{\mathrm{4}}\end{array}\right]\end{array}$$

with

$$\begin{array}{}\text{(A26)}& \begin{array}{rl}& {\mathbf{M}}_{\mathrm{II}}=\mathrm{2}\\ & \left[\begin{array}{cc}\sin \left(\mathit{\beta }\right)\sinh \left(\mathit{\beta }\right)& e^{\mathit{\beta }}\sin \left(\mathit{\beta }\right)-\cos \left(\mathit{\beta }\right)\sinh \left(\mathit{\beta }\right)\\ -\frac{\sqrt{\mathrm{2}}}{\mathrm{2}}\left[\sin \right(\mathit{\beta }\left)\cosh \right(\mathit{\beta })+\cos (\mathit{\beta }\left)\sinh \right(\mathit{\beta }\left)\right]& -\frac{\sqrt{\mathrm{2}}}{\mathrm{2}}{e}^{\mathit{\beta }}-\sin (\mathit{\beta }+\frac{\mathit{\pi }}{\mathrm{4}})\sinh \left(\mathit{\beta }\right)\end{array}\right].\end{array}\end{array}$$

The determinant of M_II is

$$\begin{array}{}\text{(A27)}& {\displaystyle }\mathrm{\Delta }{\displaystyle }=\sqrt{\mathrm{2}}\left[\mathrm{2}{\sin }^{\mathrm{2}}\right(\mathit{\beta })-\cosh (\mathrm{2}\mathit{\beta })+\mathrm{1}]\end{array}$$

so that Δ<0 for β>0, and ${lim}_{\mathit{\beta }\to \mathrm{0}}\mathrm{\Delta }=\mathrm{0}$. Therefore, the system in Eq. (A20) admits a unique solution as long as L is non-zero, and L_D finite. Solving Eq. (A25 and substituting into Eq. (A24)) leads to the solution to Eq. (A20), that can be written

$$\begin{array}{}\text{(A28)}& \mathit{c}=\mathbf{M}\mathit{r}\end{array}$$

where the coefficients of the matrix M are

$$\begin{array}{ll}\text{(A29a)}& {\displaystyle }{M}_{\mathrm{11}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{Q}}{e}^{-\mathrm{2}\mathit{\beta }}\left[\sqrt{\mathrm{2}}\sin \right(\mathrm{2}\mathit{\beta }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}})-e^{-\mathrm{2}\mathit{\beta }}]{\displaystyle }{M}_{\mathrm{12}}& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{Q}}{e}^{-\mathit{\beta }}\left[\sqrt{\mathrm{2}}\cos \right(\mathit{\beta }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}})\\ \text{(A29b)}& {\displaystyle }& {\displaystyle }+e^{-\mathrm{2}\mathit{\beta }}(\mathrm{3}\sin \mathit{\beta }-\cos \mathit{\beta })]\text{(A29c)}& {\displaystyle }{M}_{\mathrm{13}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}Q}}{e}^{-\mathrm{2}\mathit{\beta }}\left[\sin \right(\mathrm{2}\mathit{\beta })-\mathrm{1}+e^{-\mathrm{2}\mathit{\beta }}]\text{(A29d)}& {\displaystyle }{M}_{\mathrm{14}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}Q}}{e}^{-\mathit{\beta }}[\cos \mathit{\beta }-e^{-\mathrm{2}\mathit{\beta }}(\mathrm{2}\sin \mathit{\beta }+\cos \mathit{\beta }\left)\right]\text{(A29e)}& {\displaystyle }{M}_{\mathrm{21}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{Q}}{e}^{-\mathrm{2}\mathit{\beta }}[-\sqrt{\mathrm{2}}\cos (\mathrm{2}\mathit{\beta }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}})+\mathrm{2}-e^{-\mathrm{2}\mathit{\beta }}]{\displaystyle }{M}_{\mathrm{22}}& {\displaystyle }={\displaystyle \frac{\sqrt{\mathrm{2}}}{Q}}{e}^{-\mathit{\beta }}[-\sin (\mathit{\beta }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}})\\ \text{(A29f)}& {\displaystyle }& {\displaystyle }+e^{-\mathrm{2}\mathit{\beta }}\cos (\mathit{\beta }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}})]\text{(A29g)}& {\displaystyle }{M}_{\mathrm{23}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}Q}}{e}^{-\mathrm{2}\mathit{\beta }}[\mathrm{1}-\cos (\mathrm{2}\mathit{\beta }\left)\right]\text{(A29h)}& {\displaystyle }{M}_{\mathrm{24}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}Q}}{e}^{-\mathit{\beta }}\sin \mathit{\beta }[\mathrm{1}-e^{-\mathrm{2}\mathit{\beta }}]\text{(A29i)}& {\displaystyle }{M}_{\mathrm{31}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{Q}}[-\mathrm{1}+\sqrt{\mathrm{2}}{e}^{-\mathrm{2}\mathit{\beta }}\cos (\mathrm{2}\mathit{\beta }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\left)\right]{\displaystyle }{M}_{\mathrm{32}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{Q}}{e}^{-\mathit{\beta }}[\mathrm{3}\sin \mathit{\beta }+\cos \mathit{\beta }\\ \text{(A29j)}& {\displaystyle }& {\displaystyle }-\sqrt{\mathrm{2}}{e}^{-\mathrm{2}\mathit{\beta }}\sin (\mathit{\beta }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}})]\text{(A29k)}& {\displaystyle }{M}_{\mathrm{33}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}Q}}[-\mathrm{1}+e^{-\mathrm{2}\mathit{\beta }}(\sin \left(\mathrm{2}\mathit{\beta }\right)+\mathrm{1}\left)\right]\text{(A29l)}& {\displaystyle }{M}_{\mathrm{34}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}Q}}{e}^{-\mathit{\beta }}[\cos \mathit{\beta }-\mathrm{2}\sin \mathit{\beta }-e^{-\mathrm{2}\mathit{\beta }}\cos \mathit{\beta }]\text{(A29m)}& {\displaystyle }{M}_{\mathrm{41}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{Q}}[\mathrm{1}+e^{-\mathrm{2}\mathit{\beta }}(\sqrt{\mathrm{2}}\sin (\mathrm{2}\mathit{\beta }+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}})-\mathrm{2}\left)\right]\text{(A29n)}& {\displaystyle }{M}_{\mathrm{42}}& {\displaystyle }=M_{\mathrm{22}}\text{(A29o)}& {\displaystyle }{M}_{\mathrm{43}}& {\displaystyle }=M_{\mathrm{23}}\text{(A29p)}& {\displaystyle }{M}_{\mathrm{44}}& {\displaystyle }=M_{\mathrm{24}}\end{array}$$

with

$$\begin{array}{}\text{(A30)}& Q=-\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}\left[\right(\mathrm{1}-e^{-\mathrm{2}\mathit{\beta }}{)}^{\mathrm{2}}+\mathrm{2}{e}^{-\mathrm{2}\mathit{\beta }}\left(\cos \right(\mathrm{2}\mathit{\beta })-\mathrm{1})].\end{array}$$

The coefficients of M are implicit functions of $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$, and we note that the coefficients M_1j are even-symmetric to the coefficients M_3j, and the coefficients M_2j are odd-symmetric to the coefficients M_4j. This reproduces the respective evenness and oddness of f₁ and f₃, and f₂ and f₄. The leading exponential terms for all the coefficients M_1j and M_2j ensure that the deflection does not diverge for large floes.

The homogeneous solution to Eq. (A1) is then

$$\begin{array}{}\text{(A31)}& w_{\mathrm{h}}\left(x\right)=\sum _j{c}_j{f}_j\left(x\right)\end{array}$$

with the coefficients of f given from Eq. (A28).

A1.4 Summary

The solution to the BVP Eq. (A1) is given by the sum

$$\begin{array}{}\text{(A32a)}& {\displaystyle }w\left(x\right)& {\displaystyle }=w_{\mathrm{h}}\left(x\right)+w_{\mathrm{p}}\left(x\right)\text{(A32b)}& {\displaystyle }& {\displaystyle }=\sum _{j=\mathrm{1}}^{\mathrm{4}}{c}_j{f}_j\left(x\right)+\sum _{j=\mathrm{1}}^{N_{\mathrm{f}}}\mathrm{Im}\left[{\stackrel{\mathrm{̃}}{a}}_j{e}^{\mathrm{i}{\widehat{k}}_{j}x}\right].\end{array}$$

The definitions of the functions f_j are given in Sect. A1.1, the definitions of the coefficients ${\stackrel{\mathrm{̃}}{a}}_j$ as well as the complex wavenumbers ${\widehat{k}}_j$ are given in Sect. A1.2, and the definitions of the coefficients c_j in Sect. A1.3. The integer N_f is the number of frequency bins used to discretise a wave spectrum. The deflection is entirely determined by the elastic length L_D, the floe length L, and N_f tuples of amplitude a, wavenumber k, attenuation number α, and phase ϕ. Assuming independence of these quantities, the solution is parametrised by 2+4N_f real numbers. All of these, at the exceptions of the phases taking values in $(-\mathit{\pi },\mathit{\pi }]$, are positive. They can be further constrained to physically realistic ranges.

A2 Elastic energy

A2.1 Introduction

The elastic energy of a bent floe is defined as

$$\begin{array}{}\text{(A33)}& E={\displaystyle \frac{D}{\mathrm{2}h}}\underset{\mathrm{0}}{\overset{L}{\int }}({\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}w}{\mathrm{d}{x}^{\mathrm{2}}}}{)}^{\mathrm{2}}\mathrm{d}x.\end{array}$$

We introduce the floe curvature $\mathit{\kappa }\left(x\right):=\frac{{\mathrm{d}}^{\mathrm{2}}w}{\mathrm{d}{x}^{\mathrm{2}}}$. From Eq. (A30), we have

$$\begin{array}{}\text{(A34)}& {\displaystyle }{\mathit{\kappa }}_{\mathrm{h}}& {\displaystyle }=\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}(-c_{\mathrm{1}}{f}_{\mathrm{2}}+c_{\mathrm{2}}{f}_{\mathrm{1}}+c_{\mathrm{3}}{f}_{\mathrm{4}}-c_{\mathrm{4}}{f}_{\mathrm{3}})\text{(A35)}& {\displaystyle }{\mathit{\kappa }}_{\mathrm{p}}& {\displaystyle }=-\sum _j\mathrm{Im}\left[{\widehat{k}}_j^{\mathrm{2}}{\stackrel{\mathrm{̃}}{a}}_j{e}^{\mathrm{i}\widehat{k}x}\right].\end{array}$$

Let us define $b_j:=\left|{\widehat{k}}_j^{\mathrm{2}}{\stackrel{\mathrm{̃}}{a}}_j\right|$, ${\mathit{\beta }}_j:=\mathrm{Ang}{\widehat{k}}_j^{\mathrm{2}}{\stackrel{\mathrm{̃}}{a}}_j$. We can then rewrite

$$\begin{array}{}\text{(A36)}& {\mathit{\kappa }}_{\mathrm{p}}=-\sum _j{b}_j{e}^{-{\mathit{\alpha }}_{j}x}\sin (k_{j}x+{\mathit{\beta }}_j).\end{array}$$

Finally, we introduce the quantities

$$\begin{array}{}\text{(A37)}& {\displaystyle }{E}_{\mathrm{h}}=\underset{\mathrm{0}}{\overset{L}{\int }}{{\mathit{\kappa }}_{\mathrm{h}}}^{\mathrm{2}}\mathrm{d}x,E_{\mathrm{p}}=\underset{\mathrm{0}}{\overset{L}{\int }}{{\mathit{\kappa }}_{\mathrm{p}}}^{\mathrm{2}}\mathrm{d}x,E_{\mathrm{q}}=\underset{\mathrm{0}}{\overset{L}{\int }}{\mathit{\kappa }}_{\mathrm{h}}{\mathit{\kappa }}_{\mathrm{p}}\mathrm{d}x,\end{array}$$

which are the contribution to the elastic energy of respectively the homogeneous part of the displacement, the inhomogeneous part of the displacement, and their quadratic interaction.

A2.2 Homogeneous contribution

We start by expanding E_h as

$$\begin{array}{}\text{(A38)}& \begin{array}{rl}{E}_{\mathrm{h}}& ={c_{\mathrm{1}}}^{\mathrm{2}}{I}_{\mathrm{2}}+{c_{\mathrm{2}}}^{\mathrm{2}}{I}_{\mathrm{1}}+{c_{\mathrm{3}}}^{\mathrm{2}}{I}_{\mathrm{4}}+{c_{\mathrm{4}}}^{\mathrm{2}}{I}_{\mathrm{3}}\\ & +\mathrm{2}(-c_{\mathrm{1}}{c}_{\mathrm{2}}{I}_{\mathrm{12}}-c_{\mathrm{1}}{c}_{\mathrm{3}}{I}_{\mathrm{24}}+c_{\mathrm{1}}{c}_{\mathrm{4}}{I}_{\mathrm{23}}+c_{\mathrm{2}}{c}_{\mathrm{3}}{I}_{\mathrm{14}}\\ & -c_{\mathrm{2}}{c}_{\mathrm{4}}{I}_{\mathrm{13}}-c_{\mathrm{3}}{c}_{\mathrm{4}}{I}_{\mathrm{34}}),\end{array}\end{array}$$

with

$$\begin{array}{}\text{(A39)}& I_j=\underset{\mathrm{0}}{\overset{L}{\int }}{f_j}^{\mathrm{2}}\mathrm{d}x,I_{jn}=\underset{\mathrm{0}}{\overset{L}{\int }}{f}_j{f}_n\mathrm{d}x.\end{array}$$

These integrals evaluate to

$$\begin{array}{}\text{(A40a)}& {\displaystyle }{I}_{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{e^{\mathrm{2}\mathit{\beta }}\left(\sqrt{\mathrm{2}}\sin \right(\mathrm{2}\mathit{\beta }+\frac{\mathit{\pi }}{\mathrm{4}})+\mathrm{2})-\mathrm{3}}{\mathrm{8}\stackrel{\mathrm{̃}}{\mathit{\sigma }}}}\text{(A40b)}& {\displaystyle }{I}_{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{e^{\mathrm{2}\mathit{\beta }}(-\sqrt{\mathrm{2}}\sin (\mathrm{2}\mathit{\beta }+\frac{\mathit{\pi }}{\mathrm{4}})+\mathrm{2})-\mathrm{1}}{\mathrm{8}\stackrel{\mathrm{̃}}{\mathit{\sigma }}}}\text{(A40c)}& {\displaystyle }{I}_{\mathrm{3}}& {\displaystyle }={\displaystyle \frac{-e^{-\mathrm{2}\mathit{\beta }}\left(\sqrt{\mathrm{2}}\cos \right(\mathrm{2}\mathit{\beta }+\frac{\mathit{\pi }}{\mathrm{4}})+\mathrm{2})+\mathrm{3}}{\mathrm{8}\stackrel{\mathrm{̃}}{\mathit{\sigma }}}}\text{(A40d)}& {\displaystyle }{I}_{\mathrm{4}}& {\displaystyle }={\displaystyle \frac{e^{-\mathrm{2}\mathit{\beta }}\left(\sqrt{\mathrm{2}}\cos \right(\mathrm{2}\mathit{\beta }+\frac{\mathit{\pi }}{\mathrm{4}})-\mathrm{2})+\mathrm{1}}{\mathrm{8}\stackrel{\mathrm{̃}}{\mathit{\sigma }}}}\text{(A40e)}& {\displaystyle }{I}_{\mathrm{12}}& {\displaystyle }={\displaystyle \frac{-\sqrt{\mathrm{2}}{e}^{\mathrm{2}\mathit{\beta }}\cos (\mathrm{2}\mathit{\beta }+\frac{\mathit{\pi }}{\mathrm{4}})+\mathrm{1}}{\mathrm{8}\stackrel{\mathrm{̃}}{\mathit{\sigma }}}}\text{(A40f)}& {\displaystyle }{I}_{\mathrm{13}}& {\displaystyle }={\displaystyle \frac{L}{\mathrm{2}}}+{\displaystyle \frac{\sin \left(\mathrm{2}\mathit{\beta }\right)}{\mathrm{4}\stackrel{\mathrm{̃}}{\mathit{\sigma }}}}\text{(A40g)}& {\displaystyle }{I}_{\mathrm{14}}& {\displaystyle }={\displaystyle \frac{{\sin }^{\mathrm{2}}\mathit{\beta }}{\mathrm{2}\stackrel{\mathrm{̃}}{\mathit{\sigma }}}}\text{(A40h)}& {\displaystyle }{I}_{\mathrm{23}}& {\displaystyle }={\displaystyle \frac{{\sin }^{\mathrm{2}}\mathit{\beta }}{\mathrm{2}\stackrel{\mathrm{̃}}{\mathit{\sigma }}}}\text{(A40i)}& {\displaystyle }{I}_{\mathrm{24}}& {\displaystyle }={\displaystyle \frac{L}{\mathrm{2}}}-{\displaystyle \frac{\sin \left(\mathrm{2}\mathit{\beta }\right)}{\mathrm{4}\stackrel{\mathrm{̃}}{\mathit{\sigma }}}}\text{(A40j)}& {\displaystyle }{I}_{\mathrm{34}}& {\displaystyle }={\displaystyle \frac{-\sqrt{\mathrm{2}}{e}^{-\mathrm{2}\mathit{\beta }}\sin (\mathrm{2}\mathit{\beta }+\frac{\mathit{\pi }}{\mathrm{4}})+\mathrm{1}}{\mathrm{8}\stackrel{\mathrm{̃}}{\mathit{\sigma }}}}.\end{array}$$

The products of c_j,c_n pairs simplify little, and can be evaluated numerically.

A2.3 Particular contribution

We can expand E_p as

$$\begin{array}{}\text{(A41)}& E_{\mathrm{p}}=\sum _{j=\mathrm{1}}^{N_{\mathrm{f}}}[{b_j}^{\mathrm{2}}{I}_j^p+\mathrm{2}\sum _{n=j+\mathrm{1}}^{N_{\mathrm{f}}}{b}_j{b}_n{I}_{j,n}]\end{array}$$

with

$$\begin{array}{ll}\text{(A42a)}& {\displaystyle }{I}_j^p& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{L}{\int }}{e}^{-\mathrm{2}{\mathit{\alpha }}_{j}x}[\mathrm{1}-\cos (\mathrm{2}{k}_{j}x+\mathrm{2}{\mathit{\beta }}_j\left)\right]\mathrm{d}x{\displaystyle }{I}_{j,n}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{L}{\int }}{e}^{-{\mathit{\alpha }}_{j,n}x}\left[\cos \right(k_{j,n}^{-}x+{\mathit{\beta }}_{j,n}^{-})\\ \text{(A42b)}& {\displaystyle }& {\displaystyle }-\cos (k_{j,n}^{+}x+{\mathit{\beta }}_{j,n}^{+})]\mathrm{d}x\end{array}$$

where

$$\begin{array}{}\text{(A43)}& {\mathit{\alpha }}_{j,n}:={\mathit{\alpha }}_j+{\mathit{\alpha }}_n,k_{j,n}^{\pm }:=k_j\pm k_n,{\mathit{\beta }}_{j,n}^{\pm }:={\mathit{\beta }}_j\pm {\mathit{\beta }}_n.\end{array}$$

Let us define $K_j:=\left|{\widehat{k}}_j\right|$, ${\mathit{\theta }}_j:=\mathrm{Ang}{\widehat{k}}_j$. Assuming α_j≠0, we can then evaluate

$$\begin{array}{}\text{(A44)}& \begin{array}{rl}& I_j^p={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}[{\displaystyle \frac{\mathrm{1}-e^{-\mathrm{2}{\mathit{\alpha }}_{j}L}}{{\mathit{\alpha }}_j}}\\ & +{\displaystyle \frac{\sin (\mathrm{2}{\mathit{\beta }}_j-{\mathit{\theta }}_j)-e^{-\mathrm{2}{\mathit{\alpha }}_{j}L}\sin (\mathrm{2}{k}_{j}L+\mathrm{2}{\mathit{\beta }}_j-{\mathit{\theta }}_j)}{K_j}}].\end{array}\end{array}$$

Similarly, we define $K_{j,n}^{\pm }:=|k_{j,n}^{\pm }+\mathrm{i}{\mathit{\alpha }}_{j,n}|$ and ${\mathit{\theta }}_{j,n}^{\pm }:=\mathrm{Ang}(k_{j,n}^{\pm }+\mathrm{i}{\mathit{\alpha }}_{j,n})$, which leads to

$$\begin{array}{}\text{(A45)}& \begin{array}{r}{I}_{j,n}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}[{\displaystyle \frac{\sin ({\mathit{\beta }}_{j,n}^{+}-{\mathit{\theta }}_{j,n}^{+})-e^{-{\mathit{\alpha }}_{j,n}L}\sin (k_{j,n}^{+}L+{\mathit{\beta }}_{j,n}^{+}-{\mathit{\theta }}_{j,n}^{+})}{K_{j,n}^{+}}}\\ -{\displaystyle \frac{\sin ({\mathit{\beta }}_{j,n}^{-}-{\mathit{\theta }}_{j,n}^{-})-e^{-{\mathit{\alpha }}_{j,n}L}\sin (k_{j,n}^{-}L+{\mathit{\beta }}_{j,n}^{-}-{\mathit{\theta }}_{j,n}^{-})}{K_{j,n}^{-}}}].\end{array}\end{array}$$

A2.4 Quadratic interaction contribution

We can expand E_q as

$$\begin{array}{ll}{\displaystyle }{E}_{\mathrm{q}}& {\displaystyle }=\underset{\mathrm{0}}{\overset{L}{\int }}-\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}(-c_{\mathrm{1}}{f}_{\mathrm{2}}+c_{\mathrm{2}}{f}_{\mathrm{1}}+c_{\mathrm{3}}{f}_{\mathrm{4}}-c_{\mathrm{4}}{f}_{\mathrm{3}})\\ \text{(A46)}& {\displaystyle }& {\displaystyle }\sum _j\mathrm{Im}\left[{\widehat{k}}_j^{\mathrm{2}}{\stackrel{\mathrm{̃}}{a}}_j{e}^{\mathrm{i}\widehat{k}x}\right]\text{(A47)}& {\displaystyle }& {\displaystyle }=-\mathrm{2}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}\sum _{j=\mathrm{1}}^{N_{\mathrm{f}}}{b}_j[-c_{\mathrm{1}}{I}_{\mathrm{2},j}^q+c_{\mathrm{2}}{I}_{\mathrm{1},j}^q+c_{\mathrm{3}}{I}_{\mathrm{4},j}^q-c_{\mathrm{4}}{I}_{\mathrm{3},j}^q]\end{array}$$

with

$$\begin{array}{}\text{(A48)}& I_{n,j}^p=\underset{\mathrm{0}}{\overset{L}{\int }}{f}_n\left(x\right)\exp (-{\mathit{\alpha }}_{j}x)\sin (k_{j}x+{\mathit{\beta }}_j)\mathrm{d}x.\end{array}$$

We define

$$\begin{array}{}\text{(A49a)}& {\displaystyle }{Q}_j^{++}& {\displaystyle }=(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+{\mathit{\alpha }}_j{)}^{\mathrm{2}}+(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+k_j{)}^{\mathrm{2}}\text{(A49b)}& {\displaystyle }{Q}_j^{+-}& {\displaystyle }=(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+{\mathit{\alpha }}_j{)}^{\mathrm{2}}+(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-k_j{)}^{\mathrm{2}}\text{(A49c)}& {\displaystyle }{Q}_j^{-+}& {\displaystyle }=(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-{\mathit{\alpha }}_j{)}^{\mathrm{2}}+(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+k_j{)}^{\mathrm{2}}\text{(A49d)}& {\displaystyle }{Q}_j^{--}& {\displaystyle }=(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-{\mathit{\alpha }}_j{)}^{\mathrm{2}}+(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-k_j{)}^{\mathrm{2}},\end{array}$$

noting $Q^{--}$ is non-zero under the same condition Eq. (A19) that $\stackrel{\mathrm{̃}}{a}$ exists.

$$\begin{array}{}\text{(A50a)}& \begin{array}{rl}{I}_{\mathrm{2},j}^p& =e^{-\mathit{\beta }}\mathit{\left\{}{K}_j\right[\sin ({\mathit{\theta }}_j-{\mathit{\beta }}_j)({\displaystyle \frac{\mathrm{1}}{Q^{-+}}}-{\displaystyle \frac{\mathrm{1}}{Q^{--}}})\\ & +e^{-{\mathit{\alpha }}_{j}L}[{\displaystyle \frac{\sin \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+k_j)L-{\mathit{\theta }}_j+{\mathit{\beta }}_j)}{Q^{-+}}}\\ & +{\displaystyle \frac{\sin \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-k_j)L+{\mathit{\theta }}_j-{\mathit{\beta }}_j)}{Q^{--}}}\left]\right]\\ & +\sqrt{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}[-{\displaystyle \frac{\sin ({\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{-+}}}-{\displaystyle \frac{\cos ({\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{--}}}\\ & +e^{-{\mathit{\alpha }}_{j}L}[{\displaystyle \frac{\sin \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+k_j)L+{\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{-+}}}\\ & -{\displaystyle \frac{\sin \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-k_j)L-{\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{--}}}\left]\right]\mathit{\}}\end{array}\end{array}$$

$$\begin{array}{}\text{(A50b)}& \begin{array}{rl}{I}_{\mathrm{1},j}^p& =e^{-\mathit{\beta }}\mathit{\left\{}{K}_j\right[\cos ({\mathit{\theta }}_j-{\mathit{\beta }}_j)({\displaystyle \frac{\mathrm{1}}{Q^{-+}}}+{\displaystyle \frac{\mathrm{1}}{Q^{--}}})\\ & +e^{-{\mathit{\alpha }}_{j}L}[-{\displaystyle \frac{\cos \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+k_j)L-{\mathit{\theta }}_j+{\mathit{\beta }}_j)}{Q^{-+}}}\\ & +{\displaystyle \frac{\cos \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-k_j)L+{\mathit{\theta }}_j-{\mathit{\beta }}_j)}{Q^{--}}}\left]\right]\\ & +\sqrt{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}[{\displaystyle \frac{\cos ({\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{-+}}}-{\displaystyle \frac{\sin ({\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{--}}}\\ & +e^{-{\mathit{\alpha }}_{j}L}[-{\displaystyle \frac{\cos \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+k_j)L+{\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{-+}}}\\ & +{\displaystyle \frac{\cos \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-k_j)L-{\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{--}}}\left]\right]\mathit{\}}\end{array}\end{array}$$

$$\begin{array}{}\text{(A50c)}& \begin{array}{rl}{I}_{\mathrm{4},j}^p& =e^{-\mathit{\beta }}\mathit{\left\{}{K}_j\right[-\sin ({\mathit{\theta }}_j-{\mathit{\beta }}_j)({\displaystyle \frac{\mathrm{1}}{Q^{++}}}-{\displaystyle \frac{\mathrm{1}}{Q^{+-}}})\\ & +e^{-{\mathit{\alpha }}_{j}L}[-{\displaystyle \frac{\sin \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+k_j)L-{\mathit{\theta }}_j+{\mathit{\beta }}_j)}{Q^{++}}}\\ & +{\displaystyle \frac{\sin \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-k_j)L+{\mathit{\theta }}_j-{\mathit{\beta }}_j)}{Q^{+-}}}\left]\right]\\ & +\sqrt{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}[-{\displaystyle \frac{\cos ({\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{++}}}-{\displaystyle \frac{\sin ({\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{+-}}}\\ & +e^{-{\mathit{\alpha }}_{j}L}[{\displaystyle \frac{\cos \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+k_j)L+{\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{++}}}\\ & -{\displaystyle \frac{\cos \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-k_j)L-{\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{+-}}}\left]\right]\mathit{\}}\end{array}\end{array}$$

$$\begin{array}{}\text{(A50d)}& \begin{array}{rl}{I}_{\mathrm{3},j}^p& =e^{-\mathit{\beta }}\mathit{\left\{}{K}_j\right[-\cos ({\mathit{\theta }}_j-{\mathit{\beta }}_j)({\displaystyle \frac{\mathrm{1}}{Q^{++}}}-{\displaystyle \frac{\mathrm{1}}{Q^{+-}}})\\ & +e^{-{\mathit{\alpha }}_{j}L}[{\displaystyle \frac{\cos \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+k_j)L-{\mathit{\theta }}_j+{\mathit{\beta }}_j)}{Q^{++}}}\\ & +{\displaystyle \frac{\cos \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-k_j)L+{\mathit{\theta }}_j-{\mathit{\beta }}_j)}{Q^{+-}}}\left]\right]\\ & +\sqrt{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\sigma }}[-{\displaystyle \frac{\sin ({\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{++}}}-{\displaystyle \frac{\cos ({\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{+-}}}\\ & +e^{-{\mathit{\alpha }}_{j}L}[{\displaystyle \frac{\sin \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}+k_j)L+{\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{++}}}\\ & -{\displaystyle \frac{\sin \left(\right(\stackrel{\mathrm{̃}}{\mathit{\sigma }}-k_j)L-{\mathit{\beta }}_j+\frac{\mathit{\pi }}{\mathrm{4}})}{Q^{+-}}}\left]\right]\mathit{\}}.\end{array}\end{array}$$

Appendix B: Comparison to another mechanical model

In this Section, we compare the solution to floe bending issued from the model described in this publication, SWIIFT, to results issued from a model that solves for wave scattering from ice floes, hereafter referred to as WISIB (Mokus and Montiel, 2022). The main difference is that in the latter case, floe deflection is derived from the interface between fluid and floe, assuming that the floe conforms exactly to the fluid; the solution is sought assuming harmonic forcing of the plate, and incorporates forward and backward travelling waves, while SWIIFT only represents forward travelling waves. A consequence is the possibility for WISIB to create constructive interference locally increasing the deformation of the plate or its curvature. Interactions between floe lengths and wavenumbers can locally lead to resonances.

In Sect. B1, we look at curvature envelopes. In Sect. B2, we look at potential elastic energy derived from these curvatures.

B1 Curvature envelopes

In Fig. B2, we compare curvature envelopes derived from SWIIFT or WISIB. We do so for various ice thicknesses, wave periods, and floe lengths. By curvature envelope, we mean the maximum curvature attainable at any location along the length of a floe; the actual curvature oscillates and reaches it at its positive antinodes. In the case of WISIB, curvature is the sum of forward and backward travelling modes, forward and backward damped modes, and evanescent modes. As a consequence, the envelope itself oscillates. When the wavelength is long enough compared to the floe, these oscillations disappear.

When the wavelength gets significantly longer than the floe (for example, T=10 and 12 s and h=25 cm, most floe lengths), the curvature envelopes are different near the edges: SWIIFT has damped terms which oscillates with spatial frequency $\frac{\mathrm{1}}{\sqrt{\mathrm{2}}{L}_{\mathrm{D}}}$, while WISIB has damped terms which are complex solutions to the dispersion relation Eq. (19). The former depends only on ice thickness, while the latter depends primarily on the wave period.

B2 Energy from spectral forcing

More than the curvature itself, what matters to our energy-based fracture parametrisation is the potential elastic energy of a deformed floe. To compare these energies between the two models, we calculate them from the curvatures derived from both mechanical models. We do so for a spectral forcing, corresponding to a Pierson–Moskowitz spectrum discretised onto 47 frequency bins, between 0.05 to 0.52 Hz, with a width of 1 Hz. We use Latin hypercube sampling to generate an ensemble (size 289) of ice thicknesses, significant wave heights (parametrising the spectrum), and floe lengths, from uniform distributions on respectively 25 to 100 cm, 2 to 4 m, and 50 to 400 m. The spectra, integrated on the discretised frequency axis, show a median relative error to the targets $\frac{H_{\mathrm{S}}^{\mathrm{2}}}{\mathrm{16}}$ of $\mathrm{3.3}\times {\mathrm{10}}^{-\mathrm{3}}$. To each frequency bin, we associate a phase randomly sampled from 0 to 2π, to build an incoherent wave field, from which floe curvature, and eventually, elastic energy, is derived. We show the results in Fig. B3.

The resulting energies show a dependency to each of the three chosen variables, highlighted by the regression lines on the top row. The energy derived from SWIIFT is generally higher than the energy derived from WISIB, and exhibit more spread. However, the trends are similar for the energies derived from both models, so that the ratio of the two does not show any dependency to either of the variables (the regression lines on bottom left panel of Fig. B3 are mostly horizontal). In the bottom right panel of Fig. B3, we show the distribution of these ratios. It is right-skewed, and has mean 1.02 (geometric mean 0.69).

In Fig. B1, we show the correlation between our three input variables and the resulting energies. Because the trends exhibited on the top panel of Fig. B3 are non-linear, we use the Spearman correlation coefficient, which quantifies the monotonicity of a relationship. As can be seen from Fig. B3, and particularly the regression lines, the energy computed from either model are very mildly negatively correlated with thickness, mildly positively correlated with significant wave height, and clearly positively correlated with floe length. The correlations are stronger when using WISIB, which produces less scattered results. The energy ratio, however, is at most very lightly correlated with any of the variables.

Figure B1 Correlation matrix between input parameters and energy derived from SWIIFT and WISIB. We use the Spearman correlation coefficient.

Figure B2 Comparison of curvature envelopes derived from SWIIFT or WISIB for different periods of wave forcing (row) and different ice thicknesses (columns). Within each panel, from top to bottom (and darker to lighter hue) are floe lengths of 50, 100, 200, and 400 m, and the solid line is the SWIIFT solution, while the dashed line is the WISIB solution. The x-axes are normalised with respect to floe lengths, and individual curvatures are normalised with respect to the maximal curvature computed with SWIIFT. We display only the positive branch of the envelopes, which are symmetrical with respect to each corresponding y-axis. The origin of these y-axes is shown with a thin horizontal black line.

Figure B3 Comparison of potential elastic energy derived from SWIIFT or from WISIB. Top row: energy as a function of three parameters, with coloured dots for SWIIFT and black dots for WISIB. Lowess regression lines are superimposed. Bottom row: on the right, ratios of energy, as described in the text; the colours correspond to the variables of the top row, which were normalised to fit on the same axis. On the right, distribution of the energy ratio.