Modelling wave-induced sea ice break-up in the marginal ice zone

A model of ice floe break-up under ocean wave forcing in the marginal ice zone (MIZ) is proposed to investigate how floe size distribution (FSD) evolves under repeated wave break-up events. A three-dimensional linear model of ocean wave scattering by a finite array of compliant circular ice floes is coupled to a flexural failure model, which breaks a floe into two floes provided the two-dimensional stress field satisfies a break-up criterion. A closed-feedback loop algorithm is devised, which (i) solves the wave-scattering problem for a given FSD under time-harmonic plane wave forcing, (ii) computes the stress field in all the floes, (iii) fractures the floes satisfying the break-up criterion, and (iv) generates an updated FSD, initializing the geometry for the next iteration of the loop. The FSD after 50 break-up events is unimodal and near normal, or bimodal, suggesting waves alone do not govern the power law observed in some field studies. Multiple scattering is found to enhance break-up for long waves and thin ice, but to reduce break-up for short waves and thick ice. A break-up front marches forward in the latter regime, as wave-induced fracture weakens the ice cover, allowing waves to travel deeper into the MIZ.


Introduction
The Arctic marginal ice zone (MIZ) that separates open ocean from the interior pack ice is experiencing rapid changes as a result of high-latitude climate change. During summer, for example, its extent relative to the total sea ice area is expanding [1], suggesting an of sources and sinks that describe the effects of thermodynamics (melting and freezing), lead opening, ridging and fragmentation, noting the latter phenomenon is parametrized in a highly simplified manner (uniformly redistributed floe sizes). Horvat and Tziperman provide the most advanced approach, by considering the joint floe size and thickness distribution and processinformed parametrizations of the source and sink terms, particularly wave-induced floe fracture, which accounts for the strain field generated by a wave spectrum. Also note that Herman [20] showed that the Pareto distribution used to fit observational FSD data (discussed earlier) emerges as a stable solution to a stochastic model of FSD evolution based on the generalized Lotka-Volterra equation.
Here, the three-dimensional phase-resolving scattering model of wave energy attenuation in the MIZ reported by Montiel et al. [9] is enhanced with a floe break-up model, allowing us to investigate the two-way wave-MIZ coupling in an idealized setting. The MIZ is constructed as an array of circular elastic floes with prescribed FSD and the forcing field is approximated by a monochromatic plane wave. The solution to the wave interaction problem provides a full description of (i) the wave field throughout the MIZ and (ii) the bending experienced by each floe. The latter information is used to derive a measure of elastic deformation in each floe which, if larger than a critical value, results in floe fracture. The post-break-up updated FSD is then fed back into the geometrical description of the MIZ, leading to a new solution of the wave interaction, which is in turn used to approximate ice floe break-up. Running this feedback loop simulation a sufficient number of times, we reach a steady-state FSD, which depends on the ice and wave parameters. The main goals of this investigation are to (i) study the evolution of the FSD towards its steady state under repeated wave break-up events, (ii) determine the effect of multiple scattering on the steady-state FSD, and (iii) examine the dependence of the FSD on the ice and wave parameters of this model.
A key novel feature of the flexural failure model proposed here is that it accounts for the two-dimensional stress field defined over the surface of each deformed floe. It is based on the twodimensional Mohr-Coulomb (MC) stress criterion, which assesses mechanical failure from the combined level of tensile and compressive deformations at each point of the floe. The MC stress criterion has been used to estimate fracture in the elasto-brittle rheological sea ice model neXtSIM under horizontal deformations and wave-induced flexure [26]. The latter fracture model is a simplified one-dimensional version of the MC criterion used here, as wave-induced ice flexure is approximated using an elastic beam model of ice floe. This is in line with previous flexural failure models [22][23][24]27,29,30], in which deformation at each point of the beam is simply quantified by the curvature of its vertical displacement function. To our knowledge, the wave-induced break-up model of ice floes considered here is the first one to account for the additional spatial dimension.
We do not attempt to compare our model with experimental measurements in this investigation. Ice floe break-up by ocean waves in the MIZ has been reported via either in situ, e.g. [13,[30][31][32], or remote-sensing, e.g. [12,33,34], observations. These papers describe qualitative or quantitative changes in the FSD after a large wave event breaks up the MIZ. It is unclear, however, which physical processes have contributed to the observed changes in the FSD, as it is not possible to isolate the effect of waves. For this reason, we focus our analysis on gaining a theoretical understanding of how ocean waves may influence sea ice break-up in the MIZ and the associated FSD.

Preliminaries
We consider a three-dimensional seawater domain with constant finite depth h and infinite horizontal extent. Cartesian coordinates x = (x, y, z) are used to locate points in the water domain, such that the planes z = 0 and z = −h coincide with the unperturbed free surface and the flat impermeable seabed, respectively. The seawater is approximated as an inviscid and homogeneous incompressible fluid with density ρ 0 ≈ 1025 kg m −3 .
A finite array of N f compliant sea ice floes is assumed to be freely floating at the equilibrium surface of the water domain. Each floe is circular with uniform thickness D and Archimedean draught d = (ρ/ρ 0 )D, where ρ ≈ 922.5 kg m −3 is the density of sea ice. The radius of floe i is a i and its centre has coordinates in the horizontal plane (x i , y i ). The horizontal region of seawater covered by a floe is defined by for any i ∈ I, where I = {1, 2, . . . , N f }. We further denote their union as Ω = Ω 1 ∪ Ω 2 ∪ · · · ∪ Ω N f and the horizontal region covered by a free surface as Ω 0 = R 2 \ Ω. We consider time-harmonic perturbations in the water with prescribed radian frequency ω. Assuming the flow is irrotational, the velocity field in the water domain is expressed as (∇, ∂ z ) Re{(g/ iω)φ(x) e − iωt }, where ∇ = (∂ x , ∂ y ) and g ≈ 9.81 m s −1 is the acceleration due to gravity. The complex-valued potential field φ then satisfies Laplace's equation in the water domain The condition of no normal flow on the seabed yields, in addition, the Neumann boundary condition We assume that the perturbations in the water induce a flow characterized by a vertical displacement at the free surface that is small compared with the horizontal characteristic length of the flow. The linearized free surface boundary condition then takes the form We introduce a coupling between the vertical deformations experienced by each ice floe and the flow in the water domain. Horizontal motions of the floes are neglected, a valid assumption in the regime considered here [35]. We model the ice floe vertical deformations using the Kirchhoff-Love thin-elastic plate theory. This model is valid provided (i) ice floe diameters are large compared with the thickness and (ii) vertical deformations are small compared with the thickness. The boundary condition on the underside of the ice floes is then given by where F = ED 3 /12(1 − ν 2 ) is the flexural rigidity of the floe, which depends on the effective flexural modulus E [21] and Poisson's ratio ν. The values E ≈ 6 GPa and ν ≈ 0.3 are commonly used for sea ice. The requirement that each floe i ∈ I has no horizontal motion is written as and the free edge conditions of zero bending moment and zero vertical shear stress at the edge are, respectively, where (r i , θ i ) are the local polar coordinates with origin at the centre of floe i, defined by ( has also been introduced. An ambient flow in the water domain is prescribed with potential φ am (x) = ψ 0 (z) A am (τ ) e ik 0 (x cos τ +y sin τ ) dτ , function ψ 0 (z) = cosh k 0 (z + h)/ cosh k 0 h, where k 0 denotes the wavenumber of travelling waves in the open ocean (defined later).
We seek a solution of the boundary-value problem (2.2)-(2.8) of the form φ = φ am + φ S , where φ S is the potential of the scattered wave field due to the presence of the ice floes in response to the ambient wave potential φ am . In the far field, the scattered wave potential satisfies the Sommerfeld radiation condition √ r(∂ r − ik)φ S → 0 as r → ∞, (2.10) where r = x 2 + y 2 .
3. Wave-scattering model (a) Single floe scattering We decompose the potential in the exterior open water region adjacent to any floe i ∈ I (i.e. for in is the local incident wave potential generated by sources away from the floe and φ (i) sc is the scattered wave potential generated due to the presence of floe i. Standard cylindrical eigenfunction expansions are used to express these potentials (e.g. [36]) The potential φ ≡ φ (i) int in the interior region to floe i is expanded as in which the summation over m starts at −2 to account for the contribution of two vertical wave modes not present in the open water region. In (3.1), J n and H n denote the Bessel and Hankel functions of the first kind of order n, respectively. The eigenfunctions describing the fluid flow in the vertical direction in the open water and ice-covered regions are given by respectively. The quantities k m , m ≥ 0, are solutions of the open water dispersion relation We denote by k 0 the only positive real root of (3.2), and by k 1 , k 2 , k 3 , . . . the infinite number of purely imaginary roots with positive imaginary part ordered such that Im(k m+1 ) > Im(k m ) for all m ≥ 1. The real root k 0 is the wavenumber of horizontally travelling wave modes at the free surface of the open water region, while the imaginary roots k m , m ≥ 1, are associated with horizontally evanescent wave modes decaying exponentially faster from their source for increasing values of m. The quantities κ m , m ≥ −2, introduced in (3.1c) are solutions of the ice-covered dispersion relation The scaled elastic constant β = F/ρ 0 g has been introduced in (3.3). We denote by κ 0 the only positive real root of (3.3), which is the wavenumber associated with horizontally travelling wave modes at the water-ice interface. In addition, κ 1 , κ 2 , κ 3 , . . . designate the infinitely many imaginary roots with increasingly large imaginary parts associated with evanescent waves, and κ −2 and κ −1 denote the two remaining complex roots with positive imaginary part that are associated with damped travelling wave modes. A relationship exists between the unknown coefficients a m,n and c (i) m,n of the eigenfunction expansions given in (3.1), as a consequence of the boundary conditions prescribed on the surface r i = a i . In addition to the condition of no horizontal motions (2.6), continuity of fluid pressure and normal velocity is imposed at the interface between the open water and ice-covered regions, i.e.
A numerical solution is obtained by approximating the series expansions in (3.1) as partial sums, such that m ≤ M and |n| ≤ N, where M and N are convergence parameters. We apply the eigenfunction matching method (EMM) proposed by Montiel et al. [35]. It provides the following matrix relations between the coefficients: where a (i) , b (i) and c (i) are the column vectors containing the coefficients a

(b) Multiple scattering
We use a self-consistent method to resolve wave interactions with the array of N f ice floes, in which the incident field φ (i) in on each floe i ∈ I is given by the coherent superposition of the prescribed ambient field and the field scattered by all the other floes. This is expressed as sc for all i ∈ I. (3.6) This system of N f equations can be solved after (i) writing the truncated cylindrical eigenfunction expansion of φ am in the local coordinate system associated with floe i, (ii) applying Graf's addition theorem [37] to express φ (j) sc in the local coordinate system associated with floe i, and (iii) using the exterior DTM of floe i to express the expansion coefficients of φ (j) sc in terms of those of φ (i) in . This procedure yields a coupled system for the coefficients a (i) m,n , which can be inverted directly [35,38,39] or solved iteratively [40,41]. The local scattered field coefficients b Numerical issues discussed by Montiel et al. [9,42] arise when the number of floes becomes larger than O(100) and/or when the floes are closely spaced. As a remedy, Montiel et al. proposed an algorithm that combines the direct approach with a domain decomposition technique, referred to as the slab-clustering method (SCM), in order to resolve wave interactions with O(10 4 -10 5 ) floes. We only give a brief summary of the key steps of the SCM here and the reader is referred to the original papers [9,42] for further details.
(i) We cluster the array of floes into N s slab regions parallel to the y-axis, so that each slab q is bounded by x = ξ q−1 and x = ξ q , with ξ q−1 < ξ q , and contains the centre of N (q) floes. In each slab q, we apply the direct method summarized above to obtain the matrix mapping between the vectors a q containing the coefficients of the locally incident field φ of the ambient field and the field scattered by adjacent slabs expressed in the local polar coordinates of each floe in slab q.
(ii) We decompose the potential at each interface q is a field propagating or decaying in the positive/negative x-direction. Neglecting the vertical evanescent wave modes associated with the imaginary roots of (3.2), we approximate these components as The validity of this approximation was confirmed by Montiel et al. [9]. We have introduced the unknown amplitude functions A ± q (χ ) and the directional parameter χ ∈ Λ, where Λ is the integration contour which extends into the complex plane. It is defined by A value of χ ∈ Λ r corresponds to a plane wave travelling at angle χ with respect to the x-axis, while χ ∈ Λ ± i corresponds to an evanescent wave decaying exponentially with x. Such evanescent wave components are generated by wave sources of the form (3.1b) from floes present in the adjacent slabs. (iii) The amplitude functions A ± q−1 (χ ) and A ± q (χ ), respectively, defined on the left and right boundary of slab q, are related through the following integral scattering relationships: where R ± q (χ : τ ) and T ± q (χ : τ ) are the so-called reflection and transmission kernels of slab q. Semi-analytical expressions for these kernels were derived in [9]. They are obtained by combining (3.7) with mappings between cylindrical wave fields and plane wave fields. (iv) We approximate numerically the scattering relationships (3.9) by discretizing the truncated integration contourΛ =Λ − i ∪ Λ r ∪Λ + i , whereΛ ± i = ±π/2 ∓ i(0, δ) for some δ ≥ 0, sampling the amplitude and kernel functions at N Λ discrete χ and τ values, and integrating this equation numerically (composite trapezoidal rule). We then obtain the following matrix equations: where A ± q are column vectors containing the sampled values of A ± q (χ ) and R ± q and T ± q are square matrices containing sampled values of the reflection and transmission kernels and the quadrature weights. (v) We solve the set of 2N s matrix equations defined by (3.10) for the amplitude vectors A ± q , q = 1, . . . , N s , using the iterative S-matrix method of Ko & Sambles [43]. This requires initialization of the forcing amplitudes as The local scattered wave fields in each slab q are obtained after transforming the plane wave forcing fields φ + q−1 and φ − q into cylindrical regular wave fields with amplitudes contained in f q and successively applying (3.7) and (3.5).
Convergence of the numerical method described here depends on the truncation parameters M and N in approximating cylindrical series expansions (3.1), and δ and N Λ in approximating the plane wave expansions (3.8). For the computations carried out in this study, we set M = 0, N = 15, δ = 1.2 and N Λ = O(100) (depending on the wave frequency), allowing us to compute scattered wave coefficients with three-digit accuracy. The choice M = 0 simply means that we are ignoring the evanescent vertical wave modes. A comprehensive convergence analysis is conducted in [9] to justify these values.

Floe break-up criterion (a) Mohr-Coulomb stress
The thin elastic plate model of an ice floe considered here, with its underlying assumption of plane stress, allows us to write the Cauchy stress and strain tensors at each point as respectively, for any floe i ∈ I. The tensorial components with subscript r and θ are the normal stresses/strains in the radial and azimuthal direction, respectively, while the component with subscript rθ denotes the shear stress/strain. We can express the components of the strain tensor as [44] Note that a typographical error in the expression of the shear strain given in [44] has been corrected here. We have defined the vertical displacement w (i) ≡ w (i) (r i , θ i , t) of floe i, which can be related to the potential on the undersurface of the floe through the kinematic condition At each point of this undersurface and at a given time t, we then compute the components of the strain tensor from (4.2) after using (4.3) and (3.1c) to express the vertical displacement. Note that asymptotic formulae must be used for ε θ and ε rθ in the limit r i → 0 [45].
The components of the stress tensor S are related to those of the strain tensor E through We can diagonalize the stress tensor as S =ṼDṼ T , wherẽ The eigenvalues σ 1 ≡ σ 1 (r i , θ i , t) and σ 2 ≡ σ 2 (r i , θ i , t) of S are referred to as the principal stresses. The corresponding normalized eigenvectorsṽ 1 (r i , θ i , t) andṽ 2 (r i , θ i , t) are orthogonal and define the so-called principal directions for which the shear stress vanishes in the polar coordinate frame centred at (r i , θ i ). In the Cartesian frame with origin at the centre of floe i, the matrix of eigenvectors becomes We now derive a criterion for the break-up of an ice floe which incorporates some distinctive mechanical properties of sea ice. There exists a range of models for the mechanical failure of a large variety of materials under different types of loads ( [46], ch. 2). They take the form F (σ 1 , σ 2 ) ≥ 0, where F = 0 denotes a curve in the (σ 1 , σ 2 )-plane that corresponds to the onset of mechanical failure, and is often referred to as the yield curve. Here, we associate the 'yielding' of an ice floe with its fracture, which is reasonable for sea ice experiencing strain rates associated with wave periods of 5-20 s where plastic yield is negligible [47].
To our knowledge, no study has attempted to determine an appropriate model of waveinduced thin plate failure for sea ice, among the range of existing models in the literature of fracture mechanics. For this reason, we choose the simplest model of mechanical failure, applicable to both fracture and yield, and used for materials that exhibit significantly different values of tensile and compression strengths, defined as the maximum tensile and compression  stresses that sea ice can experience before fracturing, respectively. In the case of sea ice, the tensile strength σ t is typically an order of magnitude lower than the compression strength σ c [48]. The failure model we use here is referred to as the Mohr-Coulomb (MC) criterion [46]. The yield curve is defined by where ) as the MC stress at the point (r i , θ i ) and time t, and σ MC as the MC critical stress. The yield curve F = 0 is depicted in figure 1a. Neglecting the effect of fatigue, F is assumed to be stationary. We now define the criterion for floe break-up as follows: a floe i ∈ I fractures into two smaller floes if the condition is satisfied. We refer to σ (i) br as the potential break-up stress of floe i. The break-up criterion depends on the compressive and tensile strength of sea ice, which we estimate to be σ c = 3 MPa and σ t = 0.5 MPa, respectively. These values were chosen from empirical formulae reported by Timco & Weeks [48], assuming a brine volume fraction of approximately 0.04 and a strain rate of 2 × 10 −5 s −1 consistent with a loading at a mid-range wave period of 10 s. These values are chosen to be realistic for sea ice in the MIZ, as opposed to attempting to replicate a particular observed field condition, which is likely to be associated with much variability for these parameters. The sensitivity of the break-up simulations conducted here with respect to the sea ice mechanical properties is not within the scope of the present work. It is noted that flexural strength of sea ice may be a better approximation for σ t than the tensile strength used here for flexurally induced fracture. The value calculated from the empirical formula given in [48], however, is close to that for the tensile strength, so the distinction is actually of no practical importance.  We devise a sensitivity test to assess the potential for break-up of a single ice floe (i.e. N s = 1 and N (1) = 1) for a range of radii and wave periods. We prescribe a unidirectional plane wave with unit amplitude travelling in the positive x direction, by setting A am (τ ) = δ(τ ), where δ denotes the Dirac delta. We fix the water depth to h = 200 m. We assume, without loss of generality, that the floe has its centre coinciding with the origin (x, y) = (0, 0) of the horizontal Cartesian coordinate system. We compute the potential break-up stress σ br for wave periods T = 5-20 s, floe radii a = 5-500 m and floe thicknesses D = 1, 2 and 4 m. For all three thicknesses considered, we generally observe that the potential break-up stress increases rapidly with the floe size before plateauing for floe radii greater than a certain value. This behaviour is seen at all wave periods for D = 1 m and for periods greater than approximately 6 and 9 s for D = 2 and 4 m, respectively. For shorter waves, σ br oscillates with respect to a, suggesting resonances periodically induce large stress values, as the floe size approximately equals an integer multiple of the wavelength. We expect that similar oscillations would be observed for D = 1 m at wave periods less than 5 s.
According to our break-up criterion (4.10), fracture occurs for e ≥ 1, corresponding to the second contour in figure 2, i.e. the one separating the second and third darkest shades of blue. For all wave periods, it is seen that break-up occurs for all floe radii larger than a critical radius denoted by a MC and referred to as the MC radius. It is shown as a function of wave period in figure 3 for the three ice thicknesses considered (dashed lines). We observe that it reaches a minimum at T = 7 and 10 s for D = 2 and 4 m, respectively, while it seems to approach a minimum near T = 5 s for D = 1 m. These correspond to resonant frequencies, when the floe diameter approximately coincides with the open water wavelength and half the ice-covered wavelength.
In figure 3, we further indicate by solid lines the radius corresponding to the contour defined by e = 0.5 in figure 2. We denote it by a crit and refer to it as the critical break-up radius. Although break-up does not occur for these radii in the single floe scattering simulations conducted here with a unit amplitude plane incident wave, numerical experiments showed that multiple interacting floes may generate more energetic wave forcing with the ability to fracture smaller floes as a result of constructive interference. In all subsequent simulations, we set a crit = a crit (T) as the minimum floe radius below which break-up cannot occur, unless otherwise discussed. This approximation allows us to perform multiple scattering simulations at a manageable computational cost. It should be noted that our arbitrary choice of a crit is conservative in most  cases in the sense that the probability that a floe with a < a crit will fracture is small compared with the probability that a floe with a > a crit does not.

Break-up model
We The algorithm used for our break-up simulations is outlined as follows.
(i) Compute the exterior and interior DTMs associated with each floe radius a (p) , for p = 1, . . . , N r , using the method discussed in §3a, and store them (the size of each DTM is O(10)).
where t i and Ω (2) i as shown in figure 1b. We then define a (1) i and a (2) i as the radii of the discs with the area of regions Ω (1) i and Ω (2) i , respectively. Their expressions are Note that if either a (1) i or a (2) i is not equal to one of the radii a (l) , l = 1, . . . , N r , we round it to the nearest one. The break-up model described here should be seen as a new method to generate an FSD from repeated wave-induced floe break-up events. No attempt was made to replicate the fracture mechanism of ice floes in the MIZ with realistic shapes. From this perspective, the gross approximation of generating two circular floes from the break-up of a circular floe is acceptable, as we are only interested in the size of newly created floes as opposed to their shape. This approach is analogous to that of [19], in which floe size in the MIZ was measured by calculating the diameter of discs with the same area as that of the observed floes.
In step (v), we only include floes with radius a (l) ≥ a crit in V (FSD) s to generate the updated random array. Neglecting the influence of smaller floes on the break-up simulations improves the efficiency of the scattering computation in step (ii), noting that these small floes are still counted in the vector V (FSD) s+1 . It should be further noted that, although the break-up algorithm described here is intended to preserve the ice concentration (defined as the fraction of ice-covered ocean surface), rounding the radius of the two floes generated in step (iv) introduces a small change of concentration. However, simulations (not shown here) have revealed that these changes average out over a large number of break-up events, so the concentration actually remains quasi-constant.  The floe break-up algorithm described in §5 is then used to simulate the evolution of the FSD for N br = 50 break-up events, which we find is generally sufficient to reach a steady state. The ice concentration is set to 50%, so that the random array of floes generated after each break-up event is enclosed in a square region with side length √ 2πa max . The outputs of the break-up algorithm are the vectors V In figure 4, the floe size PDF is plotted after s = 5, 10, 20 and 50 break-up events for the wave periods T = 6 s (blue lines), 10 s (red lines) and 14 s (green lines) and the three ice thicknesses D = 1 m (figure 4a-d), 2 m (figure 4e-h) and 4 m (figure 4i-l). The first striking feature is that the distributions obtained from the break-up simulations clearly cannot be identified as power-law distributions, as the number of floes decreases to zero for smaller and smaller radii. Instead, the distributions after 50 break-up events look either unimodal or multimodal, as the distributions contain one or multiple maxima. Closely spaced successive maxima, such as those seen for T = 6 s and D = 2 m around a = 20 m, are likely to be an artefact of the floe radius sampling chosen for the simulations, so that they actually represent a single mode of the corresponding continuous distribution. We regard these oscillations as background noise. Accounting for this, we propose that all the PDFs are either unimodal or bimodal. Unimodality is observed for D = 1 m at all wave periods and D = 2 m for T = 6 s, while bimodality seems to manifest itself in all other cases. The second mode for D = 4 m at T = 6 and 10 s is located at a ≈ 55 m for both wave periods,   although the magnitude of its peak is only slightly larger than the background noise, so caution is recommended in interpreting them as modes. It should be further noted that the unimodal distributions are all positively skewed, suggesting they may be the superposition of two closely spaced unimodal distributions. The bimodality may be explained by the fact that each floe with sufficient stress breaks into two floes only, so that the repeated break-up of the ice cover results in two dominant floe sizes that correspond approximately to the break-up of the smallest floe that can fracture for a given thickness and wave period. We do not attempt to analyse the bimodal property of the PDFs further. For each parameter configuration considered in figure 4, we also observe a convergence of the PDFs with respect to the number of break-up events, as all distributions obtained for 20 events are almost identical to those obtained after 50 events, indicating that the break-up has ceased and a steady state has been reached after 20 events. The convergence seems to occur faster for longer waves and thicker floes, which is a consequence of fewer floe break-ups occurring for increasing values of these parameters.

Results
To understand better the convergence of the FSDs, figure 5 shows the evolution of the mean and standard deviation (s.d.) of the PDFs, and the number of floes per square kilometre through the 50 break-up events, for each wave period and thickness considered in figure 4. In the first few break-up events (i.e. up to s ≈ 5), we observe that the statistics of the distribution change exponentially fast, corresponding to the regime in which all the floes in the array fracture, so that the number of floes doubles after each event. Interestingly, the mean floe size decreases independently of the wave period and floe thickness in this regime, while the standard deviation consistently remains higher at T = 6 s than at the other periods for the three floe thicknesses considered, suggesting that shorter waves generate a broader FSD under intense wave-induced break-up as the ice cover fractures very quickly.
For s > 5, the statistics of the FSD quickly reach a steady-state regime. The transition between the exponential break-up and steady-state regimes is very sharp for T = 14 s, i.e. within five more break-up events, while it is longer for shorter waves, e.g. at T = 6 s for which it can take more  For D = 1 and 2 m, the mean and standard deviation smoothly vary with wave period with a general increasing trend for longer waves. This can also be explained by the fact that floe breakup diminishes for longer waves, which tends to enhance the formation of bimodal distributions. Note the minimum reached by the standard deviation at T = 13 s for D = 2 m. Inspection of the PDF obtained for this parameter configuration (not displayed here) shows that the bimodal shape of the distribution is not apparent, in contrast with T = 12 and 14 s for which it clearly exists. We could not further explain this feature.
For D = 4 m, the mean floe size reaches a minimum at T = 7 s. It should be noted that this is smaller than T ≈ 9 s for which the minimum of the critical radius a crit (T) discussed in figure 3 occurs for this thickness. This is likely to be a consequence of multiple wave scattering enhancing the stress field in the floes within the array at smaller wave periods, and break-up ensuing due to constructive interference. This further explains the need to take values of a crit for the break-up simulations smaller than those computed in §4b, as discussed at the beginning of this section.

(b) Array of floes
We now investigate how wave scattering by an array of floes influences the break-up process and associated evolution of its FSD. We consider the following initial configurations:  break up more, resulting in an increased computational cost required to solve the corresponding wave interaction problem with many floes). All the floes have initial radius a max = 200 m and the ice concentration is 50%, so that arrays in configurations (i) and (ii) have approximate horizontal extents of 1.5 × 5 km and 2.5 × 10 km, respectively. The unique floe radii and critical radii a crit for each thickness are the same as those chosen for the single floe break-up simulations in §6a. We perform the array break-up simulations for N br = 50 break-up events. As a result of the significantly higher computational cost associated with the array break-up simulations compared with the single floe break-up simulations, no ensemble averaging was performed here, so all the results in this subsection are obtained from a single random realization for each wave period and floe thickness. Additional random realizations performed on a few selected cases showed remarkable consistency of the resulting FSD, suggesting very little variability exists in the stochastic process described here.
The mean, standard deviation and number of floes per square kilometre obtained after 50 break-up events for each wave period and floe thickness considered are plotted in figure 7, where they are compared with the steady-state statistics of the corresponding single floe break-up simulations. For D = 1 m, the steady-state statistics obtained by breaking up the 10 × 3 array are very similar to those obtained from the break-up of a single floe, over the range of wave periods. It should be noted that the mean is consistently slightly smaller, while the number of floes is slightly larger, suggesting more break-up takes place for the array simulations, probably as a result of enhanced floe break-up due to constructive interference caused by multiple scattering. A similar observation can be made for D = 2 m thick ice in the range of wave periods T ≥ 7 s. For shorter waves, there is a clear deviation from the single floe break-up statistics, as the mean and standard deviation increase and the number of floes decreases as T decreases, which is the opposite trend to what happens for single floe break-up. This indicates that scattering is sufficiently dominant in this regime to prevent some floes from fracturing. The large increase in the standard deviation further suggests that the FSD spreads over a larger range of floe radii. Results obtained for D = 4 m reinforce our explanations of the role of multiple scattering in break-up through the large array, with the observations of two regimes, i.e. T < 10 s (short waves) and T ≥ 10 s (long waves). In the short-wave regime, scattering generates sufficient wave energy attenuation to prevent floe breakup at some level of penetration in the ice-covered domain. In the long-wave regime, scattering enhances floe break-up due to constructive interference of the wave fields radiated by the freely floating individual ice floes.
We seek more insight about the FSD obtained from the break-up of large arrays by plotting the PDFs of the distributions after 50 break-up events in figure 8. We can identify the two regimes discussed previously, i.e. enhanced break-up and reduced break-up, as a result of multiple wave scattering by a large array of floes. Based on our observations from figure 7, we associate enhanced   break-up with the PDFs depicted in figure 7a,b,c,e,f ,i, corresponding to cases for which scattering is not significant (i.e. long waves and/or thin ice). The PDFs show the presence of a larger number of small floes compared with the single floe break-up simulations, while the presence of larger floes decreases. Although this can be viewed as a shift of the PDF towards lower floe radii, the shape of the distribution also changes, particularly for thicknesses D = 2 and 4 m. In these cases, the bimodality observed in the distributions obtained from the break-up of a single floe is not a persistent feature of the PDFs associated with the break-up of an array, as the second mode (i.e. the mode corresponding to larger floes) is damped or removed. For D = 1 m, the PDFs are similar to those of the single floe break-up, being only slightly shifted towards smaller floe sizes. The second regime, characterized by reduced break-up in the array, corresponds to the PDFs shown in figure 7d,g, h. The large spread of these distributions discussed above is clearly observed, with an increased presence of smaller floes (compared with the single floe break-up), as is the case in the long-wave/thin ice regime, as well as larger floes. Important wave scattering occurring in the front slabs of the array prevents wave energy from causing break-up deeper in the array. In the extreme case D = 4 m and T = 6 s depicted in figure 7g, we can see the presence of floes with radius a = 200 m, which is the initial radius of all floes. Wave energy attenuation due to scattering is sufficiently strong in this case that some floes located deep enough in the array do not break.
To   The reader is reminded that the results presented in this section are based on a single random realization of the break-up simulations, so that an ensemble average may be necessary to resolve the small trends associated with the effect of penetration in the array on floe break-up. For the two short-wave cases (T = 6 s) considered in figure 8, we clearly observe the effect of scattering and associated wave energy attenuation in preventing floe break-up at some distance in the array. For D = 2 m (figure 9a-c), the FSD in the first row converges quickly to a steady state, while that in row 5 also converges but at a slower rate as the steady state seems to be approximately reached towards the end of the simulation and with a significantly smaller number of floes than in row 1. This suggests that wave attenuation has an effect in reducing break-up in the first few rows. The evolution of the FSD deeper into the array is different, as we observe break-up taking place during the first four events but then suddenly stopping. Break-up then resumes in rows 10 and 15 after eight and 26 events, respectively, while it does not in row 20. It is hypothesized that large floes in these rows initially fracture, as they do not require much energy to reach the critical break-up stress, as suggested in figure 2. After a few break-up events, the floes become too small to break under a wave field strongly attenuated by the front rows. As break-up in these front rows continues, however, wave attenuation due to scattering by smaller and smaller floes also decreases, so that more wave energy propagates further into the array with the result that floes are gradually broken up deeper and deeper in the ice field. In other words, we have shown in our model that wave-induced break-up has the capacity to reduce the structural integrity of the MIZ, enabling waves to travel further and cause break-up there, further weakening the ice cover. This positive feedback process is often used as a motivational concept for observational studies on wave-sea ice interactions. Note that a steady state is not reached after 50 break-up events in this simulation, so it is possible that break-up starts to resume in row 20 after the rows in front have experienced sufficient break-up.
In the case D = 4 m (figure 9g-i), break-up in rows 1 and 5 behaves similarly to that for D = 2 m, but rows 10 onwards do not experience break-up during the 50 events simulated here. Because thicker ice tends to increase the degree of wave attenuation, this is a plausible outcome of the simulation. Inspection of the evolution of the FSD in rows 6, 7, 8 and 9 (not shown here) indicates that break-up is ongoing after 50 events, so that rows 10 onwards may gradually start breaking after a larger number of break-up events.

(c) Wave spectrum
We now consider the break-up of the arrays used in the previous section under a unidirectional wave spectrum, in part acknowledging that the FSD observed in ice-covered oceans is the result of break-up by a sea state composed of a range of frequency and directional components. We do not seek to reproduce the break-up response to an observed random sea state in the present model, but instead we conduct a controlled numerical experiment to assess qualitatively how the FSD induced by an ad hoc wave spectrum differs from that obtained under monochromatic wave forcing. For this purpose, we sample wave components from a two-parameter Bretschneider spectrum defined by the spectral density function where T p is the peak wave period and H s is the significant wave height [49]. Sampling the period spectrum at integer wave periods between 5 and 15 s, the amplitude of the ambient field at each centre period of 1 s bandwidth is given by At this resolution, the incident wave amplitude is close to 1 m for each sampled period, being slightly larger around the peak period and slightly smaller away from it, and hence we expect break-up to occur for all periods, which is our intent in using this ad hoc model. Choosing a finer resolution would reduce the amplitude of each component and may not cause sufficient break-up to analyse the FSD and compare it with that obtained under monochromatic forcing. For our break-up simulations, we select one wave period T ∈ [5,15] randomly at each event and compute the break-up induced by a unidirectional plane wave of period T with amplitude given by (6.2). An ensemble of 10 random realizations of the break-up simulation is computed for each ice thickness. Although this approach is potentially different from break-up by a wave field composed of multiple frequencies, it is conjectured that the randomization of both wave period and array in conjunction with ensemble averaging provides a legitimate approximation. We set the parameters of the spectrum to T p = 10 s and H s = 2 m, which corresponds to a typical swell observed in the Southern Ocean [5]. Figure 10 shows the evolution of the mean, standard deviation and number of floes per square kilometre of the FSD for different rows of the array and for the entire array (figure 10a-c,e-g,i-k). We observe a clear convergence trend towards a steady state for all rows of the array and for each thickness. The front rows consistently converge faster than the back ones and to a lower mean and a higher number of floes, so that the degree of break-up decreases with the level of penetration in the array. Interestingly, the evolution of the mean and number of floes of the FSD for the entire array almost coincides with those associated with the middle row of the array. This suggests a simple linear dependence of the FSD statistics with respect to the row number, which will be demonstrated later. The standard deviation of the FSD for the entire array is consistently at least as large as that of the last row, which in turn is the largest of all the rows. This reflects the larger spread of floe sizes on a large scale (entire array) compared with the local scale (each row).   We further display the PDF of the FSD for all cases discussed above in figure 10d,h,l. Observe first that the distributions look nearly normal. This is probably the result of effectively averaging over the wave periods. The bimodality seen in previous distributions does not completely disappear, however, as a shoulder-type feature can be seen on the right of the peak (i.e. for large radii), particularly for D = 1 and 2 m. This suggests the existence of a second mode for large radii but with a small peak.
The row dependence of the mean, standard deviation and number of floes after 50 break-up events is analysed further in figure 11. As indicated above, we observe a linear increase in the mean of the FSD and a linear decrease in the number of floes with respect to row number. This allows us to estimate the rate of change of mean floe size with respect to penetration in the array, which is important for large-scale modelling studies of the MIZ. Fitting a straight line through the curves generated and extracting the slope, we find that the mean floe radius increases at a

Conclusion
A new model of ice floe break-up under ocean wave forcing in the MIZ has been developed. It combines the time-harmonic multiple scattering theory for a finite array of floating elastic discs proposed by [9] with a parametrization of flexural failure causing an ice floe to fracture into two floes, provided that the stress field satisfies a particular break-up criterion, the so-called MC criterion. We derived a quantity, referred to as the MC stress, that uniquely defines the level of stress at each point of the surface of an ice floe, allowing us to test simply if break-up is expected to take place at any point. A numerical experiment was then conducted to analyse the MC stress experienced by a single ice floe under a unit amplitude unidirectional wave forcing and determine the regime in which we expect floe break-up. It was found that (i) a minimum floe diameter exists for each thickness below which break-up cannot occur, and (ii) this critical diameter depends on the wave period in a way that it reaches a minimum at a resonant wave period, for which the floe diameter is approximately equal to the open water wavelength and half the ice-covered wavelength.
A closed-loop feedback algorithm has been proposed to model the evolution of the FSD in the MIZ under a sustained wave event. Each loop consists of (i) computing the MC stress in all the floes, (ii) breaking up each floe satisfying the MC criterion, and (iii) generating a new array of floes from the updated FSD, which is then used as the geometry of the wave interaction problem in the next loop. We conducted a number of numerical experiments to determine the evolution of the FSD towards a steady state through 50 break-up events (i.e. loops), for different wave and ice configurations. We simulated the break-up of (i) a single large floe for a monochromatic unitamplitude plane wave forcing, (ii) an array of large floes for the same monochromatic forcing, and (iii) an array of large floes for a Bretschneider spectrum forcing. Key findings are summarized below.
(i) Break-up of a single large floe causes the emergence of a bimodal FSD for most wave and ice parameters considered. Larger values of wave period and ice thickness correspond to FSDs with larger floe sizes and more separated peaks of the associated bimodal distribution. The convergence of the FSD towards its steady state under repeated breakup events is very quick for long waves and slower for short waves. Increasing values of ice thickness also tends to decrease the rate of convergence of the FSD, suggesting that multiple wave scattering within the array of broken floes enhances break-up and therefore influences the steady-state FSD. (ii) Break-up of an array of large floes under monochromatic forcing provides additional insight into the effect of multiple scattering on the FSD. First, the bimodality of the FSD observed for the single floe break-up simulations is either damped or removed, as the population of large floes associated with the second mode is consistently redistributed to smaller floes. Second, results of our simulations indicate that two scattering regimes exist, i.e. long waves/thin ice and short waves/thick ice. In the former regime, waveinduced ice break-up is enhanced compared with the single floe break-up, which is likely to be a consequence of constructive interference of wave fields radiated by the individual floes. In the second regime, multiple scattering causes sufficient wave energy attenuation through the array to prevent some ice floes from fracturing, resulting in broader FSDs compared with those obtained from the repeated break-up of a single floe. array with a small upward trend with distance from the ice edge. In the short-wave regime, we observe the positive feedback between wave-induced ice break-up and iceinduced wave attenuation. Break-up originally only takes place in the front rows as waves are attenuated by scattering deeper in the array. After sufficient break-up, waves are less attenuated, carry energy at a higher level of penetration and cause break-up there. The overall effect is the observation of a break-up front marching forward in the MIZ as the structural integrity of the ice cover reduces under wave-induced break-up. (iv) Break-up of an array under a Bretschneider spectrum forcing generates near-normal FSDs for all ice thicknesses, which is likely to be a consequence of averaging with respect to the wave period. As opposed to the simulations with monochromatic forcing, breakup decreases with the level of penetration in the array, such that the mean floe size increases linearly with distance from the ice edge. The rate of floe size increase is larger for thicker ice.
An important outcome of our investigation is that no power-law FSD was generated from the simulations. Our model consistently predicts that the number of floes decreases to zero for smaller floe sizes in an MIZ with uniform ice thickness. This results from the fact that small floes are less prone to elastic deformations than large floes, i.e. they behave similarly to a rigid body. Analysis of the data generated in §4b shows that the wave energy required to cause break-up in small floes increases exponentially fast as floe size decreases below a critical size, so that flexural failure by ocean waves is unlikely to be responsible for the observed increasing number of small floes in the MIZ [19], acknowledging that small floes with size O(1 m) cannot be resolved accurately by current observation techniques, so that it is unclear whether the power law extends to very small scales. In any case, our analysis suggests that wave-induced floe break-up alone does not create or preserve the observed power-law FSD or the apparent increase in small floe numbers, but may partially contribute to it. Other processes acting on longer time and/or length scales (e.g. thermodynamics, internal stress or collisions) must be considered to explain this feature of the FSD. Although recent modelling work has shown the emergence of a power-law FSD [27,28,50], it is still unclear how much each process contributes to the observed result. It should also be noted that we attempted to fit a power-law curve in the large floe regime of the FSDs obtained from our simulations, i.e. for radii larger than the peak radius, but the limited extent of this regime (i.e. spanning less than one order of magnitude in floe radius) did not allow us to obtain statistically significant results.
Although our findings provide much theoretical understanding of the wave-induced ice break-up process in the MIZ, the underlying model was constructed based on a number of simplifying assumptions, which may influence certain results. Specifically, the FSD resulting from the break-up of a more heterogeneous ice cover, i.e. governed by a floe shape and ice thickness distribution and possibly different ice types, under a realistic multidirectional random sea state, may spread the FSD over a wider range of floe sizes to the point where a power law could be fitted to a portion of the curve. In addition, the validity of our break-up model, which involves breaking a circular floe into two circular floes and time-harmonic wave forcing, is unclear and requires further investigation. It would be difficult to relax these assumptions in the context of the three-dimensional wave-scattering model considered here. We may envisage a simpler two-dimensional model, however, in which the one-dimensional ice cover is initialized as a semiinfinite beam and forced by a transient incident wave generated from a frequency spectrum. More theoretical investigation into this physical process is needed to be able to provide large-scale ice-ocean wave forecasting models with a quantitative parametrization of the two-way coupling between ocean waves and sea ice. Monitoring the evolution of the FSD in an area of ice-covered ocean during a wave break-up event either in the field or in a controlled laboratory setting, with the ability to resolve small floes accurately, would be needed to provide a clearer picture of the complicated processes investigated here. Data accessibility. The Matlab codes and data files created for the numerical investigation reported here can be accessed through http://www.maths.otago.ac.nz/files/icebreakup/Data_breakup.zip.