Articles | Volume 16, issue 10
The Cryosphere, 16, 4447–4472, 2022
The Cryosphere, 16, 4447–4472, 2022
Research article
20 Oct 2022
Research article | 20 Oct 2022

Wave-triggered breakup in the marginal ice zone generates lognormal floe size distributions: a simulation study

Wave-triggered breakup in the marginal ice zone generates lognormal floe size distributions: a simulation study
Nicolas Guillaume Alexandre Mokus and Fabien Montiel Nicolas Guillaume Alexandre Mokus and Fabien Montiel
  • Department of Mathematics and Statistics, University of Otago, Dunedin, New Zealand

Correspondence: Nicolas Guillaume Alexandre Mokus (


Fragmentation of the sea ice cover by ocean waves is an important mechanism impacting ice evolution. Fractured ice is more sensitive to melt, leading to a local reduction in ice concentration, facilitating wave propagation. A positive feedback loop, accelerating sea ice retreat, is then introduced. Despite recent efforts to incorporate this process and the resulting floe size distribution (FSD) into the sea ice components of global climate models (GCMs), the physics governing ice breakup under wave action remains poorly understood and its parametrisation highly simplified. We propose a two-dimensional numerical model of wave-induced sea ice breakup to estimate the FSD resulting from repeated fracture events. This model, based on linear water wave theory and visco-elastic sea ice rheology, solves for the scattering of an incoming time-harmonic wave by the ice cover and derives the corresponding strain field. Fracture occurs when the strain exceeds an empirical threshold. The geometry is then updated for the next iteration of the breakup procedure. The resulting FSD is analysed for both monochromatic and polychromatic forcings. For the latter results, FSDs obtained for discrete frequencies are combined following a prescribed wave spectrum. We find that under realistic wave forcing, lognormal FSDs emerge consistently in a large variety of model configurations. Care is taken to evaluate the statistical significance of this finding. This result contrasts with the power law FSD behaviour often assumed by modellers. We discuss the properties of these modelled distributions with respect to the ice rheological properties and the forcing waves. The projected output can be used to improve empirical parametrisations used to couple sea ice and ocean wave GCM components.

1 Introduction

Sea ice is a distinctive feature of both polar oceans and has a profound influence on our climate. It blankets a significant fraction of the Earth, is hard to reach, and offers particularly harsh fieldwork conditions. Consequently, numerical modelling is a valuable tool not only for forecasting ice extent evolution, but also to gain insights, at a global scale, into the physical processes shaping this evolution. Hindcast results straying away from observations (Stroeve et al.2007) hint at not fully understood internal climate variability (Zhang et al.2018; Castruccio et al.2019) or missing physics, such as the effect of waves on the ice cover (Squire2020). Global climate models (GCMs) have typically overlooked this impact, even if advances were made in recent years (Roach et al.2018; Boutin et al.2020). Waves can break the ice, especially as thinner ice becomes prevalent. For instance, the 2012 record-low Arctic sea ice extent was amplified by wave activity (Parkinson and Comiso2013). With thinner and weaker first-year ice becoming dominant in the Arctic (Kwok et al.2009; Kwok2018), the sea ice grows more vulnerable to flexure-induced failure.

The marginal ice zone (hereafter MIZ), a belt of loosely to densely packed ice floes, serves as a buffer between the ice-free open ocean and the pack ice; it is a region notably affected by waves (Dumont et al.2011). Tracking these individual, floating pieces of ice and characterising their dynamics at the basin scale is not possible. Since the pioneering work of Rothrock and Thorndike (1984), researchers have taken interest in describing the floe size distribution (hereafter FSD) and its effect on the climate system. Of particular interest here, fragmentation caused by ocean waves makes the floes more sensitive to melt (Steele1992; Perovich and Jones2014), even for larger ones (Horvat et al.2016), locally decreasing the ice concentration and allowing waves to propagate further into the MIZ. It leads to more fragmentation, thus introducing an ice–wave feedback loop (Asplin et al.2012; Thomson and Rogers2014).

Remote sensing observations of floe sizes (e.g. Rothrock and Thorndike1984; Toyota et al.2006; Steer et al.2008; Toyota et al.2011; Wang et al.2016) have led to the rooted conception that the FSD follows a power law: an – often truncated – Pareto distribution. However, it is unclear how much waves contribute to the emergence of this power law, as wave conditions prior to or during observations, as well as ice properties, are not always reported (Herman et al.2021). A variety of other processes (such as failure from wind or internal stress, lateral melting or growth, ridging, rafting, or welding) are thought to be likely to alter the FSD (Rothrock and Thorndike1984). Additionally, a broad spectrum of acquisition techniques, areas and times of studies, and a priori assumptions (e.g. the parametric form of the distribution tail) led to parametrisations of this power law covering a large span of exponents. Stern et al. (2018) exposed that the widespread distribution fitting technique used, least squares regressions in log–log space, is likely to have led to significant bias in these exponent estimates.

For the last decades, wave–ice interaction has been an active field of research, and numerous theoretical models, gradually benefiting from advances in computing power, have been developed. Classically, the goal was to understand the attenuation imposed on the waves by the ice cover (Wadhams1973). A component of this attenuation is scattering, which results from inhomogeneities of the ice cover such as cracks and leads (Williams and Porter2009), changes in thickness as caused by pressure ridges (Vaughan and Squire2007), and floes themselves (Meylan and Squire1994; Kohout and Meylan2008; Bennetts and Squire2011; Montiel et al.2016). The thin plate model, which represents floes as floating elastic plates, has been commonly used to retrieve an attenuation rate, which has a functional dependence on frequency (e.g. Kohout and Meylan2008). Dissipative terms can be included in these formulations (Squire and Fox1992; Fox and Squire1994). Alternatively, rheological models seek to represent the ice cover as a viscous continuum whose properties are representative of a discontinuous field (Keller1998; Wang and Shen2010; Santi and Olla2017). They can be adapted to account for grease ice conditions, rather than for discrete floes, and can also include pancake ice, but it is unclear whether the benefits they bring outshine the introduced complexity (Mosig et al.2015). Concurrently, an extensive body of observational research (e.g. Squire and Moore1980; Wadhams et al.1988; Meylan et al.2014; Montiel et al.2018, 2022) has been conducted and used to study attenuation.

The reciprocal response of the ice to the waves is unsatisfactorily understood, as direct observations of wave-induced floe breakup are scarce and localised in both time and space. Numerical assessment of the feedbacks between breakup and wave propagation were pioneered by Dumont et al. (2011) and built upon by Williams et al. (2013). Various models have implemented a breakup parametrisation, either to investigate the FSD (Montiel and Squire2017; Herman2017) or to evaluate the impact of its introduction on other quantities such as ice thickness or concentration (Roach et al.2018; Bateson et al.2020; Boutin et al.2021). These parametrisations are usually based on a fracture criterion involving either stress (Williams et al.2017; Montiel and Squire2017) or strain (Kohout and Meylan2008; Williams et al.2013; Horvat and Tziperman2015; Boutin et al.2018) or a combination of both (Dumont et al.2011). When these quantities exceed a critical value, breakup is triggered. These models cover a large span of complexities, from ad hoc configuration with simplified geometry to inclusion in a global sea ice model run in stand-alone mode (Bennetts et al.2017; Williams et al.2017; Roach et al.2018; Bateson et al.2020) or coupled to other GCM components (Roach et al.2019; Boutin et al.2021). The inclusion of the FSD in sea ice models is therefore actively being developed.

Zhang et al. (2015) proposed an FSD theory treating breakup as a stochastic process redistributing floe sizes and abstracting the wave forcing into a model parameter depending on the local wind, ice concentration, floe size, and ice thickness. The model was calibrated to reproduce observations, assuming a power law distribution of floe sizes. Boutin et al. (2018) generalised the formulations developed by Dumont et al. (2011) and Williams et al. (2013) to adapt them to the spectral wave model WAVEWATCH III® (WW3). Wave attenuation is determined by a mean floe size, computed assuming a power law FSD. Their work was extended by coupling WW3 to the sea ice models LIM3 (Boutin et al.2020) and neXtSIM (Boutin et al.2021). The former assumed a power law FSD with fixed exponent, while the latter let the exponent vary based on ice conditions. Bateson et al. (2020) represent the FSD as a power law in CICE and allow the exponent to vary seasonally or under the influence of the ice concentration. Therefore, it appears that most representations of wave-induced sea ice breakup in large-scale wave or sea ice models are based on the assumption that the FSD follows a power law. The prognostic model of Roach et al. (2018), based on the FSD modelling framework proposed by Horvat and Tziperman (2015), does not make this assumption, at the expense of computational cost. Horvat and Roach (2022) were able to reduce the additional cost attributable to the FSD evolution by training a neural network with the results from the original parametrisation.

Recent numerical experiments have been conducted to investigate the wave effect on the FSD without a priori assumptions on the distribution shape. Montiel and Squire (2017) extended the 3D linear wave scattering model of wave attenuation in the MIZ proposed by Montiel et al. (2016) by including a stress-based failure criterion. They investigated the FSD obtained after repeated breakup events and found that near-normal or bimodal distributions emerged for a wide range of wave and ice conditions. However, computational constraints limited their ability to perform simulations on sufficiently large scales to conduct robust statistical analyses of these distributions. Herman (2017) coupled a non-hydrostatic, non-linear wave model to sea ice represented as bonded grains, hence relaxing assumptions inherent to potential flow theory and allowing for the computation of a transient solution, at the expense of computational efficiency, limiting the usability of the model to smaller-scale configurations. The resulting FSDs were narrow, bounded distributions governed by the grain sizes. Both approaches rely on some binning of the floe sizes, either directly (Montiel and Squire2017) or indirectly through the use of discrete elements (Herman2017), effectively ensuing discrete distributions. Recent field observations of floes directly impacted by waves (Dumas-Lefebvre and Dumont2021; Herman et al.2021) and laboratory experiments (Herman et al.2018; Dolatshah et al.2018; Passerotti et al.2022) also suggest contrasting distributions.

In this study, we model the wave-induced breakup process with the aim of quantifying the resulting FSD. Breakup happens on timescales shorter than other processes affecting the FSD, such as thermodynamics (Collins et al.2015), allowing us to neglect these in our model. We use the thin plate formalism in a two-dimensional setting (one horizontal, one vertical), relying on an established scattering theory, augmented with an energy-dissipating process. We then incorporate a strain-based breakup parametrisation. We let an FSD emerge by repeatedly breaking off floes from a semi-infinite ice cover, and we link the resulting distribution to the ice properties and the wave forcing. This approach is similar to the work presented by Montiel and Squire (2017); however, we simplify the geometry, stripping out one horizontal dimension, hence reducing the numerical costs to allow for the generation of more fragments. We observe that under a realistic wave forcing, our model generates FSDs appropriately described by lognormal distributions; this holds in a large span of model configurations. We discuss the effects of the wave and the ice properties on the distribution parameters. Even though we acknowledge that any parametric distribution is likely to be an inaccurate depiction of a real ice cover, they have the advantage of efficiently encoding the information to be exchanged between GCM components (Horvat and Tziperman2015).

This paper unwinds as follows. In Sect. 2, we introduce notations and the underlying mathematical formulations. In Sect. 3, we describe the main components of our model, scattering and breakup. In Sect. 4 we present direct results from monochromatic simulations, upon which we build in Sect. 5 to suggest a polychromatic parametrisation. We discuss these results in Sect. 6.

2 Preliminaries

We consider surface gravity waves propagating in a two-dimensional fluid domain of constant, finite depth H associated with a Cartesian coordinate system (x,z), where x and z are the horizontal and vertical coordinates, respectively. Translational invariance is assumed in the second horizontal direction. We assume the fluid to be inviscid and incompressible with density ρw. The flow is assumed to be irrotational so that the fluid velocity can be described by the gradient of a scalar potential Φ, which satisfies Laplace's equation:

(1) 2 Φ = 0 .

We place an array of Nf+1 non-overlapping ice floes, modelled as floating visco-elastic plates, in the domain; two adjacent floes are separated by open water. Their mechanical behaviour is determined by their density ρ, thickness h, flexural rigidity D=Yh3121-ν2 (where Y and ν are respectively Young's modulus and Poisson's ratio), and viscosity γ; their draught is d=ρρwh. Floe j, where j0,,Nf, is located in space by the horizontal coordinate of its left edge xj (ordered so that xj<xj+1) and its length Lj, with LNf being infinite. At rest, the fluid region covered by floe j is encompassed in the sub-domain Ωjf=[xj,xj+Lj]×[-H,-d]. The interface with the ice Ωjf is on z=-d. We denote Ωf=0NfΩjf and Ωf=0NfΩjf. The ice-free sub-domain left of the floe-covered sub-domain Ωjf is Ωjw so that at rest,

(2) Ω j w = [ x j - 1 + L j - 1 , x j ] × [ - H , 0 ] , j > 0 ( - , x 0 ] × [ - H , 0 ] , j = 0 .

The horizontal bottom boundary ∂ΩH is at z=-H, and the interface Ωjw between the atmosphere and Ωjw is at z=0 when the fluid is at rest. We define Ωw and ∂Ωw in the same way as Ωf and ∂Ωf. The whole fluid domain is Ω=ΩwΩf. Our notations and the geometry of the model are summarised in Fig. 1.

The system is forced by a monochromatic plane wave of angular frequency ω propagating in the positive x direction. The wave amplitude a is assumed to be small compared to the wavelength. The perturbed top boundary of the fluid is located at z=η(x,t), whether it is an interface with the atmosphere (in which case η≈0) or with an ice floe (in which case η-d). We further set ​​​​​​​Φ=Reϕ(x,z)exp(-iωt) with ϕ a time-independent, complex-valued function.

Figure 1Geometry of the model at rest. Wave forcing would alter the fluid boundaries around z=0 and z=-d.


We consider the seabed to be impervious, hence not allowing for normal flow, so that on ∂ΩH

(3) ϕ z = 0 .

The small amplitude forcing allows us to use linear surface wave theory in Ωw, leading to the boundary condition

(4) ϕ z = ω 2 g ϕ

on ∂Ωw, where g is the acceleration due to gravity. We model the flexural motion of the ice floes using the modified Kirchhoff–Love plate theory, introduced by Robinson and Palmer (1990). Coupling to the fluid motion in Ωf yields

(5) D ρ w g 4 x 4 + 1 - ω 2 d g - i γ ω ρ w g ϕ z = ω 2 g ϕ

on ∂Ωf. Equations (4) and (5) stem from the assumptions that under small-amplitude wave forcing, waves do not break, and fluid is at all times in contact with the bottom of the ice. Details on the derivations can be found in, for example, Fox and Squire (1994) or Williams et al. (2013).

We also neglect the surge motion of the floes, meaning that

(6) ϕ x = 0

on xj,xj+Lj×[-d,0].

We finally add the free-edge conditions

(7) 3 ϕ x 2 z = 0 , 4 ϕ x 3 z = 0

on xj,xj+Lj×-d, which assume that bending moment and vertical stress vanish at the floe boundaries, respectively.

The boundary conditions given in Eqs. (3)–(7) together with Eq. (1) in its time-independent form, 2ϕ=0, complete our boundary value problem.

3 Methods

The solution to the boundary value problem described in Sect. 2 is based on eigenfunction expansions and multiple wave scattering theory. A closely related problem (no viscous dissipation and zero-draught approximation) was considered by Kohout and Meylan (2008) and forms the basis of the wave attenuation formulation in Horvat and Tziperman (2015); it was revisited by Montiel et al. (2012) with non-zero draught. Although similar in some aspects, our method differs from these two studies in the ways scattering by a single floe edge is treated and multiple scattering is resolved.

3.1 Wave scattering

In any sub-domain Ωjw or Ωjf, the velocity potential is decomposed as the superposition of forward-travelling and backward-travelling plane waves, using the ansatz c+eikx-θ++c-e-ikx-θ-ζ(z), where c+,c-C are coefficients to be determined, and θ+,θ-[0,2π) are phase shifts introduced to simplify analytical derivations and improve numerical stability. Solving the boundary value problem described in Sect. 2 in all sub-domains, the potential is expanded into series of travelling and evanescent wave modes

(8) ϕ = n = 0 c j - 1 , n w + e i k n w x - θ j - 1 , n w + + c j , n w - e - i k n w x - θ j , n w - ζ n w ( z ) , ( x , z ) Ω j w n = - 2 c j , n f + e i k n f x - θ j , n f + + c j , n f - e - i k n f x - θ j , n f - ζ n f ( z ) , ( x , z ) Ω j f

where superscripts w and f are related to an open-water and floe-covered sub-domain.

The wavenumbers {knwn0} are the roots of the dispersion relation

(9) k tanh k H = ω 2 g

such that k0w is a positive real number (therefore associated with left- and right-propagating wave modes in Ωw), while {knwn>0} are purely imaginary numbers with a positive imaginary part, sorted by ascending imaginary part (associated with exponentially decaying evanescent wave modes).

Likewise, the wavenumbers {knfn-2} are the roots of the dispersion relation

(10) D ρ w g k 4 + 1 - ω 2 d g - i γ ω ρ w g k tanh k ( H - d ) = ω 2 g

such that {knfn>-2} are complex numbers in the first quadrant of the complex plane, sorted by ascending imaginary part for n>0, and k-2f is a complex number in the second quadrant of the complex plane.

Since ORek0fOImk0f, k0f is associated with left- and right-propagating wave modes in Ωf. In contrast, OReknfOImknf for n>0: these modes are associated with exponentially decaying wave modes. Finally, OReknf=OImknf for n{-2,-1}, so these two roots are associated with attenuating, propagating wave modes. Details on these behaviours can be found in Williams et al. (2013). In the special case where γ=0 (purely elastic floes), k0f is a positive real number, k-1f is in the first quadrant of the complex plane, k-2f=-k-1f, and {knfn>0} are purely imaginary numbers with a positive imaginary part. We note that when 1-ω2dg becomes negative (large frequencies or thickness), instead of having the three distinct roots k-2f,k-1f, and k0f, Eq. (10) may admit one double root or one triple root (Williams2006, p. 39) as the complex roots may become purely imaginary. Therefore, we enforce ωgd in this study to avoid this issue. For an ice thickness of 1 m, it corresponds to a minimum admissible period of approximately 1.90 s.

For any wave mode n, floe j radiates four waves: two from its left edge (with coefficients cj,nw-, and cj,nf+) and two from its right edge (with coefficients cj,nw+, and cj,nf-), except for the leftmost, semi-infinite floe (j=Nf), whose right edge is ignored. The coefficients c-1,nw+ and cNf,nf- are prescribed and represent the right-travelling incident wave forcing in Ω0w and the absence of forcing in ΩNff, respectively: only c-1,0w+=-igωa is non-zero. The quantity θ=θ-1,0w+ is an arbitrary phase associated with the forcing. We also note that {θ-1,nw+n>0} and {θNf,nf-n-2} do not take part in the computation and are left undefined. The remaining phases are determined so that the exponential terms in Eq. (8) reduce to 1 when evaluating ϕ at the edge radiating the wave; i.e.

(11) θ j , n w + = k n w x j + L j ; θ j , n w - = k n w x j ; θ j , n f + = k n f x j ; θ j , n f - = k n f x j + L j .

Finally, the functions

(12) ζ n w ( z ) = cosh ( k n w ( H + z ) ) cosh ( k n w H ) , z [ - H , 0 ] , ζ n f ( z ) = cosh ( k n f ( H + z ) ) cosh ( k n f ( H - d ) ) , z [ - H , - d ]

are vertical basis functions in the free-surface sub-domains and the ice-covered sub-domains, respectively.

3.1.1 Scattering by one floe edge

We obtain the solution to the multiple scattering problem of the incident wave by the ice floes by imposing continuity of pressure and normal velocity across the vertical boundaries between adjacent ice-free and ice-covered sub-domains, i.e. at each floe edge. These conditions are enforced by matching ϕ and u=ϕx on both sides of each interface. The single edge matching problem is solved using an integral equation method, as described by Williams and Porter (2009) and Mosig (2018).

Considering the scattering by the left edge of floe j and assuming knowledge of cm-1,nw+ and cm,nf-, the method generates a set of scattering relations relating these incident wave mode coefficients to those associated with wave modes propagating or decaying away from the edge, cj,nf+ and cj,nw-. When truncating the series in Eq. (8) to Nv evanescent modes, the relations can be summarised by the matrix equation

(13) c j f + c j w - = T j f w R j f R j w T j w f S j w 0 0 S j f c j - 1 w + c j f - ,

where TfwCNv+3×Nv+1, TwfCNv+1×Nv+3, RfCNv+3×Nv+3, and RwCNv+1×Nv+1 are matrices, respectively describing transmission and reflection of waves in either direction through the floe edge; cjw±=cj,0w±,,cj,Nvw±T, cjf±=cj,-2f±,,cj,Nvf±T are vectors of unknown coefficients; and SwCNv+1×Nv+1 and SfCNv+3×Nv+3 are diagonal phase shift matrices.

By symmetry, the scattering by the right edge of floe j is described by

(14) c j w + c j f - = T j w f R j w R j f T j f w S j f 0 0 S j + 1 w c j f + c j + 1 w - .

The reflection and transmission matrices depend only on the quantities present in the dispersion relations, Eqs. (9) and (10).

We assess convergence of the numerical procedure by investigating energy conservation after scattering for floes of zero viscosity. Our analysis proved Nv=2 to be adequate. Therefore, we use it in the rest of this study.

3.1.2 Scattering by an array of floes

Wave fields radiated by adjacent floes are coupled, which is clearly shown by Eqs. (13) and (14). To solve the multiple scattering problem described by these equations, we take advantage of the sparsity of the matrix representing the combined linear system, using a dedicated solver (Demmel et al.1999; Virtanen et al.2020). Specifically, we solve

(15) M c = f ,

where M is a tridiagonal block matrix and c the vector of unknown potential coefficients.

An array of Nf+1 finite floes leads to the matrix

(16) M ̃ = M 0 U 1 0 0 L 1 M 1 U 2 0 L 2 M 2 0 U N f 0 0 L N f M N f ,

where each block element of M̃ is a square matrix of size 4(Nv+2), and 0 denotes 0-filled matrices, with sizes compatible with other matrices in the same rows and columns.

As the last floe here is infinite, we obtain M introduced in Eq. (15) by trimming down M̃ from its last 2(Nv+2) rows and columns. While the size of M is 2Nv+22Nf+12, it has only 4Nv+2Nf2Nv+5+12 non-zero elements: this number grows linearly with Nf, instead of quadratically.

The vector of unknown coefficients is

(17) c = c 0 c N f - 1 c N f w - c N f f + T , with c j = c j w - c j f + c j f - c j w + T ,

and the forcing term

(18) f = e i θ R 0 w c - 1 w + T 0 w f c - 1 w + 0 0 T .

Building M and f and solving for c is linear in time for 𝒪(Nf)>10. For Nf=105, these operations are done in around 10 ms​​​​​​​ on an Intel Core i5-6300U 2016 laptop. Solving Eq. (15) for c fully determines the spatial part of the potential field in Ω.

3.2 Breakup parametrisation

To parametrise the breakup, we build upon the commonly used strain-based approach (e.g. Kohout and Meylan2008; Horvat and Tziperman2015). As opposed to these authors' approaches, we account for each floe individually rather than parametrising the strain decay from the number of floes or considering the ice cover to be continuous.

When using the plane stress approximation and our symmetry assumption, the strain ε̃j undergone by ice floe j is

(19) ε ̃ j x , z = - z 2 w x 2 ,

where x,z describes a coordinate system local to the floe, defined as x=x-xj and z=z-h2-d=z-h12-ρρ0, hence setting the origin on the intersection of the floe's left edge and horizontal middle surface. Under the plane stress approximation, the vertical displacement field undergone by the floe, w(x,t), does not depend on z: w(x,t)=η(x,t). As zh2, the maximum (in absolute value) strain, εj, is located on either surface of the floe, i.e. z=±h2. It follows that

(20) ε j x = h 2 2 η x 2 = h T 4 π Re i 2 x 2 ϕ z ,

with T=2πω the wave period.

Floe j is set to break when

(21) max x 0 , L j ε j x > ε c ,

where εc is an empirically determined strain threshold.

We situate the breakup point at

(22) x b = argmax ε j

so that a floe of length Lj is turned into two floes of length xb and Ljxb. Hence, the number of floes at most doubles, if all the floes break in a single breakup simulation.

3.3 Numerical experiment set-up

The values of the parameters kept fixed across all simulations are given in Table 1. Our experiments unwind as follows:

  1. The model is initialised with a set of physical parameters as input and a single, semi-infinite floe.

  2. The scattered wave field is determined and the strain field evaluated for every floe in the domain. All the floes for which the conditions for breakup are met are split, and the domain is updated.

  3. The second step is repeated until none of the floes break or a prescribed number of iterations reached.

These steps are summarised in Fig. 2. The output is the set of coordinate-length pairs describing the final geometry.

Figure 2Outline of the numerical experiment.


For a given iteration, all floes are scrutinised for breakup first, then their positions are updated. As we neglect floe motion, to prevent them from overlapping, they are assigned order-preserving random locations, determined through a process described in Appendix A. As we do not build any energy dissipation mechanism for the fluid in Ωw, the width of the gap between the floes impacts the phase of the wave, but not its amplitude as it reaches floe j+1 after transmission by floe j. Thus, whether floe j+1 breaks or not is unlikely (insomuch as its length allows sufficient strain to be reached) to be affected by this gap, but the breakup location may be. To account for this introduced randomness, we run each simulation as an ensemble with 50 separate realisations. As a side benefit, doing so removes the possible effect of local resonances in any single realisation, as described by Kohout and Meylan (2008). We ran a sensitivity analysis (see Appendix A) to confirm that the choice of the parameters governing the random placement does not affect the resulting distribution of floe lengths.

We note that in the arrangement considered here, the transmission and reflection matrices (Eqs. 13 and 14) are the same for every floe: they need only be computed once. This stems from the experiment design, as all floes are assumed to inherit their physical properties (i.e. thickness, density, rigidity, and viscosity) from the parent semi-infinite floe used to initialise the simulation. These properties could be altered after the breakup, but we choose to keep them constant as breakup happens on timescales shorter than those of the processes that would, for example, alter the thickness. Our model can however handle a general case with floes of varying thickness, flexural rigidities, densities, and viscosities. Such simulations are outside the scope of this paper. The information carried by the positions and dimensions of the floes are encapsulated in the phase shift matrices (Eqs. 13 and 14) which are truly floe-dependent.

The final result extracted from the simulation is the set of newly formed floe lengths, excluding the semi-infinite floe on the right of the domain, considered to be a steady-state FSD.

Table 1Fixed model parameters and their values.

Download Print Version | Download XLSX

4 Monochromatic forcing

We first investigate the FSD our model generates under monochromatic forcing with prescribed wave period T, corresponding to angular frequency ω=2πT. In addition to wave frequency, we seek to characterise the effect of the ice thickness h and the strain breakup threshold εc on the FSD.

Figure 3a and b show example histograms of FSDs obtained for T=8 s and a=50 cm, while Fig. 3c and d show the influence of T on the spread of the FSD. The strain thresholds considered span the range of observed values reported from scarce field measurements (Kohout and Meylan2008).

Figure 3Impact of varying T, εc, and h on the FSD with a=50 cm. (a) Histograms of the floe length for different εc. The bin width is 1 m and T=8 s. Lognormal fits are superimposed in dashed lines. (b) Same as (a) but for different h. (c) Evolution of the median floe size, when increasing εc, for different T. The shaded area indicates the corresponding interquartile range. (d) Same as (c) but for increasing h. In (a) and (c) h=1 m; in (b) and (d) εc=4×10-5. In (a) and (b), darker hues denote weaker ice (lower strain threshold, lower thickness), and lighter hues denote stronger ice. In (c) and (d), the contrasting dots indicate the values plotted in the corresponding top panels. All plotted quantities are ensemble averages over 100 realisations.


We obtain unimodal, right-skewed histograms (Fig. 3a, b), appropriately fitted by lognormal distributions (more details on this distribution are presented in Sect. 5.2). In the following, we use the median and the interquartile ranges as measures of central tendency and spread, respectively (Fig. 3c, d). We note that for right-skewed distribution, the median is bounded below by the mode of the distribution and bounded above by its mean. In this section, and in the following, we display so-called number distributions, such that the area of one histogram bin is the fraction of the obtained floes whose dimensions fall within the bin. Alternatively, considering so-called area distributions, where observations are weighted by their area (which we assume proportional to the length squared), alters the shape so that the peak of the distribution is shallower and shifted towards larger floes, and the tail of the distribution is thicker. The previous comments on unimodality, skewness, and lognormal fits, however, still hold.

Increasing εc has only a moderate effect on the shape of the FSD (the skewness being respectively 1.48±0.15, 1.46±0.23, and 1.63±0.37 for increasing values of εc displayed in Fig. 3a), a more noticeable effect being a shift towards larger floes. However, increasing h from 50 cm to 1 m has a more dramatic influence on the FSD, shifting it towards larger floes, increasing the spread, and thickening its tail. Doubling the thickness a second time brings the same sort of changes, in a less pronounced fashion. Therefore, an increased ice mechanical resistance (through either its thickness or its strain threshold) leads to the presence of larger floes. This is expected, as the free-edge boundary conditions (Eq. 7) impose zero strain at the floe edges, thus involving longer stretches of ice to reach higher absolute strains.

Figure 3c and d show two clear trends. Both longer waves and stronger ice lead to larger median floe sizes and wider interquartile ranges, with the exception of T=8 s in panel (c), which shows a decrease in spread with increasing εc. It shows that wave properties alone do not govern the dominant floe size. However, the trend in median floe size is more ambiguous at high thickness. The apparent oscillations observed at h>1.5 m, T=4 s are a spurious consequence of the fairly low amount of breakup observed in these configurations (Fig. 4b), enhancing the apparent variability. Similarly, at T=8 s, there is an order of magnitude drop in the number of floes when h increases from 1.5 to 2 m, which may cause the tapering off of and slight decrease observed in the median floe size. Notably, median floe sizes are much smaller than the corresponding half-wavelength (respectively 12, 50, and 112 m for periods of 4, 8, and 12 s), as observed by Herman et al. (2021).

Figure 4Various metrics at the end of a simulation. (a) Number of floes for different εc. (b) Same as (a) but for different h. (c) Minimum floe size in the sample for different εc. (d) Same as (c) but for different h. (e) Breakup width for different εc. (f) Same as (e) but for different h. The breakup width is defined as the cumulated length of the broken-off floes. In (a), (c), and (e) h=1 m; in (b), (d), and (f) εc=4×10-5; T=8 s, a=50 cm for all panels. The wavelength axis corresponds to the open-water propagating wavelength at a given T. Darker hues denote weaker ice (lower strain threshold, lower thickness), and lighter hues denote stronger ice. The plotted lines are ensemble average; the shaded areas indicate 1 standard deviation.


The final number of floes (number of floes reached when the forcing wave field no longer breaks any floe during a simulation) depends sharply on T, as shown in Fig. 4a and b. Three regimes can be identified: a rapid increase with T for lower periods (higher frequencies), then a plateau phase, preceding a sudden decrease at higher periods (lower frequencies). The precise delimitations of these regimes depend on the ice properties and can be explained by the non-linear relationship between T, h, and the undergone strain. Increasing h or T explicitly increases the maximum strain (Eq. 20). However, increasing T leads to a longer wavelength, translating to a decline in magnitude of the surface curvature term, offsetting the increase. Additionally, waves propagate with a longer wavelength under the ice than in open water: the thicker the ice, the longer the wavelength becomes. Therefore, increasing h also decreases the surface curvature. Lastly, the fraction of wave energy transmitted by a floe edge, close to 1 for longer waves, drops as the period decreases (Fox and Squire1990); the precise magnitude of this dip depends on floe length. When considering multiple scattering, these reflections exponentially stack up, making a few floes an effective barrier to low-period wave propagation.

The minimum floe size, shown in Fig. 4c and d, also follows three regimes roughly delimited by the same boundaries: a sharp decrease when increasing T from the lowest periods; then a steady growth (appearing to be evolving linearly with the wavelength) from a local minimum, corresponding to the plateau; and a more abrupt increase corresponding to the drop at higher periods. We observe this qualitative behaviour independently of the values of εc or h considered. An increase in minimum floe size with the wavelength is expected. The initial decrease, for shorter wavelength, corresponds to simulations with a very small number of floes: as previously mentioned, floes effectively reflect waves of low period. Short waves are not able to further fragment the first initial pieces that broke off the semi-infinite cover, which leads to these higher minima.

The breakup width, here defined as the cumulated length of ice broken off the semi-infinite floe and shown in Fig. 4e and f, is clearly concave down and peaks for wavelengths corresponding to the plateau phase. Interestingly, if a lower strain threshold leads to wider MIZs, so does thicker ice. If thick ice, when broken up, spawns less floes than thin ice and needs a higher wavelength to be bent enough to break, it allows for breakup over a much larger wavelength band and produces longer floes. Counterintuitively, the consequence is that intermediate waves are able to propagate much further into a thicker ice cover.

5 Polychromatic forcing

Our model, as described in Sects. 2 and 3, parametrises the wave forcing with a single amplitude–frequency pair. In this section, we propose an approach to extend our results to a polychromatic wave forcing.

5.1 Wave spectrum

To estimate the effect of a developed sea on the FSD fL(l), we take the weighted average of distributions f̃L(l,ω) resulting from monochromatic model runs, with weights taken from the energy density S(ω) of a theoretical ocean spectrum over a truncated frequency interval [ωmin,ωmax] so that

(23) f L ( l ) = ω min ω max f ̃ L ( l , ω ) S ( ω ) d ω ω min ω max S ( ω ) d ω .

This FSD model assumes that different frequency components of the wave spectrum affect the FSD independently from each other. This is a strong assumption, the validity of which is discussed in Sect. 6. Several expressions exist for S; we choose to use a Pierson–Moskowitz spectrum (Pierson and Moskowitz1964) as it can be easily adapted to depend only on the significant wave height Hs (Ochi2005), giving

(24) S ( ω ) = c 1 g 2 ω 5 exp - c 2 g 2 ω 4 H s 2 ,

where c1=8.1×10-3 and c2=3.24×10-2 are non-dimensional constants. This spectrum has been used in previous wave–sea-ice interaction studies (e.g. Kohout and Meylan2008; Dumont et al.2011) and is a reduced version of the two-parameter Bretschneider spectrum, used in this context as well (e.g. Horvat and Tziperman2015, 2017; Montiel and Squire2017).

We evaluate Eq. (23) numerically on 200 linearly spaced frequency bins, setting Hs=2a. As by definition

(25) H s = 4 0 + S ( ω ) d ω ,

the bounds of integration ωmin and ωmin are set so that the tails 16Hs20ωminS(ω)dω=16Hs2ωmax+S(ω)dω=5×10-7, which captures a significant part of the spectrum. If necessary, ωmax is adjusted to ensure ωgd, as discussed in Sect. 3.1.

For each of the 200 frequencies, we draw an FSD f̃L at random from the 50 realisations of this configuration. These FSDs are combined as previously mentioned. We repeat (with replacement) this random drawing stage 500 times in order to constitute a distribution ensemble from which statistics can be derived. Proceeding like so allows us to observe variations between different fL. As we have 50 independent realisations for each of the 200 frequencies, we have virtually infinitely many ways to build fL (502006×10339).

5.2 Reference configuration

We consider a reference configuration where Hs=1 m, h=1 m, and εc=4×10-5; this gives us frequency bounds that correspond to the wave period varying from 1.90 to 9.23 s. In Sect. 5.4, we discuss variations around this scenario.

Figure 5(a) FSDs for the reference configuration: lognormal fits to separate ensemble elements (thin black lines) and average distribution as detailed in text (thicker blue line). The underlying histogram depicts the ensemble average distribution. (b) Normal quantile–quantile plots of the standardised data: the coloured area corresponds to 1 standard deviation around the mean of the quantile–quantile lines; the dashed line indicates the main diagonal. (c) Cumulative FSDs: individual empirical CDFs (thin black lines) and CDF of the average distributions represented in (a) (thicker blue line).


The resulting FSD, shown in Fig. 5a, is remarkably well fitted by a three-parameter lognormal distribution. A random variable L is said to follow such a distribution with parameters μ, σ2, and τ if log (Lτ) is normally distributed with mean μ and variance σ2:

(26) L L N μ , σ 2 , τ log L - τ N μ , σ 2 .

The associated density function fL(l) is positive for l>τ: hence τ is the smallest floe size describable by this statistical model. The scale parameter s=exp μ has the physical dimension of the random variable – in this study, a length. The median and the mode of the distribution are given by

(27) median = τ + s ; mode = τ + s exp σ - 2 .

More details on the lognormal distribution can be found in Crow and Kunio (1988). In the following, we use the notation θ=s,σ,τ for the parameter vector.

We obtain a point estimate θ^ with maximum likelihood estimation (see e.g. Azzalini1996; Crow and Kunio1988). Among the ensemble, θ^ seems to be normally distributed, with strong correlations between the three parameters and small variances. Therefore, we use the mean vector θ to parametrise an underlying representative lognormal distribution, depicted in Fig. 5a and c as thick blue lines. The linear combination of lognormal distributions does not have a simplified expression and is not, generally, a lognormal. By defining the mean distribution as a lognormal parametrised with θ, we ensure this model is preserved. This can be justified by the low spread of θ^. In that regard, our random sampling aims at providing a confidence interval on the parameters values. The transformation between σ2 and σ and between μ and s is obviously not linear. However, for the range of values considered here, we find that averaging before or after taking the transformation leads to an absolute relative difference of less than 1 % (median for s: 4×10-2%; median for σ: 4×10-2%).

We outline the goodness of fit with a quantile–quantile plot shown in Fig. 5b. We standardised the data using Eq. (26) before deriving the quantiles, to ease the comparison between the 500 ensemble elements and use the symmetry property of the normal distribution. Therefore, an ideal match would have the data lying on the main diagonal. We observe departure from this line for larger floes, the shallower slope suggesting that the lognormal parametrisation overpredicts large floes that are not generated by the numerical model. However, as 68 % of the theoretical distribution belongs between the ±1 ticks of the horizontal axis, and 95 % belongs between the ±2 ticks, we deem the fit to be excellent. Figure 5c is another visualisation of the goodness of fit, comparing the (strictly speaking, complementary) empirical cumulative distribution functions (CDFs) to the CDF of the theoretical average lognormal distribution. They diverge significantly only from the 10−2 tick, indicating agreement for more than 99 % of the range of the data. As shown in Fig. 5c, the CDF of a random variable following a lognormal distribution could easily be misinterpreted as piecewise straight lines when represented on a log–log plot. This kind of graph is often used for field observations and is at the root of the ingrained power law or split power law conclusion.

5.3 Sub-domain FSD evolution

We do not expect the FSD to have the same shape all across an ice pack or even across the MIZ. To illustrate this effect, we analyse the evolution of the distribution when considering subsets of the domain.

Our experiment design generates distributions in parallel, for various periods in the spectrum. Therefore, there is no unique definition of the breakup width. To circumvent this issue, we use a sliding window whose bounds are relative to the local breakup width for each period used in the discretisation of the spectrum. We estimate the density f̃L;binfbsup, with 0binf<bsup1, not on all floe lengths Ljj0,,Nf-1 but on the subset

(28) L j b inf < Λ j Λ b sup ,

where Λj=m=0jLm is the cumulated floe size up to floe j, and Λ=ΛNf-1 is the total length of finite ice in the domain. We ignore the open-water gaps in the definition of the breakup width. These monochromatic densities are then combined as detailed in Sect. 5.1. The difference bsup-binf gives the width of the sliding window, which we fix at 0.5. The results of this procedure, applied to the reference configuration introduced in Sec. 5.2, are presented in Fig. 6a.

Figure 6(a) FSDs, fitted to lognormal distributions, in the whole domain and various sub-domains, as described in the text. The area under the partial FSDs is scaled by the number of observations, so that summing the first and the last windows, covering non-overlapping halves of the domain, yields the whole domain FSD. The effective number of floes, relative to the total number of floes in the domain, decreases steadily from 0.64 to 0.36. Densities averaged over 500 realisations. (b–d) Variation in the estimated parameters among the 500 realisations. In each panel, the white leftmost box corresponds to the whole domain; the coloured boxes follow the colour bar in panel (a).


The distributions remain remarkably lognormal, with the distribution parameters following regular trends for most positions of the sliding window (Fig. 6b–d). Our breakup parametrisation generates floes that tend to get smaller as they get further away from the semi-infinite floe marking the right boundary of the domain. This is true across all periods. As a consequence, the prevalence of larger floes grows for increasing binf values, shifting the distribution mode towards larger floes while thickening the distribution tail. This behaviour is similar to the effect of increasing the thickness, presented in Fig. 3b.

5.4 Forecast based on fitted parameters

We expand the analysis conducted in Sect. 5.2 to other combinations of Hs, h, and εc. For simplicity, we focus here on investigating the effect of varying one variable at a time from their reference values. Some alternative behaviour may emerge from multivariate simulations, which are outside the scope of this study. Hence, when the value of one variable is specified, the other variables assume their reference values, stated in Sect. 5.2. We assess the suitability of fitting the lognormal model though exploratory analysis, as in Sect. 5.2. Histograms of these simulations are presented in Appendix C.

We observe a remarkably good fit over most of the parametric space explored, with a few notable limitations. The smallest waves (significant wave heights between 40 and 60 cm) give rise to different patterns, which we do not analyse further. This is mostly due to them causing small amounts of breakup, leading to a limited number of floes. The empirical distributions obtained with thicker ice (h>1.4 m) are less skewed and have a more pronounced peak and thinner tails than the fitted lognormals. Across most configurations, some limitations arise in the tails, with fits for stronger ice overpredicting larger floes, while fits for higher waves underpredict them. We report the estimated parameters, averaged over 500 realisations for each configuration, in Fig. 7a–i. We note in Fig. 7g that τ^<0 for large enough wave height, which would allow the model to generate negative floe sizes. This issue will be discussed in the next section. Additionally, we compute the Kolmogorov–Smirnov statistic to further qualify the goodness of fit. The evolution of this statistic is presented in Appendix D.

Figure 7(a–i) Box-and-whisker plots of the estimated lognormal parameters for different model configurations. (j–l) Box-and-whisker plots of the modal floe size derived from the estimated parameters (Eq. 27). The plots show the variation across 500 realisations. The reference configuration is highlighted with a contrasting colour, where the size of the boxes allows.


As can be seen in Fig. 7, the lognormal fit parameters have fairly simple dependences on the physical variables and can be interpolated between computed values, with the exception of the low-amplitude outliers (Fig. 7d). We observe in Fig. 7k and l that the mode of the FSD (see Eq. 27), or modal floe size, grows with stronger ice, i.e. thicker floes or a larger strain threshold. This behaviour is analogous to the monochromatic case, as reported in Fig. 3. More surprisingly, the modal floe size first decreases with larger wave heights, reaching a local minimum for Hs=1.2 m before increasing with wave height. We believe this may emerge from the repartition of spectral energy between wave periods. Small significant wave heights do not induce much breakup, both because small-period waves, dominating the spectrum, are effectively reflected and because the smaller amplitudes cause lower strains. In contrast, higher significant wave heights have a larger component of high-period waves, leading to a smaller curvature of the ice floes, causing lower strains as well. This forecast mode has a strong dependence on the fitted shape parameter (Eq. 27). The sample mode is not captured accurately by the fit for Hs in the 1.4 to 2.2 m range, but the same behaviour is visible on the histograms (Fig. C1). This non-monotonic evolution does not support the notion that the wavelength alone governs the dominant floe size (as the peak propagating wavelength is, in the Pierson–Moskowitz framework, proportional to the significant wave height), as has been conjectured and observed (Dumas-Lefebvre and Dumont2021; Herman et al.2021).

The distribution parameters do not have a clear physical significance by themselves. Beyond the estimation of summary statistics, their main interest is the generation of floe sizes samples without the numerical cost of running the physical model. We illustrate such forecasts, and the associated errors, in Fig. 8. We use the mean distributions, parametrised by θ for each model realisation, to determine the ranges of floe sizes that would be predicted by the lognormal model. Every vertical slice in Fig. 8a–c is a representation of the predicted FSD for the relevant physical variable on the horizontal axis.

Figure 8(a–c) Ranges of forecast distributions, generated with mean estimated parameters. The lightest-colour area in each panel is the interquartile range; the dotted line denotes the median. The colour areas are symmetrical around the median. (d–f) Errors, with respect to the experimental quantiles, in the forecast quantiles. The represented quantiles bound the colour areas in panels (a)(c). The 25th, 50th, and 75th percentiles are respectively the first quartile, the median, and the third quartile. Negative values suggest overprediction, while positive values suggest underprediction. The reference configuration is highlighted by coloured ticks on the horizontal axes.


Figure 8b and c again show the dependence of the FSD on h, and εc leads to a behaviour similar to the monochromatic case. Increasing the ice thickness shifts the median floe size towards larger values and increases the spread, the distribution covering a larger span, especially for floes sizes beyond the 75th percentile. The smaller floe categories, below the 50th percentile, are shifted upwards as well and do not contribute to the increasing spread. Increasing the strain threshold sparks a steady increase in median floe size without much effect on the spread for εc3×10-5. The fraction of floes in the smallest size categories tends to increase. Increasing the significant wave height (hence, indirectly, the peak wavelength) leads to a more nuanced behaviour. The general trend is an increase in spread for both small and large floe extremes. There seems to be an inflection point around Hs=1.2 m for the growth of the 99.5th percentile, this wave height corresponding to the mode local minimum (see Fig. 7j).

Figure 8d–f point out the limitation of the lognormal model by showing the differences between numerical results and statistical predictions for chosen percentiles. As expected, differences arise for more extreme quantiles, corresponding to the long right tail of the distribution. These differences are more pronounced for small wave height and large thickness configurations. The prevalence of extreme floe sizes in the numerical results is, by definition, low. The errors in the three quartiles are smaller than 1 m on all the spans of the studied domains. The bulk of the FSD is hence well characterised by the lognormal model.

6 Discussion and conclusions

The emergence of a lognormal FSD from repeated wave-induced breakup is the key outcome of this paper. It contrasts with the power laws often assumed in modelling studies (Williams et al.2013; Bennetts et al.2017; Boutin et al.2020) or with more narrow distributions from process-based sea ice breakup modelling (Herman2017; Montiel and Squire2017). Contrary to the assumption made by Williams et al. (2013), we systematically obtain floes larger than half the wavelength (monochromatic simulations) or half the peak wavelength (polychromatic simulations). Most floes (about 95 %) are, however, smaller than that threshold. Anecdotally, the lognormal distribution has been reported for the size of brash ice pieces in navigation channels (Huang1988; Bonath et al.2019). One of its earliest applications was the description of particle sizes from repeated fragmentation events (Kolmogoroff1941). In a companion paper, we attempt to make a connection between the repeated fragmentation theory of Kolmogorov and the breakup of a sea ice cover caused by a wave event (Montiel and Mokus2022b).

The lognormal model does come with limitations. Field observations show the extensive spatial variability in the FSD. For instance, Paget et al. (2001) and Inoue (2004) both report an increase in the relative number of small floes when going towards the ice edge, respectively in the Antarctic and the Arctic. Our modelling results mirror this trend, highlighting the difficulty to settle on an all-around FSD parametrisation. We purposely ignore thermal and internal stress effects to focus on the effect of waves on the FSD, so validation with observational data would only be appropriate for a MIZ post-wave-induced breakup. As detailed in Sect. 5.4 and illustrated in Fig. 8d–f, the distribution struggles to capture the behaviour of the most extreme floe sizes in our simulations. We note, however, that for such floe sizes wave-induced breakup is not likely to be the dominant mechanism governing the evolution of the FSD (Roach et al.2018). We further found that the lognormal fit is not valid for small waves. The range of thicknesses we analysed is at the limit of validity, suggesting that even in this simplified setting the lognormal is not universal. Again, in such regimes, thermal and internal stress effects are likely to dominate over wave-induced breakup. Another point of concern is that we have τ^<0 for Hs≥2 m, meaning the probability of sampling negatively sized floe is not 0. Constraining the maximum likelihood estimation to yield positive estimates led to poor fit performance. This issue, which concerns a small fraction (<10-5) of the floes, can be easily circumvented by artificially bounding and rescaling the density function.

Our polychromatic forcing simulations (Sect. 5) lie on the underlying assumption that waves of different periods act independently on the ice. Although it is coherent with linear wave theory, which does not resolve non-linear interactions between wave components, it is unclear whether the resulting FSD describes reality appropriately. Other polychromatic parametrisations can be used to simulate the repeated breakup. For a given geometry, we could also compute the strain field obtained by linear superposition of strain fields at individual wave periods (with incident amplitude sampled from the prescribed spectrum and random phase). We compared the two approaches in the conference proceedings of Mokus and Montiel (2022a) and showed that they typically yield different results, with the strain superposition method generally leading to less breakup and larger floes. Qualitatively, however, the distributions follow similar shapes. It is not clear which one of these two approaches is physically more justifiable, as the strain-averaged approach assumes steady state to be reached by waves of all periods at the same time, even though dispersive effects tend to separate in time different spectral components of the sea state. Alternatively, the period–amplitude pairs can be sampled randomly between two model iterations (Montiel and Squire2017). Ultimately, the breakup of an ice cover in response to a wave event is a transient process. Consequently, validations of spectrum-generated FSDs will need to be sought against time-dependent simulations or experiments that control the wave forcing. Recent work by Dumas-Lefebvre and Dumont (2021) and Passerotti et al. (2022) may provide the necessary datasets to conduct such validation studies.

The Pierson–Moskowitz spectrum chosen in our polychromatic simulations was selected for its simplicity. To make sure our results are not qualitatively sensitive to the choice of spectrum, we conducted additional simulations of FSD generation for a range of different spectra, including symmetric ones. The simulations are described, and results are discussed in Appendix B. In short, we find the emergence of right-skewed FSDs is consistent across the spectra considered, suggesting that the lognormal FSD model is not an artefact of the choice of forcing spectrum.

Additionally to the wave height, the ice thickness, and the strain threshold, we analysed the effect of varying the ice viscosity γ (not shown here). We observed a significant contrast between γ=0 (purely elastic ice) and γ>0. The simulations we run with purely elastic floes are the only ones that reached the maximum number of iterations, set at 1000. It seems that multiple scattering alone is not effective enough at attenuating wave energy, leading to a rapid and sustained growth of the number of floes. However, marginal differences exist between ensuing distributions as long as some viscosity is introduced (γ in 1–100 Pa s m−1). Williams et al. (2013) used the same dissipation scheme with γ=13Pasm-1 derived from a 1979 campaign in the Bering Sea, while Mosig et al. (2015) fitted γ=6.9Pasm-1 to a 2012 Southern Ocean dataset (Kohout and Williams2015). Massom et al. (2018) derived γ=13.5Pasm-1 from the same dataset. We used γ=20Pasm-1 throughout, which is a bit more conservative but, as stated, does not significantly impact the results. Although associated with the ice, this parameter can be thought of as a parametrisation of the collection of all wave dissipation effects (Squire and Montiel2016).

A framework to model the evolution of the ice thickness distribution (ITD) has been introduced by Thorndike et al. (1975). The ITD does not have a preferred functional form (Dupont et al.2022) and is usually represented at the sub-grid level in sea ice models, such as CICE (Hunke et al.2021), by various thickness categories. Horvat and Tziperman (2015) extended this framework to include the floe size through a joint floe size and thickness distribution, which evolves under the action of separate physical processes such as thermodynamics, ridging, and wave-induced breakup. Therefore, the authors did not have to make any assumption on the shape of the FSD. This method was adapted and ported to CICE (Roach et al.2018), then coupled to a wave model (Roach et al.2019). The model setup of Roach et al. (2019) was then used to train a neural network model in order to accelerate the source terms' estimation in the FSD evolution model presented in Horvat and Roach (2022), which hence aims to replicate it. Their formulation for wave-induced breakup relies on generating a discrete probability density of floe sizes, dependent on the wave field (controlling surface elevation) and floe thickness (controlling wave attenuation, when not handled by a wave model), by deriving strain from the attenuated wave field and using it to populate a histogram of floe sizes. However, any other way to generate that density may be substituted, including parametric forms. If more evidence were to point towards the lognormal being a compelling choice for FSDs observed in nature, with parameters that could be linked to combinations of wave (significant wave height, peak period or wavelength) and ice (thickness) properties, as our study suggests, it would be a straightforward candidate. The (truncated) power law is another obvious candidate for such a substitution, and hindcasts comparing the two approaches could be a way to weigh in favour of one or the other.

Clauset et al. (2009) analysed 14 empirical datasets of different continuous variables, originally modelled with power laws, coming from a mixture of research areas. They observe that the power law is statistically appropriate for 8 of these datasets, while the lognormal holds for 13 of them. They use a relative goodness-of-fit test to show that the lognormal is more suitable than the power law in 12 cases, and significantly so in 4 cases, concluding the following: “In general, we find that it is extremely difficult to tell the difference between log-normal and power-law behaviour” (Clauset et al.2009). Stern et al. (2018) recommended a procedure for analysing floe size data in order to raise awareness of better fit methods, considering an alternative distribution as an optional step. We believe this study heads in that direction and that the lognormal distribution should be considered to be a viable alternative. Revisiting some studies tabulated by Stern et al. (2018) with this hypothesis could shed some light on its validity. Such results are presented in Montiel and Mokus (2022b).

Parametric distributions may never be flawless descriptions of the quantities they model. In this study, we show the relevance of the lognormal distribution when considering wave-induced breakup. We describe the evolution of the distribution shape for a range of ice properties, under various wave forcings. These results aim at being a step towards the parametrisation of wave action in FSD-evolving models.

Appendix A: Post-breakup floe positioning

This section details the process of redistributing the floes after a breakup event. The length of the gap between two floes does not matter inasmuch as the fluid is inviscid, and no wave energy is lost to it. The two main constraints are preserving the order of the floes and ensuring they do not overlap.

Floes are positioned from left to right and localised by their left edge, whose location is drawn from a uniform distribution. The leftmost floe is placed in an interval of fixed width, symmetric around its location at the previous iteration. For subsequent floes, the left bound corresponds to the right edge of the last positioned floe (on their left). The right bound corresponds to the previous right bound, augmented by the length of the last positioned floe. Locations are drawn between these two bounds; an illustration is given in Fig. A1. As the width of that interval quickly tends to 0, we enforce a minimal length (δmin) for our random draw. It is set to 1 cm throughout the paper. It does not mean that floes cannot get arbitrarily close to one another but that, if so needed, the right bound is moved further right to enforce this minimum width. The room allocated to the first floe (δinit) is set to 100 m throughout the paper.

We ran simulations with alternative values to ensure these values do not have any impact on our results. We used the case presented in Sect. 4, with T=8 s, a=50 cm, h=1 m, and γ=20Pasm-1. Summary statistics and CDFs are presented in Table A1 and Fig. A2.

Table A1Distribution sensitivity to various parameters of random positioning. All quantities, except the sample size, are expressed in metres. The statistics are averages computed over 50 realisations. The support is truncated to the 0.5th and 99.5th percentiles.

Download Print Version | Download XLSX

Figure A1Illustration of the floe repositioning method. The top row shows current floes with identified breakup location marked by vertical bars and the resulting lengths. Successive rows show the iterative positioning as described in the text. Below each row, a segment shows the interval from which a location will be randomly drawn for the next floe; the cross marks that location.


Figure A2Complementary CDFs of floe lengths, averaged over 50 model realisations, with varied positioning parameters and model parameters mentioned in the text.


Appendix B: Sensitivity to the choice of spectrum

We ran comparisons using alternative weighting functions, represented in Fig. B1. Descriptions of these functions can be found in Table B1.

The lognormal density function has, qualitatively, a shape similar to the Pierson–Moskowitz spectrum expressed as a function of frequency. However, we obtain similar, skewed unimodal densities with a range of symmetrical weighting functions, as displayed in Fig. B2.

The one case that stands out, with a lot of secondary peaks, is n_md. It corresponds to a function giving more weights to high-frequency waves (5 to 15 m), which is unrealistic in the context of our model.

Table B1Details of several weighting methods.

Download Print Version | Download XLSX

Figure B1Representation of functions used as alternative weights. They are all normalised on the positive real half-line; because of the finite range of frequencies supported by our model, truncations imply that some of them integrate to less than unity. Description of the legend entries can be found in Table B1.


Figure B2Densities obtained when combining monochromatic model runs with different weighting functions, as represented in Fig. B1.


Appendix C: Histograms

Figure C1Histograms and average lognormal fits, as described in Sect. 5.2, for varying significant wave height.


Figure C2Histograms and average lognormal fits, as described in Sect. 5.2, for varying ice thickness.


Figure C3Histograms and average lognormal fits, as described in Sect. 5.2, for varying strain threshold.


Appendix D: Kolmogorov–Smirnov statistic

The Kolmogorov–Smirnov statistic DKS is the largest difference between an empirical cumulative distribution function (CDF) and a reference CDF (Massey1951). By definition of the CDF, DKS is bounded by 0 and 1. When the distribution parameters have been estimated from the data, DKS can be used to run a Lilliefors test (Lilliefors1967). By comparing DKS to a critical value, depending on sample size and a chosen confidence level, one uses this test to reject a distribution hypothesis – not to confirm it. However, the power of this test, and others, notoriously increases with the sample size, making them able to detect trivial deviations from a reference distribution. This is a simple consequence of the fact that a model cannot perfectly fit the data. Hence, these tests only give a binary answer, not taking into account the usefulness of an imperfect model. A more purposeful alternative consists of studying the relative goodness of fit between different models, which we do not explore in detail here.

Instead of rejecting the lognormal hypothesis at an arbitrarily chosen confidence level, we report DKS as an indicative performance metric to compare our different configurations. More specifically, for each fitted lognormal with estimated parameters θ^ (that is, 500 different realisations per model configuration), we generate a random sample of size Neff from the distribution. We use Kish's effective sample size (Kish1965), the rounded-up ratio of the squared sum of weights to the sum of the squared weights, as Neff. We use maximum likelihood estimation to estimate θ^b from this sample, and we compute DKS for this sample and the distribution parametrised by θ^b. We then define ΔDKS as the difference between DKS from the random sample and DKS from our data. We repeat these three steps 1000 times to derive the distribution of ΔKS for each model configuration. This is analogous to the bootstrapping method described by Clauset et al. (2009). It follows that ΔDKS is bounded by −1 and 1, with ΔDKS>0 indicating cases where our data are fitted by the lognormal model better than data actually lognormally sampled, in terms of distance between the CDFs.

We report the results in Figs. D1D3. This procedure quantitatively illustrates the conclusions derived from analysing histograms and quantile–quantile plots.

Figure D1Distribution of ΔDKS, as defined in the text, for varying significant wave height. The boxes are bounded by the first and third quartiles, and the black lines are medians. The whiskers' lengths are 1.5 times the interquartile range. Black circles represent outliers.


Figure D2Same as Fig. D1 for varying ice thickness.


Figure D3Same as Fig. D1 for varying strain threshold.


Code and data availability

The model code, software tools developed for analysis, and the resulting preprocessed output presented in this paper are publicly available at (Mokus and Montiel2021)​​​​​​​. The raw output is available upon request.

Author contributions

NGAM and FM designed the numerical experiments. NGAM developed the model code, ran the simulations, and conducted the analysis. NGAM prepared the manuscript with significant input from and under the supervision of FM.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors thank Christopher Horvat for proposing the method of combining FSDs obtained at individual frequencies in order to form the FSD resulting from the breakup by a wave spectrum (Eq. 23).

Financial support

This research has been supported by the University of Otago (doctoral scholarship), the Marsden Fund (grant nos. 18-UOO-216 and 20-UOO-173), and Antarctica New Zealand (Antarctic Science Platform Project 4).

Review statement

This paper was edited by Yevgeny Aksenov and reviewed by three anonymous referees.


Asplin, M. G., Galley, R., Barber, D. G., and Prinsenberg, S.: Fracture of summer perennial sea ice by ocean swell as a result of Arctic storms, J. Geophys. Res.-Oceans, 117, C06025​​​​​​​,, 2012. a

Azzalini, A.: Statistical inference: based on the likelihood, in: Monographs on statistics and applied probability, 1st edn., 68, Chapman & Hall/CRC, Boca Raton, New York, ISBN 9780412606502, 1996. a

Bateson, A. W., Feltham, D. L., Schröder, D., Hosekova, L., Ridley, J. K., and Aksenov, Y.: Impact of sea ice floe size distribution on seasonal fragmentation and melt of Arctic sea ice, The Cryosphere, 14, 403–428,, 2020. a, b, c

Bennetts, L. G. and Squire, V. A.: On the calculation of an attenuation coefficient for transects of ice-covered ocean, P. Roy. Soc. A-Math. Phy., 468, 136–162,, 2011. a

Bennetts, L. G., O'Farrell, S., and Uotila, P.: Brief communication: Impacts of ocean-wave-induced breakup of Antarctic sea ice via thermodynamics in a stand-alone version of the CICE sea-ice model, The Cryosphere, 11, 1035–1040,, 2017. a, b

Bonath, V., Zhaka, V., and Sand, B.: Field measurements on the behavior of brash ice, in: Proceedings of the 25th International Conference on Port and Ocean Engineering under Arctic Conditions, Delft, The Netherlands, 9–13 June 2019, ISSN 0376-6756, (last access: 12 September 2022​​​​​​​), 2019. a

Boutin, G., Ardhuin, F., Dumont, D., Sévigny, C., Girard-Ardhuin, F., and Accensi, M.: Floe Size Effect on Wave-Ice Interactions: Possible Effects, Implementation in Wave Model, and Evaluation, J. Geophys. Res.-Oceans, 123, 4779–4805,, 2018. a, b

Boutin, G., Lique, C., Ardhuin, F., Rousset, C., Talandier, C., Accensi, M., and Girard-Ardhuin, F.: Towards a coupled model to investigate wave–sea ice interactions in the Arctic marginal ice zone, The Cryosphere, 14, 709–735,, 2020. a, b, c

Boutin, G., Williams, T., Rampal, P., Olason, E., and Lique, C.: Wave–sea-ice interactions in a brittle rheological framework, The Cryosphere, 15, 431–457,, 2021. a, b, c

Castruccio, F. S., Ruprich-Robert, Y., Yeager, S. G., Danabasoglu, G., Msadek, R., and Delworth, T. L.: Modulation of Arctic Sea Ice Loss by Atmospheric Teleconnections from Atlantic Multidecadal Variability, J. Climate, 32, 1419–1441,, 2019. a

Clauset, A., Shalizi, C. R., and Newman, M. E. J.: Power-Law Distributions in Empirical Data, SIAM Review, 51, 661–703,, 2009. a, b, c

Collins, C. O., Rogers, W. E., Marchenko, A., and Babanin, A. V.: In situ measurements of an energetic wave event in the Arctic marginal ice zone, Geophys. Res. Lett., 42, 1863–1870,, 2015. a

Crow, E. L. and Kunio, S.: Lognormal distributions: theory and applications, Statistics, textbooks and monographs, v. 88, Marcel Dekker New York, ISBN 9780824778033. 1988. a, b

Demmel, J. W., Eisenstat, S. C., Gilbert, J. R., Li, X. S., and Liu, J. W. H.: A Supernodal Approach to Sparse Partial Pivoting, SIAM J. Matrix Anal. A., 20, 720–755,, 1999. a

Dolatshah, A., Nelli, F., Bennetts, L. G., Alberello, A., Meylan, M. H., Monty, J. P., and Toffoli, A.: Letter: Hydroelastic interactions between water waves and floating freshwater ice, Phys. Fluids, 30, 091702,, 2018. a

Dumas-Lefebvre, E. and Dumont, D.: Aerial observations of sea ice break-up by ship waves, The Cryosphere Discuss. [preprint],, in review, 2021. a, b, c

Dumont, D., Kohout, A., and Bertino, L.: A wave-based model for the marginal ice zone including a floe breaking parameterization, J. Geophys. Res., 116, C04001,, 2011. a, b, c, d, e

Dupont, F., Dumont, D., Lemieux, J.-F., Dumas-Lefebvre, E., and Caya, A.: A probabilistic seabed–ice keel interaction model, The Cryosphere, 16, 1963–1977,, 2022. a

Fox, C. and Squire, V. A.: Reflection and transmission characteristics at the edge of shore fast sea ice, J. Geophys. Res., 95, 11629,, 1990. a

Fox, C. and Squire, V. A.: On the oblique reflexion and transmission of ocean waves at shore fast sea ice, Philos. T. Roy. Soc. A, 347, 185–218,, 1994. a, b

Herman, A.: Wave-induced stress and breaking of sea ice in a coupled hydrodynamic discrete-element wave–ice model, The Cryosphere, 11, 2711–2725,, 2017. a, b, c, d

Herman, A., Evers, K.-U., and Reimer, N.: Floe-size distributions in laboratory ice broken by waves, The Cryosphere, 12, 685–699,, 2018. a

Herman, A., Wenta, M., and Cheng, S.: Sizes and Shapes of Sea Ice Floes Broken by Waves – A Case Study From the East Antarctic Coast, Front. Earth Sci., 9, 655977,, 2021. a, b, c, d

Horvat, C. and Roach, L. A.: WIFF1.0: a hybrid machine-learning-based parameterization of wave-induced sea ice floe fracture, Geosci. Model Dev., 15, 803–814,, 2022. a, b

Horvat, C. and Tziperman, E.: A prognostic model of the sea-ice floe size and thickness distribution, The Cryosphere, 9, 2119–2134,, 2015. a, b, c, d, e, f, g

Horvat, C. and Tziperman, E.: The evolution of scaling laws in the sea ice floe size distribution, J. Geophys. Res.-Oceans, 122, 7630–7650,, 2017. a

Horvat, C., Tziperman, E., and Campin, J.-M.: Interaction of sea ice floe size, ocean eddies, and sea ice melting, Geophys. Res. Lett., 43, 8083–8090,, 2016. a

Huang, H.-P.: Ice formation in frequently transited navigation channels, PhD thesis, The University of Iowa, ISBN 9798207503257, 1988. a

Hunke, E., Allard, R., Bailey, D. A., Blain, P., Craig, A., Dupont, F., DuVivier, A., Grumbine, R., Hebert, D., Holland, M., Jeffery, N., Jean-Francois Lemieux, Osinski, R., Rasmussen, T., Ribergaard, M., Roberts, A., Francois Roy, Turner, M., and Worthen, D.: CICE-Consortium/CICE: CICE Version 6, Zenodo [code],, 2021. a

Inoue, J.: Ice floe distribution in the Sea of Okhotsk in the period when sea-ice extent is advancing, Geophys. Res. Lett., 31, L20303,, 2004. a

Keller, J. B.: Gravity waves on ice-covered water, J. Geophys. Res.-Oceans, 103, 7663–7669,, 1998. a

Kish, L.: Survey sampling​​​​​​​, John Wiley New York, ISBN 9780471109495, 1965. a

Kohout, A. and Williams, M.: Waves in-ice observations made during the SIPEX II voyage of the Aurora Australis, 2012, Ver. 1, Australian Antarctic Data Centre [data set],, 2015. a

Kohout, A. L. and Meylan, M. H.: An elastic plate model for wave attenuation and ice floe breaking in the marginal ice zone, J. Geophys. Res., 113, C09016,, 2008. a, b, c, d, e, f, g, h

Kolmogoroff, A.: The logarithmically normal law of distribution of dimensions of particles when broken into small parts, in: CR (Doklady) Acad. Sci. URSS (NS), vol. 31, pp. 99–101, 1941. a

Kwok, R.: Arctic sea ice thickness, volume, and multiyear ice coverage: losses and coupled variability (1958–2018), Environ. Res. Lett., 13, 105005,, 2018. a

Kwok, R., Cunningham, G. F., Wensnahan, M., Rigor, I., Zwally, H. J., and Yi, D.: Thinning and volume loss of the Arctic Ocean sea ice cover: 2003–2008, J. Geophys. Res., 114, C07005,, 2009. a

Lilliefors, H. W.: On the Kolmogorov-Smirnov Test for Normality with Mean and Variance Unknown, J. Am. Stat. Assoc., 62, 399–402,, 1967. a

Massey, F. J.: The Kolmogorov-Smirnov Test for Goodness of Fit, J. Am. Stat. Assoc., 46, 68–78,, 1951. a

Massom, R. A., Scambos, T. A., Bennetts, L. G., Reid, P., Squire, V. A., and Stammerjohn, S. E.: Antarctic ice shelf disintegration triggered by sea ice loss and ocean swell, Nature, 558, 383–389,, 2018. a

Meylan, M. and Squire, V. A.: The response of ice floes to ocean waves, J. Geophys. Res., 99, 891–900​​​​​​​,, 1994. a

Meylan, M. H., Bennetts, L. G., and Kohout, A. L.: In situ measurements and analysis of ocean waves in the Antarctic marginal ice zone, Geophys. Res. Lett., 41, 5046–5051,, 2014. a

Mokus, N. and Montiel, F.: Model code and simulation results for the investigation of a wave-generated floe size distribution, Figshare [code and data set],, 2021. a

Mokus, N. and Montiel, F.: Floe size distributions in irregular sea states, in: Proceedings of the 37th International Workshop on Water Waves and Floating Bodies, Giardini Naxos, Italy, 10–13 April 2022, 106–109, ISBN 9788876170539, (last access: 12 September 2022), 2022a. a

Montiel, F. and Mokus, N.: Theoretical framework for the emergent floe size distribution in the marginal ice zone: the case for lognormality, Phil. Trans. R. Soc. A, 380, 20210257,​​​​​​​, 2022b. a, b

Montiel, F. and Squire, V. A.: Modelling wave-induced sea ice break-up in the marginal ice zone, P. Roy. Soc. A-Math. Phy., 473, 20170258,, 2017. a, b, c, d, e, f, g, h

Montiel, F., Bennetts, L., and Squire, V.: The transient response of floating elastic plates to wavemaker forcing in two dimensions, J. Fluid. Struct., 28, 416–433,, 2012. a

Montiel, F., Squire, V. A., and Bennetts, L. G.: Attenuation and directional spreading of ocean wave spectra in the marginal ice zone, J. Fluid Mech., 790, 492–522,, 2016. a, b

Montiel, F., Squire, V. A., Doble, M., Thomson, J., and Wadhams, P.: Attenuation and Directional Spreading of Ocean Waves During a Storm Event in the Autumn Beaufort Sea Marginal Ice Zone, J. Geophys. Res.-Oceans, 123, 5912–5932,, 2018. a

Montiel, F., Kohout, A. L., and Roach, L. A.: Physical Drivers of Ocean Wave Attenuation in the Marginal Ice Zone, J. Phys. Oceanogr., 52, 889–906,, 2022. a

Mosig, J. E. M.: Contemporary wave–ice interaction models, PhD thesis, University of Otago,, 2018. a

Mosig, J. E. M., Montiel, F., and Squire, V. A.: Comparison of viscoelastic-type models for ocean wave attenuation in ice-covered seas, J. Geophys. Res.-Oceans, 120, 6072–6090,, 2015. a, b

Ochi, M. K.: Ocean waves: the Stochastic Approach, vol. 6​​​​​​​, Cambridge University Press, ISBN 9780521017671, 2005. a

Paget, M., Worby, A. P., and Michael, K. J.: Determining the floe-size distribution of East Antarctic sea ice from digital aerial photographs, Ann. Glaciol., 33, 94–100,, 2001. a

Parkinson, C. L. and Comiso, J. C.: On the 2012 record low Arctic sea ice cover: Combined impact of preconditioning and an August storm, Geophys. Res. Lett., 40, 1356–1361,, 2013. a

Passerotti, G., Bennetts, L. G., von Bock und Polach, F., Alberello, A., Puolakka, O., Dolatshah, A., Monbaliu, J., and Toffoli, A.: Interactions between irregular wave fields and sea ice: A physical model for wave attenuation and ice break up, J. Phys. Oceanogr., 52, 1431–1446,, 2022. a

Passerotti, G., Bennetts, L. G., von Bock und Polach, F., Alberello, A., Puolakka, O., Dolatshah, A., Monbaliu, J., and Toffoli, A.: Interactions between Irregular Wave Fields and Sea Ice: A Physical Model for Wave Attenuation and Ice Breakup in an Ice Tank, J. Phys. Oceanogr., 52, 1431–1446,, 2022. a

Perovich, D. K. and Jones, K. F.: The seasonal evolution of sea ice floe size distribution, J. Geophys. Res.-Oceans, 119, 8767–8777,, 2014. a

Pierson, W. J. and Moskowitz, L.: A proposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii, J. Geophys. Res., 69, 5181–5190,, 1964. a

Roach, L. A., Horvat, C., Dean, S. M., and Bitz, C. M.: An Emergent Sea Ice Floe Size Distribution in a Global Coupled Ocean-Sea Ice Model, J. Geophys. Res.-Oceans, 123, 4322–4337,, 2018. a, b, c, d, e, f

Roach, L. A., Bitz, C. M., Horvat, C., and Dean, S. M.: Advances in Modeling Interactions Between Sea Ice and Ocean Surface Waves, J. Adv. Model.Earth Sy., 11, 4167–4181,, 2019. a, b, c

Robinson, N. and Palmer, S.: A modal analysis of a rectangular plate floating on an incompressible liquid, J. Sound Vib. 142, 453–460,, 1990. a

Rothrock, D. A. and Thorndike, A. S.: Measuring the sea ice floe size distribution, J. Geophys. Res., 89, 6477–6486​​​​​​​,, 1984. a, b, c

Santi, F. D. and Olla, P.: Effect of small floating disks on the propagation of gravity waves, Fluid Dyn. Res., 49, 025512,, 2017. a

Squire, V. and Fox, C.: On ice coupled waves: a comparison of data and theory, in: Advances in ice technology: Proc. 3rd Int. Conf. on Ice Technology, Computational Mechanics Publications Cambridge, MA, 269–280, 1992. a

Squire, V. A.: Ocean Wave Interactions with Sea Ice: A Reappraisal, Annu. Review Fluid Mech., 52, 37–60,, 2020. a

Squire, V. A. and Montiel, F.: Evolution of Directional Wave Spectra in the Marginal Ice Zone: A New Model Tested with Legacy Data, J. Phys. Oceanogr., 46, 3121–3137,, 2016. a

Squire, V. A. and Moore, S. C.: Direct measurement of the attenuation of ocean waves by pack ice, Nature, 283, 365–368,, 1980. a

Steele, M.: Sea ice melting and floe geometry in a simple ice-ocean model, J. Geophys. Res.-Oceans, 97, 17729–17738,, 1992. a

Steer, A., Worby, A., and Heil, P.: Observed changes in sea-ice floe size distribution during early summer in the western Weddell Sea, Deep-Sea Res. Pt. II, 55, 933–942,, 2008. a

Stern, H. L., Schweiger, A. J., Zhang, J., and Steele, M.: On reconciling disparate studies of the sea-ice floe size distribution, Elementa: Science of the Anthropocene, 6, 1–16,, 2018. a, b, c

Stroeve, J., Holland, M. M., Meier, W., Scambos, T., and Serreze, M.: Arctic sea ice decline: Faster than forecast, Geophys. Res. Lett., 34, L09501,, 2007. a

Thomson, J. and Rogers, W. E.: Swell and sea in the emerging Arctic Ocean, Geophys. Res. Lett., 41, 3136–3140,, 2014. a

Thorndike, A. S., Rothrock, D. A., Maykut, G. A., and Colony, R.: The thickness distribution of sea ice, J. Geophys. Res., 80, 4501–4513,, 1975. a

Toyota, T., Takatsuji, S., and Nakayama, M.: Characteristics of sea ice floe size distribution in the seasonal ice zone, Geophys. Res. Lett., 33, L02616,, 2006. a

Toyota, T., Haas, C., and Tamura, T.: Size distribution and shape properties of relatively small sea-ice floes in the Antarctic marginal ice zone in late winter, Deep-Sea Res. Pt. II, 58, 1182–1193,, 2011. a

Vaughan, G. L. and Squire, V. A.: Scattering of ice coupled waves by a sea-ice sheet with random thickness, Wave. Random Complex, 17, 357–380,, 2007. a

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, İ., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and SciPy 1.0 Contributors: SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nat. Methods, 17, 261–272,, 2020.  a

Wadhams, P.: Attenuation of swell by sea ice, J. Geophys. Res., 78, 3552–3563,, 1973. a

Wadhams, P., Squire, V. A., Goodman, D. J., Cowan, A. M., and Moore, S. C.: The attenuation rates of ocean waves in the marginal ice zone, J. Geophys. Res., 93, 6799–6818​​​​​​​,, 1988. a

Wang, R. and Shen, H. H.: Gravity waves propagating into an ice-covered ocean: A viscoelastic model, J. Geophys. Res., 115, C06024,, 2010. a

Wang, Y., Holt, B., Rogers, W. E., Thomson, J., and Shen, H. H.: Wind and wave influences on sea ice floe size and leads in the Beaufort and Chukchi Seas during the summer-fall transition 2014, J. Geophys. Res.-Oceans, 121, 1502–1525,, 2016. a

Williams, T. and Porter, R.: The effect of submergence on the scattering by the interface between two semi-infinite sheets, J. Fluid. Struct., 25, 777–793,, 2009. a, b

Williams, T. D., Bennetts, L. G., Squire, V. A., Dumont, D., and Bertino, L.: Wave–ice interactions in the marginal ice zone. Part 1: Theoretical foundations, Ocean Model., 71, 81–91,, 2013. a, b, c, d, e, f, g, h

Williams, T. D., Rampal, P., and Bouillon, S.: Wave–ice interactions in the neXtSIM sea-ice model, The Cryosphere, 11, 2117–2135,, 2017. a, b

Williams, T. D. C.: Reflections on ice: scattering of flexural gravity waves by irregularities in Arctic and Antarctic ice sheets, PhD thesis, University of Otago, (last access: 12 September 2022), 2006. a

Zhang, J., Schweiger, A., Steele, M., and Stern, H.: Sea ice floe size distribution in the marginal ice zone: Theory and numerical experiments, J. Geophys. Res.-Oceans, 120, 3484–3498,, 2015. a

Zhang, L., Delworth, T. L., Cooke, W., and Yang, X.: Natural variability of Southern Ocean convection as a driver of observed climate trends, Nat. Clim. Change, 9, 59–65,, 2018. a

Short summary
On the fringes of polar oceans, sea ice is easily broken by waves. As small pieces of ice, or floes, are more easily melted by the warming waters than a continuous ice cover, it is important to incorporate these floe sizes in climate models. These models simulate climate evolution at the century scale and are built by combining specialised modules. We study the statistical distribution of floe sizes under the impact of waves to better understand how to connect sea ice modules to wave modules.