Wavetriggered breakup in the marginal ice zone generates lognormal floe size distributions: a simulation study
 Department of Mathematics and Statistics, University of Otago, Dunedin, New Zealand
 Department of Mathematics and Statistics, University of Otago, Dunedin, New Zealand
Correspondence: Nicolas Guillaume Alexandre Mokus (nmokus@maths.otago.ac.nz)
Hide author detailsCorrespondence: Nicolas Guillaume Alexandre Mokus (nmokus@maths.otago.ac.nz)
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 twodimensional numerical model of waveinduced sea ice breakup to estimate the FSD resulting from repeated fracture events. This model, based on linear water wave theory and viscoelastic sea ice rheology, solves for the scattering of an incoming timeharmonic 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.
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 (Squire, 2020). 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 recordlow Arctic sea ice extent was amplified by wave activity (Parkinson and Comiso, 2013). With thinner and weaker firstyear ice becoming dominant in the Arctic (Kwok et al., 2009; Kwok, 2018), the sea ice grows more vulnerable to flexureinduced failure.
The marginal ice zone (hereafter MIZ), a belt of loosely to densely packed ice floes, serves as a buffer between the icefree 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 (Steele, 1992; Perovich and Jones, 2014), 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 Rogers, 2014).
Remote sensing observations of floe sizes (e.g. Rothrock and Thorndike, 1984; 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 Thorndike, 1984). 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 (Wadhams, 1973). A component of this attenuation is scattering, which results from inhomogeneities of the ice cover such as cracks and leads (Williams and Porter, 2009), changes in thickness as caused by pressure ridges (Vaughan and Squire, 2007), and floes themselves (Meylan and Squire, 1994; Kohout and Meylan, 2008; Bennetts and Squire, 2011; 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 Meylan, 2008). Dissipative terms can be included in these formulations (Squire and Fox, 1992; Fox and Squire, 1994). Alternatively, rheological models seek to represent the ice cover as a viscous continuum whose properties are representative of a discontinuous field (Keller, 1998; Wang and Shen, 2010; Santi and Olla, 2017). 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 Moore, 1980; 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 waveinduced 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 Squire, 2017; Herman, 2017) 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 Squire, 2017) or strain (Kohout and Meylan, 2008; Williams et al., 2013; Horvat and Tziperman, 2015; 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 standalone 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 waveinduced sea ice breakup in largescale 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 stressbased failure criterion. They investigated the FSD obtained after repeated breakup events and found that nearnormal 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 nonhydrostatic, nonlinear 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 smallerscale 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 Squire, 2017) or indirectly through the use of discrete elements (Herman, 2017), effectively ensuing discrete distributions. Recent field observations of floes directly impacted by waves (DumasLefebvre and Dumont, 2021; 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 waveinduced 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 twodimensional setting (one horizontal, one vertical), relying on an established scattering theory, augmented with an energydissipating process. We then incorporate a strainbased breakup parametrisation. We let an FSD emerge by repeatedly breaking off floes from a semiinfinite 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 Tziperman, 2015).
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.
We consider surface gravity waves propagating in a twodimensional 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:
We place an array of N_{f}+1 nonoverlapping ice floes, modelled as floating viscoelastic 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=\frac{Y{h}^{\mathrm{3}}}{\mathrm{12}\left(\mathrm{1}{\mathit{\nu}}^{\mathrm{2}}\right)}$ (where Y and ν are respectively Young's modulus and Poisson's ratio), and viscosity γ; their draught is $d=\frac{\mathit{\rho}}{{\mathit{\rho}}_{w}}h$. Floe j, where $j\in \left\{\mathrm{0},\mathrm{\dots},{N}_{f}\right\}$, is located in space by the horizontal coordinate of its left edge x_{j} (ordered so that ${x}_{j}<{x}_{j+\mathrm{1}}$) and its length L_{j}, with ${L}_{{N}_{f}}$ being infinite. At rest, the fluid region covered by floe j is encompassed in the subdomain ${\mathrm{\Omega}}_{j}^{f}=[{x}_{j},{x}_{j}+{L}_{j}]\times [H,d]$. The interface with the ice $\partial {\mathrm{\Omega}}_{j}^{f}$ is on $z=d$. We denote ${\mathrm{\Omega}}^{f}={\bigcup}_{\mathrm{0}}^{{N}_{f}}{\mathrm{\Omega}}_{j}^{f}$ and $\partial {\mathrm{\Omega}}^{f}={\bigcup}_{\mathrm{0}}^{{N}_{f}}\partial {\mathrm{\Omega}}_{j}^{f}$. The icefree subdomain left of the floecovered subdomain ${\mathrm{\Omega}}_{j}^{f}$ is ${\mathrm{\Omega}}_{j}^{w}$ so that at rest,
The horizontal bottom boundary ∂Ω_{H} is at $z=H$, and the interface $\partial {\mathrm{\Omega}}_{j}^{w}$ between the atmosphere and ${\mathrm{\Omega}}_{j}^{w}$ 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 $\mathrm{\Omega}={\mathrm{\Omega}}^{w}\cup {\mathrm{\Omega}}^{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=\mathit{\eta}(x,t)$, whether it is an interface with the atmosphere (in which case η≈0) or with an ice floe (in which case $\mathit{\eta}\approx d$). We further set $\mathrm{\Phi}=\mathrm{Re}\left[\mathit{\varphi}(x,z)\mathrm{exp}(i\mathit{\omega}t)\right]$ with ϕ a timeindependent, complexvalued function.
We consider the seabed to be impervious, hence not allowing for normal flow, so that on ∂Ω_{H}
The small amplitude forcing allows us to use linear surface wave theory in Ω^{w}, leading to the boundary condition
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
on ∂Ω^{f}. Equations (4) and (5) stem from the assumptions that under smallamplitude 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
on $\left\{{x}_{j},{x}_{j}+{L}_{j}\right\}\times [d,\mathrm{0}]$.
We finally add the freeedge conditions
on $\left\{{x}_{j},{x}_{j}+{L}_{j}\right\}\times \left\{d\right\}$, 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 timeindependent form, ∇^{2}ϕ=0, complete our boundary value problem.
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 zerodraught 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 nonzero 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 subdomain ${\mathrm{\Omega}}_{j}^{w}$ or ${\mathrm{\Omega}}_{j}^{f}$, the velocity potential is decomposed as the superposition of forwardtravelling and backwardtravelling plane waves, using the ansatz $\left[{c}^{+}{e}^{i\left(kx{\mathit{\theta}}^{+}\right)}+{c}^{}{e}^{i\left(kx{\mathit{\theta}}^{}\right)}\right]\mathit{\zeta}\left(z\right)$, where ${c}^{+},{c}^{}\in \mathbb{C}$ are coefficients to be determined, and ${\mathit{\theta}}^{+},{\mathit{\theta}}^{}\in [\mathrm{0},\mathrm{2}\mathit{\pi})$ are phase shifts introduced to simplify analytical derivations and improve numerical stability. Solving the boundary value problem described in Sect. 2 in all subdomains, the potential is expanded into series of travelling and evanescent wave modes
where superscripts w and f are related to an openwater and floecovered subdomain.
The wavenumbers {${k}_{n}^{w}\mid n\ge \mathrm{0}\mathit{\}}$ are the roots of the dispersion relation
such that ${k}_{\mathrm{0}}^{w}$ is a positive real number (therefore associated with left and rightpropagating wave modes in Ω^{w}), while {${k}_{n}^{w}\mid n>\mathrm{0}\mathit{\}}$ are purely imaginary numbers with a positive imaginary part, sorted by ascending imaginary part (associated with exponentially decaying evanescent wave modes).
Likewise, the wavenumbers $\mathit{\{}{k}_{n}^{f}\mid n\ge \mathrm{2}\mathit{\}}$ are the roots of the dispersion relation
such that $\mathit{\{}{k}_{n}^{f}\mid n>\mathrm{2}\mathit{\}}$ are complex numbers in the first quadrant of the complex plane, sorted by ascending imaginary part for n>0, and ${k}_{\mathrm{2}}^{f}$ is a complex number in the second quadrant of the complex plane.
Since $\mathcal{O}\left[\mathrm{Re}\left({k}_{\mathrm{0}}^{f}\right)\right]\gg \mathcal{O}\left[\mathrm{Im}\left({k}_{\mathrm{0}}^{f}\right)\right]$, ${k}_{\mathrm{0}}^{f}$ is associated with left and rightpropagating wave modes in Ω^{f}. In contrast, $\mathcal{O}\left[\mathrm{Re}\left({k}_{n}^{f}\right)\right]\ll \mathcal{O}\left[\mathrm{Im}\left({k}_{n}^{f}\right)\right]$ for n>0: these modes are associated with exponentially decaying wave modes. Finally, $\mathcal{O}\left[\mathrm{Re}\left({k}_{n}^{f}\right)\right]=\mathcal{O}\left[\mathrm{Im}\left({k}_{n}^{f}\right)\right]$ for $n\in \mathit{\{}\mathrm{2},\mathrm{1}\mathit{\}}$, 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), ${k}_{\mathrm{0}}^{f}$ is a positive real number, ${k}_{\mathrm{1}}^{f}$ is in the first quadrant of the complex plane, ${k}_{\mathrm{2}}^{f}=\stackrel{\mathrm{\u203e}}{{k}_{\mathrm{1}}^{f}}$, and $\mathit{\{}{k}_{n}^{f}\mid n>\mathrm{0}\mathit{\}}$ are purely imaginary numbers with a positive imaginary part. We note that when $\mathrm{1}\frac{{\mathit{\omega}}^{\mathrm{2}}d}{g}$ becomes negative (large frequencies or thickness), instead of having the three distinct roots ${k}_{\mathrm{2}}^{f},{k}_{\mathrm{1}}^{f}$, and ${k}_{\mathrm{0}}^{f}$, Eq. (10) may admit one double root or one triple root (Williams, 2006, p. 39) as the complex roots may become purely imaginary. Therefore, we enforce $\mathit{\omega}\le \sqrt{\frac{g}{d}}$ 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 ${c}_{j,n}^{w}$, and ${c}_{j,n}^{f+}$) and two from its right edge (with coefficients ${c}_{j,n}^{w+}$, and ${c}_{j,n}^{f}$), except for the leftmost, semiinfinite floe (j=N_{f}), whose right edge is ignored. The coefficients ${c}_{\mathrm{1},n}^{w+}$ and ${c}_{{N}_{f},n}^{f}$ are prescribed and represent the righttravelling incident wave forcing in ${\mathrm{\Omega}}_{\mathrm{0}}^{w}$ and the absence of forcing in ${\mathrm{\Omega}}_{{N}_{f}}^{f}$, respectively: only ${c}_{\mathrm{1},\mathrm{0}}^{w+}=i\frac{g}{\mathit{\omega}}a$ is nonzero. The quantity $\mathit{\theta}={\mathit{\theta}}_{\mathrm{1},\mathrm{0}}^{w+}$ is an arbitrary phase associated with the forcing. We also note that $\mathit{\{}{\mathit{\theta}}_{\mathrm{1},n}^{w+}\mid n>\mathrm{0}\mathit{\}}$ and $\mathit{\{}{\mathit{\theta}}_{{N}_{f},n}^{f}\mid n\ge \mathrm{2}\mathit{\}}$ 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.
Finally, the functions
are vertical basis functions in the freesurface subdomains and the icecovered subdomains, 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 icefree and icecovered subdomains, i.e. at each floe edge. These conditions are enforced by matching ϕ and $u=\frac{\partial \mathit{\varphi}}{\partial 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 ${c}_{m\mathrm{1},n}^{w+}$ and ${c}_{m,n}^{f}$, 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, ${c}_{j,n}^{f+}$ and ${c}_{j,n}^{w}$. When truncating the series in Eq. (8) to N_{v} evanescent modes, the relations can be summarised by the matrix equation
where ${\mathbf{T}}^{fw}\in {\mathbb{C}}_{{N}_{v}+\mathrm{3}\times {N}_{v}+\mathrm{1}}$, ${\mathbf{T}}^{wf}\in {\mathbb{C}}_{{N}_{v}+\mathrm{1}\times {N}_{v}+\mathrm{3}}$, ${\mathbf{R}}^{f}\in {\mathbb{C}}_{{N}_{v}+\mathrm{3}\times {N}_{v}+\mathrm{3}}$, and ${\mathbf{R}}^{w}\in {\mathbb{C}}_{{N}_{v}+\mathrm{1}\times {N}_{v}+\mathrm{1}}$ are matrices, respectively describing transmission and reflection of waves in either direction through the floe edge; ${\mathit{c}}_{j}^{w\pm}={\begin{array}{c}{c}_{j,\mathrm{0}}^{w\pm},\mathrm{\dots},{c}_{j,{N}_{v}}^{w\pm}\end{array}}^{\mathrm{T}}$, ${\mathit{c}}_{j}^{f\pm}={\begin{array}{c}{c}_{j,\mathrm{2}}^{f\pm},\mathrm{\dots},{c}_{j,{N}_{v}}^{f\pm}\end{array}}^{\mathrm{T}}$ are vectors of unknown coefficients; and ${\mathbf{S}}^{w}\in {\mathbb{C}}_{{N}_{v}+\mathrm{1}\times {N}_{v}+\mathrm{1}}$ and ${\mathbf{S}}^{f}\in {\mathbb{C}}_{{N}_{v}+\mathrm{3}\times {N}_{v}+\mathrm{3}}$ are diagonal phase shift matrices.
By symmetry, the scattering by the right edge of floe j is described by
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 N_{v}=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
where M is a tridiagonal block matrix and c the vector of unknown potential coefficients.
An array of N_{f}+1 finite floes leads to the matrix
where each block element of $\stackrel{\mathrm{\u0303}}{\mathbf{M}}$ is a square matrix of size 4(N_{v}+2), and 0 denotes 0filled 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 $\stackrel{\mathrm{\u0303}}{\mathbf{M}}$ from its last 2(N_{v}+2) rows and columns. While the size of M is ${\left[\mathrm{2}\left({N}_{v}+\mathrm{2}\right)\left(\mathrm{2}{N}_{f}+\mathrm{1}\right)\right]}^{\mathrm{2}}$, it has only $\mathrm{4}\left({N}_{v}+\mathrm{2}\right)\left[{N}_{f}\left(\mathrm{2}{N}_{v}+\mathrm{5}\right)+\frac{\mathrm{1}}{\mathrm{2}}\right]$ nonzero elements: this number grows linearly with N_{f}, instead of quadratically.
The vector of unknown coefficients is
and the forcing term
Building M and f and solving for c is linear in time for 𝒪(N_{f})>10. For N_{f}=10^{5}, these operations are done in around 10 ms on an Intel Core i56300U 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 strainbased approach (e.g. Kohout and Meylan, 2008; Horvat and Tziperman, 2015). 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 ${\stackrel{\mathrm{\u0303}}{\mathit{\epsilon}}}_{j}$ undergone by ice floe j is
where $\left({x}^{\prime},{z}^{\prime}\right)$ describes a coordinate system local to the floe, defined as ${x}^{\prime}=x{x}_{j}$ and ${z}^{\prime}=z\left(\frac{h}{\mathrm{2}}d\right)=zh\left(\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\rho}}{{\mathit{\rho}}_{\mathrm{0}}}\right)$, 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)=\mathit{\eta}(x,t)$. As $\left{z}^{\prime}\right\le \frac{h}{\mathrm{2}}$, the maximum (in absolute value) strain, ε_{j}, is located on either surface of the floe, i.e. ${z}^{\prime}=\pm \frac{h}{\mathrm{2}}$. It follows that
with $T=\frac{\mathrm{2}\mathit{\pi}}{\mathit{\omega}}$ the wave period.
Floe j is set to break when
where ε_{c} is an empirically determined strain threshold.
We situate the breakup point at
so that a floe of length L_{j} is turned into two floes of length x_{b} and L_{j}−x_{b}. Hence, the number of floes at most doubles, if all the floes break in a single breakup simulation.
3.3 Numerical experiment setup
The values of the parameters kept fixed across all simulations are given in Table 1. Our experiments unwind as follows:

The model is initialised with a set of physical parameters as input and a single, semiinfinite floe.

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.

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 coordinatelength pairs describing the final geometry.
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 orderpreserving 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 semiinfinite 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 floedependent.
The final result extracted from the simulation is the set of newly formed floe lengths, excluding the semiinfinite floe on the right of the domain, considered to be a steadystate FSD.
We first investigate the FSD our model generates under monochromatic forcing with prescribed wave period T, corresponding to angular frequency $\mathit{\omega}=\frac{\mathrm{2}\mathit{\pi}}{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 Meylan, 2008).
We obtain unimodal, rightskewed 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 rightskewed 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 socalled 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 socalled 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 freeedge 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 halfwavelength (respectively 12, 50, and 112 m for periods of 4, 8, and 12 s), as observed by Herman et al. (2021).
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 nonlinear 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 Squire, 1990); 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 lowperiod 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 semiinfinite cover, which leads to these higher minima.
The breakup width, here defined as the cumulated length of ice broken off the semiinfinite 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.
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 f_{L}(l), we take the weighted average of distributions ${\stackrel{\mathrm{\u0303}}{f}}_{L}(l,\mathit{\omega})$ resulting from monochromatic model runs, with weights taken from the energy density S(ω) of a theoretical ocean spectrum over a truncated frequency interval $[{\mathit{\omega}}_{min},{\mathit{\omega}}_{max}]$ so that
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 Moskowitz, 1964) as it can be easily adapted to depend only on the significant wave height H_{s} (Ochi, 2005), giving
where ${c}_{\mathrm{1}}=\mathrm{8.1}\times {\mathrm{10}}^{\mathrm{3}}$ and ${c}_{\mathrm{2}}=\mathrm{3.24}\times {\mathrm{10}}^{\mathrm{2}}$ are nondimensional constants. This spectrum has been used in previous wave–seaice interaction studies (e.g. Kohout and Meylan, 2008; Dumont et al., 2011) and is a reduced version of the twoparameter Bretschneider spectrum, used in this context as well (e.g. Horvat and Tziperman, 2015, 2017; Montiel and Squire, 2017).
We evaluate Eq. (23) numerically on 200 linearly spaced frequency bins, setting H_{s}=2a. As by definition
the bounds of integration ω_{min} and ω_{min} are set so that the tails $\frac{\mathrm{16}}{{{H}_{\mathrm{s}}}^{\mathrm{2}}}{\int}_{\mathrm{0}}^{{\mathit{\omega}}_{min}}S\left(\mathit{\omega}\right)\mathrm{d}\mathit{\omega}=\frac{\mathrm{16}}{{{H}_{\mathrm{s}}}^{\mathrm{2}}}{\int}_{{\mathit{\omega}}_{max}}^{+\mathrm{\infty}}S\left(\mathit{\omega}\right)\mathrm{d}\mathit{\omega}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{7}}$, which captures a significant part of the spectrum. If necessary, ω_{max} is adjusted to ensure $\mathit{\omega}\le \sqrt{\frac{g}{d}}$, as discussed in Sect. 3.1.
For each of the 200 frequencies, we draw an FSD ${\stackrel{\mathrm{\u0303}}{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 f_{L}. As we have 50 independent realisations for each of the 200 frequencies, we have virtually infinitely many ways to build f_{L} (${\mathrm{50}}^{\mathrm{200}}\approx \mathrm{6}\times {\mathrm{10}}^{\mathrm{339}}$).
5.2 Reference configuration
We consider a reference configuration where H_{s}=1 m, h=1 m, and ${\mathit{\epsilon}}_{c}=\mathrm{4}\times {\mathrm{10}}^{\mathrm{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.
The resulting FSD, shown in Fig. 5a, is remarkably well fitted by a threeparameter 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}:
The associated density function f_{L}(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
More details on the lognormal distribution can be found in Crow and Kunio (1988). In the following, we use the notation $\mathit{\theta}=\left(s,\mathit{\sigma},\mathit{\tau}\right)$ for the parameter vector.
We obtain a point estimate $\widehat{\mathit{\theta}}$ with maximum likelihood estimation (see e.g. Azzalini, 1996; Crow and Kunio, 1988). Among the ensemble, $\widehat{\mathit{\theta}}$ seems to be normally distributed, with strong correlations between the three parameters and small variances. Therefore, we use the mean vector $\stackrel{\mathrm{\u203e}}{\mathit{\theta}}$ 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 $\stackrel{\mathrm{\u203e}}{\mathit{\theta}}$, we ensure this model is preserved. This can be justified by the low spread of $\widehat{\mathit{\theta}}$. 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: $\mathrm{4}\times {\mathrm{10}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathit{\%}$; median for σ: $\mathrm{4}\times {\mathrm{10}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathit{\%}$).
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 Subdomain 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 ${\stackrel{\mathrm{\u0303}}{f}}_{L;{b}_{inf}\text{\u2013}{b}_{sup}}$, with $\mathrm{0}\le {b}_{inf}<{b}_{sup}\le \mathrm{1}$, not on all floe lengths $\left\{{L}_{j}\mid j\in \left\{\mathrm{0},\mathrm{\dots},{N}_{f}\mathrm{1}\right\}\right\}$ but on the subset
where ${\mathrm{\Lambda}}_{j}={\sum}_{m=\mathrm{0}}^{j}{L}_{m}$ is the cumulated floe size up to floe j, and $\mathrm{\Lambda}={\mathrm{\Lambda}}_{{N}_{f}\mathrm{1}}$ is the total length of finite ice in the domain. We ignore the openwater gaps in the definition of the breakup width. These monochromatic densities are then combined as detailed in Sect. 5.1. The difference ${b}_{sup}{b}_{inf}$ 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.
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 semiinfinite 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 b_{inf} 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 H_{s}, 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 $\widehat{\mathit{\tau}}<\mathrm{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.
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 lowamplitude 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 H_{s}=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 smallperiod 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 highperiod 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 H_{s} in the 1.4 to 2.2 m range, but the same behaviour is visible on the histograms (Fig. C1). This nonmonotonic 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 (DumasLefebvre and Dumont, 2021; 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 $\stackrel{\mathrm{\u203e}}{\mathit{\theta}}$ 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 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 ${\mathit{\epsilon}}_{c}\ge \mathrm{3}\times {\mathrm{10}}^{\mathrm{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 H_{s}=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.
The emergence of a lognormal FSD from repeated waveinduced 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 processbased sea ice breakup modelling (Herman, 2017; Montiel and Squire, 2017). 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 (Huang, 1988; Bonath et al., 2019). One of its earliest applications was the description of particle sizes from repeated fragmentation events (Kolmogoroff, 1941). 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 Mokus, 2022b).
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 allaround 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 postwaveinduced 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 waveinduced 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 waveinduced breakup. Another point of concern is that we have $\widehat{\mathit{\tau}}<\mathrm{0}$ for H_{s}≥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 ($<{\mathrm{10}}^{\mathrm{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 nonlinear 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 strainaveraged 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 Squire, 2017). Ultimately, the breakup of an ice cover in response to a wave event is a transient process. Consequently, validations of spectrumgenerated FSDs will need to be sought against timedependent simulations or experiments that control the wave forcing. Recent work by DumasLefebvre 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 rightskewed 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 $\mathit{\gamma}=\mathrm{13}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa}\phantom{\rule{0.125em}{0ex}}\mathrm{s}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}$ derived from a 1979 campaign in the Bering Sea, while Mosig et al. (2015) fitted $\mathit{\gamma}=\mathrm{6.9}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa}\phantom{\rule{0.125em}{0ex}}\mathrm{s}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}$ to a 2012 Southern Ocean dataset (Kohout and Williams, 2015). Massom et al. (2018) derived $\mathit{\gamma}=\mathrm{13.5}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa}\phantom{\rule{0.125em}{0ex}}\mathrm{s}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}$ from the same dataset. We used $\mathit{\gamma}=\mathrm{20}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa}\phantom{\rule{0.125em}{0ex}}\mathrm{s}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{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 Montiel, 2016).
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 subgrid 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 waveinduced 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 waveinduced 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 goodnessoffit 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 lognormal and powerlaw 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 waveinduced 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 FSDevolving models.
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 $\mathit{\gamma}=\mathrm{20}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa}\phantom{\rule{0.125em}{0ex}}\mathrm{s}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}$. Summary statistics and CDFs are presented in Table A1 and Fig. A2.
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 highfrequency waves (5 to 15 m), which is unrealistic in the context of our model.
The Kolmogorov–Smirnov statistic D_{KS} is the largest difference between an empirical cumulative distribution function (CDF) and a reference CDF (Massey, 1951). By definition of the CDF, D_{KS} is bounded by 0 and 1. When the distribution parameters have been estimated from the data, D_{KS} can be used to run a Lilliefors test (Lilliefors, 1967). By comparing D_{KS} 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 D_{KS} as an indicative performance metric to compare our different configurations. More specifically, for each fitted lognormal with estimated parameters $\widehat{\mathit{\theta}}$ (that is, 500 different realisations per model configuration), we generate a random sample of size N_{eff} from the distribution. We use Kish's effective sample size (Kish, 1965), the roundedup ratio of the squared sum of weights to the sum of the squared weights, as N_{eff}. We use maximum likelihood estimation to estimate ${\widehat{\mathit{\theta}}}_{b}$ from this sample, and we compute D_{KS} for this sample and the distribution parametrised by ${\widehat{\mathit{\theta}}}_{b}$. We then define ΔD_{KS} as the difference between D_{KS} from the random sample and D_{KS} 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 ΔD_{KS} is bounded by −1 and 1, with ΔD_{KS}>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. D1–D3. This procedure quantitatively illustrates the conclusions derived from analysing histograms and quantile–quantile plots.
The model code, software tools developed for analysis, and the resulting preprocessed output presented in this paper are publicly available at https://doi.org/10.6084/m9.figshare.17303927 (Mokus and Montiel, 2021). The raw output is available upon request.
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.
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).
This research has been supported by the University of Otago (doctoral scholarship), the Marsden Fund (grant nos. 18UOO216 and 20UOO173), and Antarctica New Zealand (Antarctic Science Platform Project 4).
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, https://doi.org/10.1029/2011jc007221, 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, https://doi.org/10.5194/tc144032020, 2020. a, b, c
Bennetts, L. G. and Squire, V. A.: On the calculation of an attenuation coefficient for transects of icecovered ocean, P. Roy. Soc. AMath. Phy., 468, 136–162, https://doi.org/10.1098/rspa.2011.0155, 2011. a
Bennetts, L. G., O'Farrell, S., and Uotila, P.: Brief communication: Impacts of oceanwaveinduced breakup of Antarctic sea ice via thermodynamics in a standalone version of the CICE seaice model, The Cryosphere, 11, 1035–1040, https://doi.org/10.5194/tc1110352017, 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 03766756, https://www.poac.com/Papers/2019/pdf/POAC19106.pdf (last access: 12 September 2022), 2019. a
Boutin, G., Ardhuin, F., Dumont, D., Sévigny, C., GirardArdhuin, F., and Accensi, M.: Floe Size Effect on WaveIce Interactions: Possible Effects, Implementation in Wave Model, and Evaluation, J. Geophys. Res.Oceans, 123, 4779–4805, https://doi.org/10.1029/2017jc013622, 2018. a, b
Boutin, G., Lique, C., Ardhuin, F., Rousset, C., Talandier, C., Accensi, M., and GirardArdhuin, F.: Towards a coupled model to investigate wave–sea ice interactions in the Arctic marginal ice zone, The Cryosphere, 14, 709–735, https://doi.org/10.5194/tc147092020, 2020. a, b, c
Boutin, G., Williams, T., Rampal, P., Olason, E., and Lique, C.: Wave–seaice interactions in a brittle rheological framework, The Cryosphere, 15, 431–457, https://doi.org/10.5194/tc154312021, 2021. a, b, c
Castruccio, F. S., RuprichRobert, 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, https://doi.org/10.1175/jclid180307.1, 2019. a
Clauset, A., Shalizi, C. R., and Newman, M. E. J.: PowerLaw Distributions in Empirical Data, SIAM Review, 51, 661–703, https://doi.org/10.1137/070710111, 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, https://doi.org/10.1002/2015gl063063, 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, https://doi.org/10.1137/s0895479895291765, 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, https://doi.org/10.1063/1.5050262, 2018. a
DumasLefebvre, E. and Dumont, D.: Aerial observations of sea ice breakup by ship waves, The Cryosphere Discuss. [preprint], https://doi.org/10.5194/tc2021328, in review, 2021. a, b, c
Dumont, D., Kohout, A., and Bertino, L.: A wavebased model for the marginal ice zone including a floe breaking parameterization, J. Geophys. Res., 116, C04001, https://doi.org/10.1029/2010jc006682, 2011. a, b, c, d, e
Dupont, F., Dumont, D., Lemieux, J.F., DumasLefebvre, E., and Caya, A.: A probabilistic seabed–ice keel interaction model, The Cryosphere, 16, 1963–1977, https://doi.org/10.5194/tc1619632022, 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, https://doi.org/10.1029/jc095ic07p11629, 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, https://doi.org/10.1098/rsta.1994.0044, 1994. a, b
Herman, A.: Waveinduced stress and breaking of sea ice in a coupled hydrodynamic discreteelement wave–ice model, The Cryosphere, 11, 2711–2725, https://doi.org/10.5194/tc1127112017, 2017. a, b, c, d
Herman, A., Evers, K.U., and Reimer, N.: Floesize distributions in laboratory ice broken by waves, The Cryosphere, 12, 685–699, https://doi.org/10.5194/tc126852018, 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, https://doi.org/10.3389/feart.2021.655977, 2021. a, b, c, d
Horvat, C. and Roach, L. A.: WIFF1.0: a hybrid machinelearningbased parameterization of waveinduced sea ice floe fracture, Geosci. Model Dev., 15, 803–814, https://doi.org/10.5194/gmd158032022, 2022. a, b
Horvat, C. and Tziperman, E.: A prognostic model of the seaice floe size and thickness distribution, The Cryosphere, 9, 2119–2134, https://doi.org/10.5194/tc921192015, 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, https://doi.org/10.1002/2016jc012573, 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, https://doi.org/10.1002/2016gl069742, 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., JeanFrancois Lemieux, Osinski, R., Rasmussen, T., Ribergaard, M., Roberts, A., Francois Roy, Turner, M., and Worthen, D.: CICEConsortium/CICE: CICE Version 6, Zenodo [code], https://doi.org/10.5281/zenodo.1205674, 2021. a
Inoue, J.: Ice floe distribution in the Sea of Okhotsk in the period when seaice extent is advancing, Geophys. Res. Lett., 31, L20303, https://doi.org/10.1029/2004gl020809, 2004. a
Keller, J. B.: Gravity waves on icecovered water, J. Geophys. Res.Oceans, 103, 7663–7669, https://doi.org/10.1029/97jc02966, 1998. a
Kish, L.: Survey sampling, John Wiley New York, ISBN 9780471109495, 1965. a
Kohout, A. and Williams, M.: Waves inice observations made during the SIPEX II voyage of the Aurora Australis, 2012, Ver. 1, Australian Antarctic Data Centre [data set], https://doi.org/10.4225/15/53266BEC9607F, 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, https://doi.org/10.1029/2007jc004434, 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, https://doi.org/10.1088/17489326/aae3ec, 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, https://doi.org/10.1029/2009jc005312, 2009. a
Lilliefors, H. W.: On the KolmogorovSmirnov Test for Normality with Mean and Variance Unknown, J. Am. Stat. Assoc., 62, 399–402, https://doi.org/10.1080/01621459.1967.10482916, 1967. a
Massey, F. J.: The KolmogorovSmirnov Test for Goodness of Fit, J. Am. Stat. Assoc., 46, 68–78, https://doi.org/10.1080/01621459.1951.10500769, 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, https://doi.org/10.1038/s4158601802121, 2018. a
Meylan, M. and Squire, V. A.: The response of ice floes to ocean waves, J. Geophys. Res., 99, 891–900, https://doi.org/10.1029/93jc02695, 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, https://doi.org/10.1002/2014gl060809, 2014. a
Mokus, N. and Montiel, F.: Model code and simulation results for the investigation of a wavegenerated floe size distribution, Figshare [code and data set], https://doi.org/10.6084/m9.figshare.17303927, 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, http://www.iwwwfb.org/Abstracts/iwwwfb37/IWWWFB37GLOBAL027.pdf (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, https://doi.org/10.1098/rsta.2021.0257, 2022b. a, b
Montiel, F. and Squire, V. A.: Modelling waveinduced sea ice breakup in the marginal ice zone, P. Roy. Soc. AMath. Phy., 473, 20170258, https://doi.org/10.1098/rspa.2017.0258, 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, https://doi.org/10.1016/j.jfluidstructs.2011.10.007, 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, https://doi.org/10.1017/jfm.2016.21, 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, https://doi.org/10.1029/2018jc013763, 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, https://doi.org/10.1175/jpod210240.1, 2022. a
Mosig, J. E. M.: Contemporary wave–ice interaction models, PhD thesis, University of Otago, http://hdl.handle.net/10523/7958, 2018. a
Mosig, J. E. M., Montiel, F., and Squire, V. A.: Comparison of viscoelastictype models for ocean wave attenuation in icecovered seas, J. Geophys. Res.Oceans, 120, 6072–6090, https://doi.org/10.1002/2015jc010881, 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 floesize distribution of East Antarctic sea ice from digital aerial photographs, Ann. Glaciol., 33, 94–100, https://doi.org/10.3189/172756401781818473, 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, https://doi.org/10.1002/grl.50349, 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, https://doi.org/10.1175/jpod210238.1, 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, https://doi.org/10.1175/jpod210238.1, 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, https://doi.org/10.1002/2014jc010136, 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, https://doi.org/10.1029/jz069i024p05181, 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 OceanSea Ice Model, J. Geophys. Res.Oceans, 123, 4322–4337, https://doi.org/10.1029/2017jc013692, 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, https://doi.org/10.1029/2019ms001836, 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, https://doi.org/10.1016/0022460x(90)90661i, 1990. a
Rothrock, D. A. and Thorndike, A. S.: Measuring the sea ice floe size distribution, J. Geophys. Res., 89, 6477–6486, https://doi.org/10.1029/jc089ic04p06477, 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, https://doi.org/10.1088/18737005/aa59e1, 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, https://doi.org/10.1146/annurevfluid010719060301, 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, https://doi.org/10.1175/jpod160118.1, 2016. a
Squire, V. A. and Moore, S. C.: Direct measurement of the attenuation of ocean waves by pack ice, Nature, 283, 365–368, https://doi.org/10.1038/283365a0, 1980. a
Steele, M.: Sea ice melting and floe geometry in a simple iceocean model, J. Geophys. Res.Oceans, 97, 17729–17738, https://doi.org/10.1029/92jc01755, 1992. a
Steer, A., Worby, A., and Heil, P.: Observed changes in seaice floe size distribution during early summer in the western Weddell Sea, DeepSea Res. Pt. II, 55, 933–942, https://doi.org/10.1016/j.dsr2.2007.12.016, 2008. a
Stern, H. L., Schweiger, A. J., Zhang, J., and Steele, M.: On reconciling disparate studies of the seaice floe size distribution, Elementa: Science of the Anthropocene, 6, 1–16, https://doi.org/10.1525/elementa.304, 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, https://doi.org/10.1029/2007gl029703, 2007. a
Thomson, J. and Rogers, W. E.: Swell and sea in the emerging Arctic Ocean, Geophys. Res. Lett., 41, 3136–3140, https://doi.org/10.1002/2014gl059983, 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, https://doi.org/10.1029/jc080i033p04501, 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, https://doi.org/10.1029/2005gl024556, 2006. a
Toyota, T., Haas, C., and Tamura, T.: Size distribution and shape properties of relatively small seaice floes in the Antarctic marginal ice zone in late winter, DeepSea Res. Pt. II, 58, 1182–1193, https://doi.org/10.1016/j.dsr2.2010.10.034, 2011. a
Vaughan, G. L. and Squire, V. A.: Scattering of ice coupled waves by a seaice sheet with random thickness, Wave. Random Complex, 17, 357–380, https://doi.org/10.1080/17455030701250467, 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, https://doi.org/10.1038/s4159201906862, 2020. a
Wadhams, P.: Attenuation of swell by sea ice, J. Geophys. Res., 78, 3552–3563, https://doi.org/10.1029/jc078i018p03552, 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, https://doi.org/10.1029/jc093ic06p06799, 1988. a
Wang, R. and Shen, H. H.: Gravity waves propagating into an icecovered ocean: A viscoelastic model, J. Geophys. Res., 115, C06024, https://doi.org/10.1029/2009jc005591, 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 summerfall transition 2014, J. Geophys. Res.Oceans, 121, 1502–1525, https://doi.org/10.1002/2015jc011349, 2016. a
Williams, T. and Porter, R.: The effect of submergence on the scattering by the interface between two semiinfinite sheets, J. Fluid. Struct., 25, 777–793, https://doi.org/10.1016/j.jfluidstructs.2009.02.001, 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, https://doi.org/10.1016/j.ocemod.2013.05.010, 2013. a, b, c, d, e, f, g, h
Williams, T. D., Rampal, P., and Bouillon, S.: Wave–ice interactions in the neXtSIM seaice model, The Cryosphere, 11, 2117–2135, https://doi.org/10.5194/tc1121172017, 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, http://hdl.handle.net/10523/8154 (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, https://doi.org/10.1002/2015jc010770, 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, https://doi.org/10.1038/s4155801803503, 2018. a
 Abstract
 Introduction
 Preliminaries
 Methods
 Monochromatic forcing
 Polychromatic forcing
 Discussion and conclusions
 Appendix A: Postbreakup floe positioning
 Appendix B: Sensitivity to the choice of spectrum
 Appendix C: Histograms
 Appendix D: Kolmogorov–Smirnov statistic
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Preliminaries
 Methods
 Monochromatic forcing
 Polychromatic forcing
 Discussion and conclusions
 Appendix A: Postbreakup floe positioning
 Appendix B: Sensitivity to the choice of spectrum
 Appendix C: Histograms
 Appendix D: Kolmogorov–Smirnov statistic
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References