Wave-ice interactions in the neXtSIM sea-ice model

In this paper we describe a waves-in-ice model which calculates ice breakage and the wave radiation stress (WRS) that is coupled to the new sea-ice model neXtSIM, which is based on the Elasto-Brittle (EB) rheology. We highlight some numerical issues involved in the coupling, and investigate the impact of the WRS, and of modifying the EB rheology to lower the stiffness of the ice in the area where the ice has broken up (the marginal ice zone, or MIZ). In experiments in the absence of wind, we find that wind waves can produce noticeable movement of the ice edge in loose ice 5 (concentration around 70%) — up to 36 km, depending on the material parameters of the ice that are used, and the dynamical model used for the broken ice. The ice edge position is unaffected by the WRS if the initial concentration is higher (&0.9). Swell waves (monochromatic waves with low frequency) do not affect the ice edge location (even for loose ice), as they are attenuated much less than the higher frequency components of a wind wave spectrum, and so consequently produce a much lower WRS (by about an order of magnitude at least). 10 In the presence of wind, we find that the wind stress dominates the WRS, which while large near the ice edge, decays exponentially away from it. This is in contrast to the wind stress which is applied over a much larger ice area. In this case (when wind is present) the dynamical model for the MIZ has more impact than the WRS, although that effect too is relatively modest. When the stiffness in the MIZ is lowered due to ice breakage, we find that on-ice winds produce more compression in the MIZ than in the pack, while off-ice winds can cause the MIZ to be separated from the pack ice. 15


Introduction
Wave-ice interactions have received a great deal of attention in recent years (e.g.Dumont et al., 2011;Kohout et al., 2014;Ardhuin et al., 2016Ardhuin et al., , 2017)), with progress in both modelling and measuring (particularly via Synthetic Aperture Radar imagery, or SAR) of waves in ice.To a large extent, this is due to climate change, with a series of record lows in both minimum and maximum Arctic sea-ice extents in the last decade (e.g.Meier, 2017).Specifically, large parts of the Arctic are becoming and are expected to become even more accessible for resource exploitation and shipping in the summer, whereas 10 years ago they weren't (e.g.Stephenson et al., 2011).Associated with this low sea-ice extent is an increased open-water fetch available for wave generation which means there are potentially more large wave events in the Arctic in summer (e.g. in the Beaufort Sea in summer 2012; Thomson and Rogers, 2014).As well as being dangerous for shipping in themselves, large waves also increase the amount of ice breakage in the marginal ice zone (MIZ), creating an extra hazard as small floes could potentially be thrown onto a ship deck for example.
Closely connected to waves in ice, but with other controlling factors apart from waves, is the concept of floe size distribution (FSD e.g.Toyota et al., 2011;Herman, 2010).This can influence both the dynamics and thermodynamics of the ice, ocean and atmosphere in the MIZ.For example, it affects sea-ice rheology (Herman, 2012;Feltham, 2005) and can increase wind/ocean drag and consequently increase the stresses applied to the ice.It can also enhance lateral melting in summer (Horvat et al., 2016;Steele, 1992).Horvat et al. (2016) showed that increased horizontal salinity gradients at the floe edges produced eddies which allowed warm water to travel under the ice floes and enhance the melting from the edges.This was true even for large floes (∼ 1 km), when the lateral-to-horizontal-surface-area ratio is quite small.(Previously, this ratio was used to compute results which indicated lateral melting was unimportant for floes larger than ∼ 100 m; Steele, 1992.)Models for full numerical FSD's (Zhang et al., 2016), where a histogram of floe size bins can evolve in time, as well as joint ice thickness and floe size distributions have been proposed (Horvat and Tziperman, 2015).In the latter model, each thickness category can have its own FSD.More parametric approaches have also been used (Dumont et al., 2011;Williams et al., 2013a;Bennetts et al., 2017).
On the sea-ice modelling side, there has been a lot of progress in making sea-ice dynamics more realistic, especially in the Arctic pack.Rampal et al. (2016) presented a validation of the neXt-generation Sea Ice Model neXtSIM, looking at sea-ice area and extent, sea-ice drift, and the spatial scaling of sea-ice deformation derived from SAR (see also Bouillon and Rampal, 2015b).The dynamical core of neXtSIM is the EB sea-ice rheology, which is a thin elastic plate model with stresses constrained by a Mohr-Coulomb failure envelope.If stresses become too large and leave this envelope in a grid cell, the ice stiffness inside that cell is reduced (in practice a parameter called the damage is increased) in order to bring the stresses back onto the failure envelope (Rampal et al., 2016, for more details).When one cell is highly damaged, the likelihood for the surrounding cells to also become damaged is increased, leading to the rapid (i.e. after a few sea-ice-model time steps) emergence of very localised lines of damaged cells where sea ice can deform almost freely.These lines of concentrated damage can accommodate large deformation (i.e., opening, ridging and shearing) in a way that is similar to the so-called linear kinematic features that are observed from satellites (Kwok, 2001).
In this paper we demonstrate the coupling of a waves-in-ice model (WIM) to neXtSIM in an idealised domain.The physical effects included in the coupling are the break-up of ice by waves, the wave radiation stress (WRS), and an additional (optional) feedback to the sea-ice model where the ice stiffness is reduced where the ice is broken (in the MIZ).We conduct experiments with waves by themselves to see the impact of the WRS on the ice edge location, and also with wind to see the relative importance of the wind stress and the WRS.In addition, we do some simulations to see the particular effects of the rheological change.
We also highlight some general numerical issues involved with coupling wave models and sea-ice models on different grids.In addition, we do some theoretical reformulations of the WIM to put the ice break-up model in the context of Mohr-Coulomb failure, and do some sensitivity tests of the sensitivity of the MIZ width to the Young's modulus in particular, as well as the small-scale "cohesion" parameter in the WIM breaking model.Its response to the Young's modulus was previously uninvestigated.
2 Sea-ice model

Evolution equations
The ice is modelled as a thin elastic plate (e.g.Fung, 1965, §16.8) with constitutive relation or in full: where σ ij and ε ij (i, j = 1, 2) are respectively the stress and strain tensors, ν is Poisson's ratio and Y * is the effective Young's modulus (depending on the concentration c and the damage d), given by where C is the compactness parameter, and Y 0 is the Young's modulus of fully-compacted, undamaged ice.
The momentum balance equation we will use is the following: here ρ i , h and u are the density, actual thickness, velocity and internal stress tensor of the ice (respectively), ∇ = (∂ x , ∂ y ) T is the horizontal gradient, and τ o and τ a are the applied stresses by the ocean and the atmosphere (respectively).These latter stresses come from quadratic drag laws.Note that we neglect the Coriolis force and the gravitational force due to the slope of the ocean surface because of our idealised domain.Also appearing in (4) are the wave radiation stress (WRS), τ w,i , and the term involving P , which is a strictly positive pressure that provides a resistance to compaction and ridging (i.e., it is only activated when the divergence ∇ • u < 0): where P * is the pressure parameter, and εmin = (0.01/86400) s −1 is the minimum divergence rate.If the ice becomes very damaged, and loses its stiffness, this term prevents the ice from piling up and becoming too thick.As a default, we use the standard value of P * = 12 kPa, as suggested by Thorndike et al. (1975), but we will test the sensitivity of our results to C (see §5.3).C = 20 is commonly used in the standard sea-ice models using a Viscous Plastic (VP) rheology, so the pressure drops by a factor of about 55 when the open water fraction increases from 0 to 20%.So, for example, increasing C to 40 means the open water fraction only needs to be 10% for the pressure to reduce by 55.
We also have equations for evolution of any conserved quantity φ: φ could be concentration (c, also requiring c ≤ 1), volume (ch) or variables relating to the damage (retrieved from (1 − d) −1 ).
The terms S φ are thermodynamic source/sink terms which are switched off for this paper, since the simulations are in an idealised setting and run for short durations.In an Eulerian frame of reference, but since we work in a Lagrangian frame the relationship is simply Dφ/Dt = dφ/dt.The −φ(∇ • u) term represents the conserved quantity decreasing if the divergence is positive e.g. if a triangle in the finite element mesh increased in area then φ should drop in that triangle.
Like Williams et al. (2013a), we will parameterise the floe size distribution in terms of the maximum floe size, D max (see §3.3), which we wish to advect like a tracer: D(D max )/Dt = 0.In the Lagrangian framework, advection is usually exact, unless a local remeshing is required.This happens if the triangles of the mesh become too deformed, and requires (local) interpolation of the advected variable.Details on the remeshing procedure in the neXtSIM model can be found in Rampal et al. (2016).
Additional (global) interpolation is required to obtain D max on the fixed grid of the WIM (see §4).We found that transporting and interpolating D max itself led to some errors, which were reduced by transporting an auxiliary variable or to progress from (neXtSIM) time step n to n + 1, N floes should change according to floes /c (n) , and being interpolated when either regridding or communication with the WIM is required.
The evolution of stress and damage from time step n to n + 1 is done via an intermediate stress calculation: where Φ d is a thermodynamic source term (again not used here), while Ψ (0 < Ψ ≤ 1) is a factor determined from the position of the stress vector relative to the Mohr-Coulomb failure envelope, described in section 2.3.There is no continuous version of (9) since fracturing is an extremely rapid process, well below our typical time step ∆t.

Uncoupled neXtSIM simulation
Since the damage variable d is probably unfamiliar to most readers, we include here an example simulation illustrating its main role in the EB rheology.Figure 1 shows four fields after a 2-day simulation.The wind stress plotted is calculated from the quadratic drag law where ρ a = 1.3 kg m −3 and u a are the density and 10-m-velocity of the air, while C d,a = 7.6 × 10 −3 is the drag coefficient of the wind on the ice.The gradient in the wind stress comes from the differences in relative velocity.We have plotted this stress as a reference for when we discuss the WRS.Initially, the concentration was relatively low, so the internal stress was also low (see the formulae for Y * and P in equations 3, 5), meaning the ice was almost in free drift, being compressed against the right hand boundary.As the concentration 5 increased, the internal stress increased causing it to fail (increase d) in localised regions.Comparing the damage with the concentration and thickness, it can be seen that the regions of high compression and thickening correspond to the regions where the damage is highest.This is the usual role (without waves) of the damage -to produce localised deformation and features such as thicker regions (under shearing or convergent conditions, such as in the current simulation) and leads (under shearing or divergent conditions).We note here that the initial combination of c = 0.7 (loose ice) with no damage is not inconsistent since the damage only increases if the concentration is high, although the reasons it is usually initialised to zero are: (i) for simplicity and (ii) since it is not an observable variable.It then evolves with the other variables in response to the applied forcings.

Mohr-Coulomb failure
Let σ 1 and σ 2 be the principal stresses, with compressions corresponding to positive stresses.Then a stress state is within the Mohr-Coulomb failure envelope if the conditions are satisfied (Schulson et al., 2006;Dansereau et al., 2016;Rampal et al., 2016), where and τ 0 is the cohesion, and µ is the internal friction coefficient.See Figure 2(a) for some example envelopes (τ 0 = 629 and 989 kPa).The lines σ 2 = σ c + qσ 1 and σ 1 = σ c + qσ 2 in the space of principal stresses (Schulson et al., 2006;Dansereau et al., 2016) correspond to the lines |τ | − µσ N = τ 0 , so the material fails when the applied shear force |τ | reaches the sum of the frictional force inside the material (µσ N ) and the cohesion of the material (τ 0 ).Now (Schulson et al., 2006), where ϑ is the angle between the maximum principal stress (taken as the most compressive stress), σ 1 , and the failure plane.This reaches its maximum value when tan(2ϑ) = 1/µ, so if µ = 0.7, the failure plane is oriented at about 27.5 • from the direction of σ 1 .Equation ( 12) also lets us derive the expressions for q and σ c .
The conditions (11b) are less certain since there are fewer measurements in pure tension or compression.In particular, extending the Coulomb branches into the third quadrant in principal stress space (see Figure 2 of Dansereau et al., 2016, who instead apply tensile failure criteria σ 1 , σ 2 ≥ −σ c /q) could be seen as theoretically suspect (since there should be no friction under tension), but the observations of Weiss et al. (2007; see Figure 2) seem to support this approach.In practice, using σ N ≥ σ N,min or σ 1 , σ 2 ≥ −σ c /q was found to make little difference to large-scale simulations.Similarly, σ N,max is set large enough that it is not reached in simulations, which is reasonable since few examples of large biaxial compressive stresses have been observed (Weiss et al., 2007).Note that Dansereau et al. (2016) chose not to close the failure envelope at all for this same reason.
Returning to (9), if σ is outside the envelope it is scaled back onto the nearest branch of the envelope by setting σ (n+1) = Ψσ , where Ψ < 1.This ensures that the stress always remains within the envelope, but the damage d is increased if this happens.Otherwise, if σ is inside the envelope, Ψ = 1 and the damage is unchanged.

Scaling of the Mohr-Coulomb envelope
Mohr-Coulomb envelopes have been observed on many different scales in rock mechanics, and has also been seen in ice.The parameter µ controls the orientation of fractures that form, while the cohesion sets the sizes of the stresses which cause any fractures, and so is more influential.
This property should scale as , where L c is the size of the defects, or "stress concentrators" (Weiss, 2013, §4.2).
Put in another way, where the additional indices 0 or 1 correspond to different scales on which fracture is occurring.sizes deduced from these envelopes, using the scaling law ( 13).(These defect sizes, or sizes of stress concentrators, are only meant to give an idea of the relative sizes compared to those corresponding to the second cohesion value which is approximated to be around 1 m which is of the same order as the ice thickness.The first defect size is of the same order as the grain size -the grains measured in the sample were columns of diameter 3.9 mm and length 1 cm.)For some additional context, we also give the value used in the reference simulation of Bouillon and Rampal (2015a).This large-scale cohesion is in contrast to our small-scale cohesion (Lc ∼ 1 cm), which we use to determine if single ice floes will fracture due to wave flexure.
Note that these values do not necessarily correspond to the breaking stress of ice since the measurements are not exactly taken at the point of fracture.The lab measurement (uni-axial compression test) should be closer since we know the ice did actually break and the scale of the measurement; the in-situ measurements are certainly underestimations since the ice did not break, and in fact the value of 1 kPa was derived from a 3-day subset of the time series which was bounded by the envelope with cohesion 40 kPa.That is, the lower in-situ value corresponds to more remote fracturing, or fracturing over a larger scale.
In their presentation of the dynamical core of the neXtSIM model (using a resolution of approximately 10 km), Bouillon and Rampal (2015a) found that the model was quite sensitive to the cohesion value when varied between 0.5 kPa and 8 kPa.
However, the results for τ L 0 = 8 kPa (the superscript 'L' here indicates it is the large-scale cohesion, as opposed to the smallscale one discussed below) and τ L 0 = 4 kPa were similar.In the follow-up paper to the aforementioned one, Rampal et al. ( 2016) used τ L 0 = 8 kPa, or L c ≈ 25 m.This gave good agreement with the deformation-scaling statistics.
In the simulations done in this paper we will use a model resolution of 4 km, so we will test a range of cohesions from 4-13 kPa to be somewhat consistent with the above choice.Also, we will discuss the ice breakage by waves (below in §3.4.1) in terms of Mohr-Coulomb failure, and define an additional small-scale cohesion τ S 0 and defect scale L c for the breaking criterion we settle on in section 3.4.2.
3 Waves-in-ice model

Attenuation
The amount of attenuation that waves in ice experience is the main factor in determining the amount of momentum transferred to the ice.However, definitive confirmation of any particular physical models for this is still lacking.Meylan et al. (2014) came up with an empirical formula fitted to Antarctic attenuation from the experiments reported by Kohout et al. (2014).Ardhuin et al. (2016) compared the creep model of Wadhams (1973) (also see Tolman et al., 2016, §2.4) with drifting buoy data from within the ice, with some success in the timing of the peaks in wave heights.Other theoretical models that have been used are a viscoelastic attenuation model (Wang and Shen, 2010), and "localisation" predicted by 1D multiple scattering models (Kohout and Meylan, 2008;Bennetts and Squire, 2012).In the wave scattering context, localisation refers to how these models predict exponential decay of waves as they travel into the ice.Or in other words, the wave energy is localised in the vicinity of the ice edge.Doble and Bidlot (2013) used the model of Kohout and Meylan (2008) in Antarctic simulations using WAM, while Williams et al. (2013a) used a theoretical result from Bennetts and Squire (2012) to investigate break-up by waves.Tolman et al. (2016, §2.4) give a full summary of waves-in-ice parameterisations implemented in Wavewatch III.
Our attenuation model is essentially model B from Williams et al. (2013a), slightly modified to allow Young's modulus to be varied.It has a scattering component determined from the expected number of floes per unit length, and a dissipative component coming from the drag model of Robinson and Palmer (1990)  As stated above, the choice of attenuation model is crucial in determining the wave radiation stress, yet physical mechanisms are still relatively uncertain.However, we can still calculate the response of the ice to waves attenuated by our model, and make conclusions which should still hold for similar ranges of the WRS.

Energy transport
A general formulation for wave energy transport is where C g = c g (cos θ, sin θ) T is the group velocity vector, c g = dω/dk, ω is the radial frequency, k is the wavenumber, and E is the spectral density function (SDF) of the variance of the wave elevation η: the SDF of the time-averaged energy is E = ρ w gE, where ρ w is the water density and g the acceleration due to gravity.
We neglect the terms S in and S nl , which represent wind generation and non-linear energy transfer between frequencies and directions (respectively).The term S nl moves energy from high frequencies to lower ones, and becomes more significant if E is larger.For example, Kohout et al. (2014) described a storm event off Antarctica (with approximate latitude 61 • S and longitude 125 • E) where the significant wave height was measured to decay linearly with distance into the ice, whereas it decayed exponentially during calmer periods.Li et al. (2015) attributed this to the effect of S nl , and the fact that lower frequencies are attenuated less than higher ones.Thus we need to remember that our results could change (e.g.waves could induce ice breakage further from the edge) if our wave forcing becomes very large.In particular, the WRS may also persist further than predicted with our linear model -however, it would also have a smaller size since the longer waves are attenuated less.
The scattering kernel K distributes energy from the incident wave among the other directions and is discussed further in the next section.Various authors (e.g.Perrie and Hu, 1996;Masson and LeBlond, 1989) have used the solution for a rigid circular floating disc to deduce an expression for K; Meylan et al. (1997) extended this to make the disc elastic, and this solution was also used by Zhao and Shen (2016); Ardhuin et al. (2016) used the simpler kernel K = α scat /(2π) to distribute the incident energy uniformly in all directions.However, due to the fact that these models conserve energy, i.e. 2π 0 the operator L scat has some zero eigenvalues.(This is most easily seen by considering the discretised version of (17) -i.e.
considering only a finite number of directions -which would state that all the columns of the matrix representing L scat add to zero.Thus the rows are linearly dependant and the matrix will have at least one zero eigenvalue.)This usually means that the solution E of (15) will usually not decay exponentially into the ice (in the absence of dissipation).(This decay depends on the eigenvector(s) corresponding to the zero eigenvalue, of course, but in general they are such that E does not decay into the ice.)As a result, the results of Ardhuin et al. (2016) which included scattering in this way were quite unrepresentative of phase-resolving multiple-scattering models such as those of Kohout and Meylan (2008) and Bennetts and Squire (2012).
Consequently, we will use K = 0 and not conserve energy, since we think that it is preferable to preserve the localisation predicted by the scattering models.

Floe size distribution
We use a parametric form of the FSD.We initially require that D max ≥ D min and that large floes (> 200 m) have a uniform floe size distribution -i.e.p(D|D max > 200 m) = δ(D − 200 m).This latter assumption is somewhat vestigial but was related to the fact that wavelengths that do breaking in the ice are usually less than about 400 m.The rest of our approximation is similar to the FSD used by Dumont et al. (2011), which was based on the renormalisation group (RG) approach to the same problem, used by Toyota et al. (2011).However, this formula made the mean floe size a discontinuous function of the maximum floe size, so we have modified it to a continuous (as opposed to discrete) FSD -a power-law-type probability density function p(D) truncated at D = D max , but with the same exponent as before: where γ = 2+log f / log ξ, f is the fragility in the RG formulation of Toyota et al. (2011), and ξ 2 is the number of pieces formed during each successive break-up in the same RG formulation.We use D min = 20 m, f = 0.9 and ξ = 2, making γ ≈ 1.84.
Results for the MIZ width (not shown) with the RG approach are similar to those with the FSD ( 18), but the momentum flux is less smooth, which could cause numerical problems.We recognise that both parameterisations are completely arbitrary, and that numerical histograms (e.g. as used by Horvat and Tziperman, 2015) are preferable in terms of being able to let the wave spectrum try to produce the FSD naturally.(They also let other factors influence the FSD more easily).However, the FSD itself is not the focus of this current paper, and these alternative models are quite costly and not trivial to implement, so we do not try them out here.modulus for an ice floe.However, for waves we are interested in the stresses that are induced by a vertical displacement η.The stresses are assumed to be confined to the horizontal plane and varying linearly with the vertical coordinate z = x 3 (z = 0 is the middle of the plate, and (Fung, 1965, §16.9).We then have the following results for stresses and strains where x 1 = x and x 2 = y.For a plane wave (travelling in the x direction with amplitude A) in a thin elastic plate, η = A cos(kx − ωt), ε 11 = k 2 zη, ε 22 = ε 12 = σ 12 = 0, and so the only non-trivial stresses are given by where σ 1 and σ 2 are the principal stresses in the horizontal plane.This meets the upper Mohr-Coulomb branch when ).Note that here the shape of the tip of the failure envelope makes a difference, since a pure tensile failure criterion would increase the lower limit on σ 1 to −σ c /q ≈ −1.04τ S 0 (which would be reached at smaller wave amplitudes).However, given the uncertainty about the failure envelope under pure tension and high compression, and so that our small-and large-scale envelopes have the same shape, we use (11) for wave failure also.

Breaking criterion
The maximum strains are produced when z = ±h/2 (at the upper and lower surfaces of the ice), and so for a plane wave for Barents Sea sea ice.)When considering flexural strength measurements, however, it is useful to remember how they are obtained.In a cantilever situation, an ice beam is subjected to a force F c at one end until it breaks at the other.The force is then converted to a stress in order to remove the effects of the beam dimensions according to the formula (Frederking and Svec, 1985), where L and b are the length and breadth of the beam respectively.(Similar formulae exist for three-and four-point-bending tests.)This conversion assumes that the beam can be modelled as an Euler-Bernoulli beam (e.g.infinitesimally thin and wide).With this model, the only non-zero stress is σ 11 = Y ε 11 which would produce Mohr-Coulomb/tensile failure when σ 11 = −σ c /q. Hence the flexural strength can be used to estimate the small scale cohesion by The lab measurement of cohesion (τ S 0 = 1.1 MPa, Schulson, 2009, also see Table 1) used a sample with v b = 0.05, so σ est f ≈ 473 kPa and τ S,est 0 ≈ 454 kPa -that is, the estimated failure stress and cohesion are too small, by a factor of approximately 2.42.A similar factor was obtained by Marchenko et al. (2014), who used a full finite element 3D solver (COMSOL) to estimate the stress at the fixed end of a cantilever at the time of breaking, and found it to be approximately 2.6σ est f .Now, the results of these simulations depends on the boundary conditions used (e.g. the properties of the spring foundation used; free surface conditions when the ice was partially submerged), and in addition some predictions were not observed (e.g. they predicted the force measured in the tests should increase when the radius of the holes drilled near the beam root increased: Marchenko et al., 2017).However, it gives further indication that σ est f could definitely be a significant underestimation for the actual breaking stress.If we wanted to be consistent with the lab scale measurement of the cohesion over a range of brine volume fractions, we could propose the relationship τ S 0 ≈ 2.42τ S,est 0 ≈ 2.33σ est f .In practice though, the sensitivity studies are conducted by varying the small scale cohesion directly, and seeing the range of MIZ widths obtained.However, more observations with regard to ice breakage by waves are needed to set a definitive breaking criterion.Some laboratory experiments to this effect are planned to occur in 2018 in the wave/ice tank in Aalto, Finland, as part of the Hydralab+ programme, but field observations would also be very useful.
When we return to our plane wave in an elastic plate, the Mohr-Coulomb criterion is equivalent to the strain criterion ) is typically thought to be about 3 − 10 × 10 −5 (e.g.Langhorne et al., 1998), but this number contains a lot of assumptions, e.g. about the value of Young's modulus and the stress at the time of breaking (see the discussion below about the flexural strength).In fact, we are not aware of any strain measurements for ice which actually broke.Langhorne et al. (2001) measured strains up to about 3.6 × 10 −6 in landfast ice which was experiencing incoming waves but which did not break.When we have a spectrum of waves, the corresponding quantity to ( 22) is related to the maximum mean square strain by If all the wave energy is travelling in one direction (which direction is not relevant since we also do not attempt to consider an anisotropic wave medium), equation ( 26) is still equivalent to the Mohr-Coulomb criterion since we still have σ 2 = νσ 1 .
However, we now have a statistical (approximately normal) distribution of strains max{ε 11 }, instead of a fixed strain amplitude.
Thus ( 26) corresponds to a condition on the probability of max{ε 11 } exceeding ε c where P c is some critical probability.An alternative to (26) could be to choose P c another way (e.g.defining it as the ratio of a breaking time scale to the mean wave period), or else P(max{ε 11 } > ε c ) could be used directly in a similar formulation to Horvat and Tziperman (2015).However, for now we use (26) so that the criterion agrees with the criterion for a plane wave (e.g. a swell wave).
When the wave energy is not unidirectional, the stresses are no longer distributed on the line σ 2 = νσ 1 , so the probability condition ( 28) is no longer equivalent to the Mohr-Coulomb criterion.A simple numerical experiment generating random waves in an ice sheet and creating an artificial time series (not shown) found that P(max{ε 11 } > ε c ) was significantly lower than the probability of the stresses leaving the failure envelope (about 45% compared to about 65% in one example).However, for now we will leave this as a caveat and attempt a fuller investigation of the Mohr-Coulomb failure in a random sea at a later date.

Ice break-up
When ( 26) is satisfied, we calculate the mean zero crossing frequency from and convert this to a wavelength λ 02 using the dispersion relation for a thin elastic plate (Williams et al., 2013a, Appendix A).Then D max is reduced to λ 02 /2 (requiring that it stays above D min = 20 m, and that it is actually reduced -i.e. it can't increase, since we don't consider thermodynamic effects in this paper).

Momentum loss due to attenuation
Following Phillips (1977, Chapter 3), we first connect the mean energy per unit area (integrated over the entire water column) for a single plane wave to the mean momentum per unit area.The mean kinetic energy density is where u w and v w are the horizontal and vertical wave orbital velocities, and Z is the water depth.In a conservative system, the mean potential energy and the mean kinetic energy are equal, so the mean energy density is simply The mean momentum per unit area is: where c p = ω/k is the phase velocity.
When we consider a complete wave spectrum, then and its flux is This quantity can then be transferred to the ice, ocean and atmosphere, according to the different attenuation mechanisms, i.e.
For this study we assume that all the momentum goes to the ice -i.e.τ w,o = τ w,a = 0. Between calls, these will have changed due to dynamic (advection) and thermodynamic processes (melting, freezing).These are interpolated from the neXtSIM mesh to the WIM grid, and D max is retrieved from N floes .After the call to the WIM, N floes is passed back onto the centres of the mesh, and the stresses τ w,i are interpolated from the grid centres onto the nodes of the mesh, and are used in the solution of the momentum equation.These stresses are kept constant until the next call to the WIM -since the mesh is moving, this requires re-interpolation at each neXtSIM time step.
In an initial more naive implementation of the coupling, N floes was computed only on the WIM grid, then interpolated back onto the mesh.However, passing this field to and fro between the mesh leads to a large amount of numerical diffusion.To solve this problem, the WIM model takes in the neXtSIM mesh, and each WIM timestep the smoother integrals m 0 , m 2 and m ε are interpolated from the grid to the mesh.This allows the breaking calculation to be done on the mesh in parallel to the one on the grid -thus N floes does not need to be interpolated back to the mesh.This also reduced the diffusion in N floes significantly.
(See Figures 7-8 below.) The directional wave spectrum is remembered from the previous call, and if necessary can be updated regularly using forcing from an external model, or as in the simulations presented in this paper, using idealised (constant) wave forcing.
We can also change the dynamics of the broken ice.The default, "R0" or rheology 0, does not change the underlying EB rheology.In an alternative, "R1" or rheology 1, we increase the damage parameter d to an arbitrary high value d broken when the ice is broken by waves.This reduces the internal stress, apart from a pressure term which resists compression, causing the ice velocity to be closer to the free drift velocity.
Alternative continuum approaches to MIZ dynamics are based on the idea of a "granular temperature" (kinetic energy associated with velocity fluctuations relative to the mean flow field).Most recently, Feltham (2005) used a binary collision model to formulate an equation for the granular temperature.Previously, Shen et al. (1986Shen et al. ( , 1987) ) had used a similar but simpler approach, where the granular temperature was approximated to be in steady state.This enabled the granular temperature to be found analytically and the constitutive relation to be directly modified without solving any other equations apart from the momentum equations.Shen et al. (1987) compared the granular temperature to field data from the MIZEX campaign of 1983 (Hibler and Leppäranta, 1984), and found it to be correlated, but found that it was an order of magnitude too small.The internal ice stresses were also very low.Feltham (2005) was able to produce some qualitative features such as ice jets in a onedimensional simulation, but no further comparisons were done.This model is now being introduced into CICE-E (Community ICE code, version E; Rynders et al., 2016).
However, in the field of 3D granular flows, different types of flow regimes have also been observed.For example, the introduction of Guo and Campbell (2016) describes a transition between an inertial collision regime to an inertial non-collisional regime where the stresses follow Bagnold's law (Bagnold, 1954) as the concentration and shear rate increase, and then a further transition to what they call the elastic regime as the concentration and shear rate increase even more.This regime is characterised by the formation of force chains at high concentrations and shear rates, which deform elastically to support the applied stresses.
There have also been a number of direct (discrete) numerical simulations of collections of floes (e.g.Herman, 2013;Rabatel et al., 2015).They have also observed phenomena similar to the force chains mentioned above, where elaborate force contact networks were observed over the full domain of simulation.To summarise, the binary collisional models represent only a small fraction of the types of granular flows observed, so there is much more work required before a complete "MIZ rheology" that could be substituted for our simple modification is ready.

Note on wave and wind forcing
In our results section we will partly use incident wind wave spectra based on the Bretschneider spectrum: where H s is the significant wave height, ω p = 2π/T p , and T p is the peak period.
Since H s and T p are not totally independent, to try to make them roughly consistent we will also use a special case of (36), the Pierson-Moskowitz spectrum which was defined as an approximation for fully-developed wind seas: where a PM = 8.1 × 10 −3 , b PM = 0.74, and ω 0 = g/U 19.5 ≈ g/(1.026U 10 ).Here U 19.5 and U 10 are the wind speeds 19.5 m and 10 m above the sea (respectively) -note that these wind speeds are linked to the incident wave parameters, and we will also try to keep them consistent when we are presenting coupled WIM-neXtSIM results.The Bretschneider parameters corresponding to the Pierson-Moskowitz parameters are: Our incident wind wave spectra will then combine a Bretschneider frequency spectrum with some directional spreading: where H is the Heaviside step function.(Note that the mean wave direction is zero, ie to the right in our model domain, which can be seen in Figure 1.)We will also look at so-called swell waves, which are not locally generated, generally quite long (wave period greater than about 10 s or longer), and are monochromatic and mono-directional:

Sensitivity of MIZ width to Young's modulus and small-scale cohesion
The purpose of this section is to test sensitivity to the Young's modulus and the small-scale cohesion, not necessarily to decide on "correct" values, which are best determined from future observations.The experiments are similar to those of Williams et al. (2013b), although the effect of the Young's modulus was not tested in that paper.This is an interesting parameter since increasing it makes the ice less compliant and easier to break (ie.a given wave amplitude produces a higher stress in the ice)potentially increasing the MIZ width -but this also increases the attenuation, which could potentially reduce the MIZ width.
The effect of the small-scale cohesion will play a similar role to the breaking strain in that paper.
The Young's modulus is typically somewhere in the range of 1-10 GPa.Williams et al. (2013a) argued for values within the interval 5-7 GPa (depending on the brine volume fraction), proposing that the effective elastic modulus, which includes a response to primary, recoverable creep, should cause it to drop somewhat from the relationship of Timco and Weeks (2010).
However, Marchenko et al. (2013) derived significantly lower values of Young's modulus (about 1.5 GPa) in Svalbard fjord ice.Marchenko et al. (2017) also measured lower values in the Barents Sea, ranging between 1-4 GPa, with no obvious dependance on the brine volume.Therefore, we do some tests of the sensitivity of the MIZ width and the maximum WRS to this parameter.
Figure 4 shows the variation of the MIZ width (panel a) and the maximum WRS (panel b) with peak period for different values of the Young's modulus.Since increasing the Young's modulus increases the attenuation, the waves lose more momentum and so the maximum radiation stress increases, and this is clearly seen in Figure 4(b).However, Figure 4(a) clearly shows that the MIZ width increases with increasing Young's modulus, so its effect on the breaking criterion clearly dominates its effect on the attenuation.The magnitude of the maximum radiation stress is of the order of 0.1-1 Pa, which is comparable to the wind stress from a 10-15 m s −1 winds (see Figure 1d).However, while stresses of this size are significant, they are very much localised around the ice edge as opposed to being applied over large areas (as wind stresses are -see Figure 1(d)).Moskowitz spectra are used for the forcing.Solid curves: Bretschneider spectra are used with the significant wave height being 4m.The concentration was 0.7, the thickness was 1 m, and the small-scale cohesion used was 629 kPa.The WIM is not coupled to neXtSIM.
The dashed curves use fully-developed seas (Pierson-Moskowitz spectra), where H s increases with T p , for wave forcing.
Although waves of higher periods are attenuated less, the increasing wave height overcomes this effect and both the MIZ width and maximum radiation stress increases monotonically with peak period.
The solid curves in Figure 4 are created using an incident wave spectrum based on a Bretschneider spectrum with a constant significant wave height of H s = 4 m.Like with the dashed curves (fully developed seas), larger values of Young's modulus cause the MIZ width to increase monotonically as peak period increases (in the plotted range of periods).However, when Y = 1 GPa, as peak period is increased, the MIZ width is initially 8km, then increasing to a maximum of 12 km as the wave frequencies with the most energy are attenuated less, before dropping down to 8 km again as the waves with the most energy, while still being attenuated less strongly, now produce less strain (see equations (22)(23)(24)(25)(26)(27)).This latter result (Y = 1 GPa, constant wave height) is similar to results for constant-amplitude swell waves, plotted in Figure 5 -very low periods are attenuated too strongly to do much breaking so the MIZ width is zero; above a certain period the MIZ width increases (with period) to a maximum then drops back down to zero when the induced strain is no longer large enough to cause breakage.For this wave height of 3 m, which is relatively large, but not unrealistic for the usual range of swell periods (ca.10-20 s), the maximum radiation stress drops from about 0.1 Pa to about 0.01 Pa showing the reduced ability of swells to produce wave drift in comparison to wind seas.
Figure 6 shows the variation of the MIZ width with the peak period and the small scale cohesion.Unlike the Young's modulus, this parameter does not change the attenuation directly, and so the maximum radiation stress is essentially the same for all values of the cohesion (notwithstanding small differences, mainly due to the different MIZ widths, since the attenuation is higher in the MIZ in our model).The concentration was 0.7, the thickness was 1 m, and the small-scale cohesion used was 629 kPa.The WIM is not coupled to neXtSIM.Pierson-Moskowitz spectra are used for the forcing.Solid curves: Bretschneider spectrum are used with the significant wave height being 4m.(b) Swell waves of height 3 m.For both plots, the concentration is 0.7, the thickness is 1 m, and the Young's modulus used was 5.49 GPa.
The WIM is not coupled to neXtSIM.
The three values chosen are 270 kPa (approximately the flexural strength when v b = 0.1, 274 kPa), 629 kPa (approximately 2.33 × 274 = 638 kPa), and 1.08 MPa (approximately 1.1 MPa, the laboratory value of the cohesion).The results for the MIZ width are significantly different but all are in the correct order of magnitude (a few tens of kilometers).Therefore we will use τ 0 = 629 kPa throughout the rest of the paper.We will also use a Young's modulus of Y 0 = Y * = 5.49 GPa (i.e. the same value in neXtSIM and the WIM).

Coupled waves-in-ice results
Figure 7 shows plots of different fields after a 2-day simulation with neXtSIM coupled to the WIM.There is no wind, only waves arriving from the left (the initial wave state is shown in Figure 7(a)), breaking the ice and pushing it to the right by about 24 km by the end of the 48-h simulation.The initial ice state is the same as in Figure 1, but with the addition of unbroken ice (D max = 300 m everywhere where c > 0), as shown in Figure 7(b).This could correspond to summer ice in the Fram Strait where there can be large floes with large gaps between them (perhaps due to smaller floes melting faster), producing a low concentration.
The resulting MIZ width is about 50 km, which is not unrealistic.Following (39), there is a cos-squared type of directional spreading applied (and 16 directions used) and the upper and lower grid cells, which contain land, act to completely absorb the waves.Therefore, in Figure 7(c), the waves are slightly lower (by about 1 m) near the coast than they are at the centre.In Figure 7(f), the x-component of the WRS is plotted -note that while it reaches 1 Pa in the vicinity of the ice edge, it decays exponentially further into the ice.This is reflected in the concentration field (Figure 7(e)), which shows that the ice is much more compact at the ice edge.Note that the WRS is not varying significantly in the y direction, showing that the boundary conditions used for the waves at the coast are not having too much influence.Also note that the pack and the MIZ, as shown in the D max field (Figure 7(d)), are separated by quite a sharp boundary.This has been preserved by doing breaking on the mesh in parallel to the breaking on the grid, as opposed to simply interpolating D max back to the mesh after doing breaking on the grid.Figure 8 shows the same plot as Figure 7(d), but with this latter, more naive, method of coupling.The sharp MIZ-pack boundary has now become extremely diffuse compared to the former scheme.
Figure 9 tests the sensitivity of the ice edge motion to the rheological parameters C and τ L 0 when the ice is subjected to steady waves of varying heights (and periods).In Figure 9(a), the damage is set to 0.9999 everywhere the ice is broken by the waves, while in Figure 9(b) the damage and cohesion are unchanged by ice breakage due to waves.Consequently in Figure 9(a) for higher concentrations the internal stress is mainly coming from the ice pressure P , while in Figure 9(b) σ also plays a role since it is not damaged.
There is a strong response to the compactness factor, C, which is used in the neXtSIM model to determine how high the concentration needs to be to increase the effective elastic stiffness and the resistance to ridging to their maximum values.
In Figure 9(a), for this initial value of concentration (70%), lowering C by 10 roughly reduces the ice movement by half.
Comparing Figure 9(b) to Figure 9(a), if C = 40, σ makes a difference of between 8-15 km; if it drops to 30, the ice edge movement is approximately reduced by half; if it drops even further to 20, then the ice edge no longer moves at all.However, the large-scale cohesion makes little difference in these simulations where the ice is not failing.Part of the reason for this is that the wave radiation stress is a compressive stress, so the stresses need to be larger to move outside the Mohr-Coulomb envelope than if they were tensile or shear stresses (see Figure 2: the tensile and shear stresses are near the points of the triangles, while compressive stresses are near their bases).Some of the runs from Figure 9 (those with C = 40 and τ L 0 = 4 kPa) were repeated with swell waves (of a single frequency and direction), with amplitude of 3 m and periods ranging from 10-14 s (recalling that the maximum WRS dropped with wave period -Figure 5(b)).These were not able to produce any movement of the ice edge though.Therefore, the main influence of swell will be due to their changing of the dynamical and thermodynamical properties of the ice through the ice break-up.As can be seen from Figures 5-6, they are attenuated less and so they can produce break-up further into the ice than wind waves.
Figure 10 shows the combined effects of wind and waves on the concentration (c) and the effective thickness (ch).For reference Figures 10(a,b) have only waves (5-m waves following a Pierson-Moskowitz spectrum) and no wind (Figure 10(a) is the same as Figure 7(e)), while Figures 10(c,d) have no waves, but only a 15 m s −1 wind from the left (as in Figure 1).This wind speed is consistent with the wind wave spectrum in Figures 10(a,b).Figures 10(e-h) have both 5-m waves and 15 m s −1 wind.
All figures with wind Figure 10(c-h) exhibit similar ice edge locations, and all show thickening at the far right "coastline", concentrated in thin "ridges".The area over which the ridging is concentrated also seems similar for all the runs.However, while the pattern of thickening between the three runs seems quite different, perturbations to certain parameters in the run with the R1 modification to the EB rheology (Figures 10(g,h)), such as d break (0.99 and 0.999 were tried), or the minimum concentration of ice required to cause attenuation (0 or 5% were tried), produce similar degrees of differences.Therefore we conclude the actual ridging patterns are not significant in themselves.The main differences therefore between the R1 run and the other two are therefore in the concentrations at the ice edge (the actual thickness, h, which is not plotted, is constant near the edge).In this run, when the damage is increased if ice breakage occurs, the ice is noticeably more concentrated in a region approximately corresponding to the MIZ.Additionally, the ice edge is more diffuse, possibly due to some feedback effect where if the ice begins to become less concentrated at the ice edge, the attenuation reduces and therefore so does the wave radiation stress, and then moves more slowly compared to the more concentrated ice which will experience a higher radiation stress -an effect enhanced by the high degree of damage which keeps the more compressed ice quite mobile (as opposed to the run where the rheology is not modified).
Figure 11(a) quantifies the results of Figure 10 with respect to the ice edge location, as well as varying the wind speed.As can be seen from the figure, the waves only increase the movement by 4 km (no damage in the MIZ due to breakage) or 8 km (damage is d break = 0.9999 in the MIZ when the ice is broken).That is, the effect of the WRS on the ice edge position is almost completely dominated by the wind stress.When the initial concentration was increased to 95%, the difference was even less (0-4 km), as then the stress and ice pressure P increased due to their e −C(1−c) factors becoming closer to 1.
To repeat what we have seen in Figure 10, when the ice was subjected to on-ice winds in addition to waves, the main effect of linking the damage to the break-up due to waves was that the MIZ region became more highly compressed than the ice immediately further in.In Figure 11(b), we see the effects of off-ice winds on ice preconditioned by swell waves.For the wind speed used in the figure shown (2 m s −1 ), the wind stress is not able to move the pack ice at all, but the MIZ, which is about 60 km wide and has damage d break = 0.9999, has started to detach from the pack.The ice edge has moved about 15 km to the left in the centre of the domain, with less movement at the coasts since there is still some friction there (due to the condition of no slip applied at the top and bottom boundaries).

Conclusions and discussions
In this paper, we have investigated the impact of the WRS on sea-ice state and drift in an idealised domain.While this stress can be quite large (∼.1-1Pa), depending on the wave conditions, it is extremely localised -decaying exponentially away from the ice edge.Probably as a consequence of this localisation, overall we found its effects on ice edge location were quite   e-h) show results when steady waves (with Hs = 5 m, Tp = 11.2 s, from the left) are applied in addition to the wind forcing.Initial ice conditions are the same as in Figure 7.In (e,f) the ice rheology is not affected by the ice breakage, but in (g,h) damage is set to dbreak =0.9999.The large-scale cohesion is τ L 0 =4 kPa, C = 40, and the small-scale cohesion is τ S 0 =629 kPa.
modest, with the most noticeable effects being seen when a wind wave spectrum was applied steadily to the ice in the absence of wind.Then, depending on the initial concentration, the rheological parameters used and the response to the ice breakage by waves, the radiation stress could produce a movement of the ice edge of between 0-36 km over two days.However, this experiment is more hypothetical since wind waves are by definition associated with wind.Indeed, in the presence of wind, the wind stress dominated the WRS with almost no difference in ice edge position between experiments with and without waves.There were differences in ridging patterns in the presence of waves but these were probably not significant.However, forcing only ("no waves"), wind and waves without changing the EB rheology in the MIZ ("no damage"), and wind and waves where the damage is set to dbreak = 0.9999 in the MIZ ("0.9999").(b) The effect of swell preconditioning on the response to off-ice winds.Initial conditions: swell of 12 s period and height 3 m is sent into the ice for 24 h, where the ice has constant concentration of 95% and thickness 1 m, breaking the ice for about 50 km.The damage is set to dbreak = 0.9999 where the ice is broken.Spatially and temporally constant off-ice wind forcing is then applied for a further 48 h, at a speed of 2 m s −1 .The large-scale cohesion is 13 kPa, and C = 40.
when we modified the damage parameter after ice breakage, additional compression was observed in the MIZ after the ice was broken.Consequently, it seems that the WRS has a very limited effect in general, although it could be a very efficient process to precondition the ice cover and its mechanical properties via the formation of a MIZ area filled with highly damaged ice.
Having said this however, there are many uncertainties regarding the WRS, and we have certainly not included all of its potential effects, especially since the wave and ice models are not coupled to the ocean yet.For example, the attenuation models are still uncertain (they determine the WRS), and how the partitioning of the WRS between the ice and the ocean should be done is also unknown.On the face of it, if less of the WRS is applied to the ice, it should have even less effect than we find in our current paper.However, perhaps it could then produce similar effects to those discussed and reported by Suzuki We also highlighted the problem of numerical diffusion of N floes due to it being modified by both neXtSIM and the WIM, and therefore having to be communicated in both directions.We presented a solution to this problem, where N floes was calculated on the neXtSIM mesh each WIM time step, after interpolating smoother wave fields.While not unfeasible, this is somewhat costly and we will continue to look for alternative solutions.
As touched on in the discussion of the WRS above, we also introduced a simple MIZ rheology by increasing the damage where ice was broken, effectively putting the MIZ into free drift, with the addition of the ice pressure which resists compression.
Under compressive wind forcing this led to increased compression in the MIZ relative to the pack ice in its vicinity.This modification also influenced the ice flow when off-ice winds were applied to ice that had previously been broken by swell waves.At lower wind speeds, the MIZ was able to be move relatively freely with the wind, while the pack was still stationary.
These effects would undoubtedly be reduced in magnitude were a rheology that represented true granular flow to be used, but could still occur.However, it is difficult to know for certain without the existence of such a rheology.Direct numerical simulations such as those done by Herman (2016) could possibly reproduce some of the effects observed here.Similarly, the granular temperature model of Feltham (2005) could be tried, although this would be limited to flow regimes where large force networks are not expected to be present.
So far we have also restricted ourselves to a simple idealised domain, and with very idealised forcings.Work to set up the current model in a pan-Arctic domain is ongoing, and perhaps studies with forcings with more realistic temporal and spatial variability could find the WRS will have more impact.In addition, the study of Horvat et al. (2016) suggests that including the thermodynamic effects of ice breakage by waves could be important.We are also currently implementing the more conservative lateral melting model of Steele (1992) in our model to include this effect to some extent.With simulations using a WIM coupled to a stand-alone version of CICE-E, which contains the model of Steele (1992), Bennetts et al. (2017) found that the concentration in the vicinity of the Antarctic ice edge could drop by a modest amount (of the order of 10%) in the summer.However, this could also change with coupling to an ocean model, as well as if a different parameterisation that reflects the increased lateral melting of larger floes were used.

Code availability
This code is not publicly available.

Data availability
This data is not publicly available.
Author contributions.The paper writing and implementation of the coupling between the WIM and neXtSIM was lead by TW, with formative discussions from PR and SB guiding the progression of the writing.PR also helped with the writing itself, and in addition SB helped implement the coupling.

Competing interests. N/A
Disclaimer.N/A

Figure 1 .
Figure 1.Results after forcing from uniform, steady wind (with speed 14.9 m s −1 , from the left) has been applied for 48 h.Initially, constant ice conditions were applied (c = 0.7, h = 1 m, d = 0) were applied to the right of the ice edge, which corresponded to approximately x = 124 km.The upper, lower and right hand boundaries are closed.Fields plotted are (a) concentration, (b) effective thickness, (c) the damage (with blue being more damaged and red less), and (d) the x-component of the wind stress τ a.There are no wave interactions considered, C=40 and τ L 0 = 4 kPa.
) here, α is the scattering per floe, while β is the imaginary part of the wave number satisfying the dispersion relation ofRobinson and Palmer (1990), calculated using the method ofWilliams et al. (2013a, Appendix A)  with drag coefficient Γ = 13 Pa s m −1 .

3. 4
Ice breakage due to waves 3.4.1 Plane strain and Mohr-Coulomb failure It is instructive to put the situation of ice breakage due to a plane wave in the context of the discussion in §2.3.We also use a thin elastic plate model, so the constitutive relation is similar to equations (1-2): σ = C(Y, ν)ε, where Y is the Young's

Figure 2 Figure 2 .
Figure 2(a) plots the failure envelopes for two values of the cohesion.The figure also shows where the line corresponding to the stress state for plane waves, σ 2 = νσ 1 , meets these Mohr-Coulomb envelopes (i.e., when σ 1 = σ (tens) 1).

)
Williams et al. (2013a) imposed a strain criterion for breaking, supposing that ice would break if ε ≥ ε est c = σ est f /Y , where σ est f is the flexural strength estimated from measurements.Timco and Weeks (2010) compiled many measurements for the flexural strength, fitting the formula10 −6 σ est f = 1.76e −5.88 √ vb ,(23)where v b is the brine volume fraction.(It should be noted however, thatKarulina et al., 2013, found a different relationship instead of using ε est c .Due to cancellation of unrelated but similar factors this is approximately the same as the breaking strain of Williams et al. (2013a) (σ est f /Y ).This (ε c ) is plotted in Fig 2(b) as a function of Y .The breaking strain for sea ice (from beam tests

Fig 2
Fig 2(b) shows the breaking strains are about the right order (5 × 10 −5 is plotted as a dotted line for reference), although higher values of the cohesion combined with lower values of Young's modulus can take them up to 10 −3 .

Figure 3 .
Figure3.Schematic showing the information passed between neXtSIM and the WIM.Note that Nfloes is modified by both the WIM and neXtSIM (which use different grids), so must be treated carefully to avoid numerical diffusion.Also input to the WIM are incident wave fields, and it also outputs diagnostic fields of the waves in the ice.Optional: the WIM may also update the damage d.

Figure 4 .
Figure 4. Variation of MIZ width (a) and maximum WRS (b) with peak wave period and Young's modulus.Dashed curves: Pierson-

Figure 5 .
Figure 5. Variation of MIZ width (a) and maximum radiation stress (b) with peak wave period and Young's modulus for swells of height 3 m.

Figure 6 .
Figure 6.Variation of MIZ width with peak wave period and small-scale cohesion for (a) wind seas and (b) swells.(a) Dashed curves:

Figure 7 .
Figure 7. Waves breaking ice in an idealized experiment (the right hand, upper, lower lines of grid cells correspond to land).The wave model, based on (Williams et al., 2013a), is coupled to the neXtSIM sea-ice model.The figure shows results after 48 h of steady pushing by a Pierson-Moskowitz wind wave spectrum with significant wave height Hs = 5 m (so the peak period Tp = 11.2 s), that is arriving from the left.It initially occupies the strip shown in (a) then travels to the right, with some directional spreading; the final wave height is shown in (c).(b,d): initial, final maximum floe size (respectively); (e,f): final sea-ice concentration and x-component of the wave radiation stress (respectively).The ice has initial conditions (constant where there is ice-see (b) for the initial ice mask): c = 0.7, h = 1 m, Dmax = 300 m, and d = 0. Also C=40, τ L 0 = 4 kPa, τ S 0 = 629 kPa, and d is increased to dbreak = 0.9999 if the ice is broken.

Figure 8 .Figure 9 .
Figure8.Same figure as 7b), but where Nfloes is simply interpolated from the neXtSIM mesh onto the WIM grid and then back again after each coupling time-step.Note the boundary between the pack ice and the MIZ has diffused over a large number of grid cells, whereas it has remained much sharper when Nfloes is calculated directly on the neXtSIM mesh.

Figure 10 .
Figure 10.(a,b) Concentration (c) and effective thickness (ch) after the same experiment as Figure 7. (c,d) Same as Figure 1(c,d): wind forcing only.Figures (e-h) show results when steady waves (with Hs = 5 m, Tp = 11.2 s, from the left) are applied in addition to the wind

Figure 11 .
Figure 11.(a) shows the maximum movement of the ice edge as a function of wind speed.The different curves show the response to wind

and
Fox-Kemper (2016)  andSuzuki et al. (2016)  in relation to overturning circulation produced by the Stokes shear force and thereby change the currents and heat fluxes acting on the ice.

Table 1 .
Table 1 shows the Mohr-Coulomb parameters, and the estimated defect sizes, which have been fitted to various time series of stress measurements.Cohesion values, internal friction coefficient from measured Mohr-Coulomb failure envelopes.Also given are approximate defect