Wave–sea-ice interactions in a brittle rheological framework

Abstract

As sea ice extent decreases in the Arctic, surface ocean waves have more time and space to develop and grow, exposing the marginal ice zone (MIZ) to more frequent and more energetic wave events. Waves can fragment the ice cover over tens of kilometres, and the prospect of increasing wave activity has sparked recent interest in the interactions between wave-induced sea ice fragmentation and lateral melting. The impact of this fragmentation on sea ice dynamics, however, remains mostly unknown, although it is thought that fragmented sea ice experiences less resistance to deformation than pack ice. Here, we introduce a new coupled framework involving the spectral wave model WAVEWATCH III and the sea ice model neXtSIM, which includes a Maxwell elasto-brittle rheology. This rheological framework enables the model to efficiently track and keep a “memory” of the level of sea ice damage. We propose that the level of sea ice damage increases when wave-induced fragmentation occurs. We used this coupled modelling system to investigate the potential impact of such a local mechanism on sea ice kinematics. Focusing on the Barents Sea, we found that the internal stress decrease of sea ice resulting from its fragmentation by waves resulted in a more dynamical MIZ, particularly in areas where sea ice is compact. Sea ice drift is enhanced for both on-ice and off-ice wind conditions. Our results stress the importance of considering wave–sea-ice interactions for forecast applications. They also suggest that waves likely modulate the area of sea ice that is advected away from the pack by the ocean, potentially contributing to the observed past, current and future sea ice cover decline in the Arctic.

1 Introduction

The interactions between ocean surface waves and sea ice have been receiving a significant amount of attention in recent years, particularly motivated by the decreasing Arctic sea ice extent resulting in larger areas of open water exposed to the wind and thus available for wave generation. As a consequence, wave events in the Arctic are expected to be more frequent and more intense with waves penetrating far into the ice cover and breaking large sea ice plates into floes of less than a few hundred metres see, e.g.. The attenuation of waves by sea ice, however, limits this fragmentation to the interface between the open ocean and the pack ice in the so-called marginal ice zone (MIZ). The MIZ is a highly dynamic area characterized by strong interactions between the ocean, sea ice and atmosphere. State-of-the-art sea ice models used in climate prediction systems have been shown to fail at representing the complexity of these interactions and have their biggest errors in the MIZ . On shorter timescales, large uncertainties remain in the forecasts of the position of the sea ice edge ; however, this information is essential for the safety of the increasing number of human activities in polar regions . These inaccuracies can certainly be attributed (at least in part) to the lack of representation of some of the processes occurring in the MIZ and the impact of the waves on sea ice dynamics is one of them.

Waves can impact sea ice dynamics in the MIZ through a variety of processes. For instance, wave attenuation transfers momentum from waves to sea ice through wave radiation stress WRS,, which acts as a force that pushes the sea ice in the direction of the incident waves. Being mostly directed on-ice, its main effect is to maintain a compact sea ice pack near the ice edge, but its importance is still being discussed. Estimating wave attenuation from synthetic aperture radar (SAR) images, found it to be as important as wind stress over the first 50 km of the MIZ in the Southern Ocean, whereas did not observe any wave-induced sea ice drift in pancake ice in the Southern Ocean from in situ measurements despite a strong wave-in-ice activity. Fragmentation is also likely to change the mechanical properties of the ice, but the evolution of dynamical and mechanical properties of sea ice cover with floe size remains poorly understood.

Intuitively, we expect that broken ice is more mobile than continuous ice e.g., having lower internal stress. This seems to be consistent with the deformation observations of collected by ship radar during the N-ICE-2015 expedition. In their observations, deformation in fragmented sea ice was an order of magnitude higher than in the pack. However, in the absence of routinely available datasets providing synoptic information on sea ice drift, wave height and floe size in the MIZ, observations do not allow us to arrive at any explicit relationship between the level of fragmentation and the ability of sea ice cover to be deformed.

A few attempts have been made to relate floe size to sea ice dynamics in theoretical models. used a collisional stress term accounting for floe size to represent sea ice behaviour in the MIZ. The fluctuations of the velocity field they obtain, however, did not reproduce observations from the MIZEX campaigns; they were too small by an order of magnitude. used a similar collisional stress but allowed the velocity fluctuations to be time dependent. He showed that it enables the generation of ice jets in the MIZ but on smaller scales than the reported observations of ice jets by . Dynamics in state-of-the-art sea ice models, however, do not account for the floe size. Instead, the region of the MIZ where sea ice behaves almost in free drift is mostly a function of sea ice concentration. The potential impact of fragmentation in compact ice is therefore neglected, whereas regions of low sea ice concentration and high wave activity do not necessarily coincide . analysed a cyclone and showed how sea ice in the Antarctic MIZ is visibly deformed on sea ice concentration maps despite being highly compact. They stated that this behaviour could not be reproduced using a concentration-only criterion to dynamically distinguish pack ice from the MIZ, stressing the need to account for other properties like the floe size. However, linking floe size, or fragmentation, to sea ice dynamics remains a challenging task, first because of the limited data available but also due to the poor understanding of the way waves propagate in the MIZ, although modelling progress has been made in this particular field.

Modelling efforts relating to waves-in-ice have progressed a lot in recent years, although the heterogeneous nature of sea ice and the wide variety of wave attenuation processes in the MIZ still make wave prediction in ice highly challenging . The importance of each wave attenuation process varies with wave and sea ice properties. Scattering, for instance, is efficient to attenuate short waves in fragmented sea ice covers made of consolidated floes , while dissipative mechanisms, like under-ice friction, are expected to dominate in the case of forming ice and long swells propagating in the pack. A sequence of reviews by and gave a more detailed history of this area of work. and have stressed the importance of the floe size on wave attenuation. These reports have motivated the implementation of interactions between waves and floe size, through fragmentation, in numerical models. The first studies used one-dimensional models to look at the feedback between ice break-up and wave attenuation . These models assume that break-up occurs if wave-induced flexural stress overcomes sea ice strength and that the resulting floe size distribution (FSD) follows a truncated power law with its upper limit (often called maximum floe size, $D_{\text{max}}$) depending on the wave field . This assumption about the shape of the FSD is based on observations made by of power-law FSD in the MIZ that they explain by successive fragmentation of floes by waves. More recently, included a parameterization in the spectral wave model WAVEWATCH III WW3: that also assumes a power-law FSD. Their parameterization enables interactions between sea ice floe size and wave attenuation processes, like scattering and inelastic dissipation, and was shown to explain well the wave height evolution during the ice break-up event reported by . evaluated this model by comparing their results to remote sensing and field measurements during a storm event in the Beaufort Sea, showing good agreement for the measured and modelled wave-in-ice attenuation and broken sea ice extent. Sea ice representation in these wave models remains, however, too simplistic to more deeply investigate the impact of waves on sea ice. It has led to the development of coupled wave–ice models, first using one-dimensional wave-in-ice models and more recently using more complex spectral wave models like WW3 . These developments were made possible by the implementation of wave–ice interactions in state-of-the-art sea ice models, particularly representations of the FSD . These FSDs have been mostly used to investigate the effects of lateral melting on sea ice properties over timescales of a few weeks to a few years in a context where lateral melt is expected to be enhanced by the wave-induced sea ice fragmentation .

The dynamical aspects of wave–ice interactions have received less attention. found that the WRS could regionally impact sea ice melt and sea surface properties in the Arctic MIZ at the end of summer. Concerning the impact of wave-induced fragmentation, suggested combining the classical elasto-visco-plastic rheology used in most sea ice models with a granular rheology in the MIZ to better represent floe–floe interactions. This granular rheology depends on the floe size. Numerical simulations with this approach show an overall increase of the sea ice drift speed in the Arctic year-round compared to a reference simulation using a standard version of the sea ice model CICE . suggested another approach to relate sea ice dynamics to wave-induced sea ice fragmentation using the sea ice model neXtSIM neXt generation Sea Ice Model, in a stand-alone set-up. The elasto-brittle (EB) rheology used by neXtSIM stores a variable for sea ice called damage, which tracks the level of mechanical damage of the sea ice over each grid cell . The higher the damage, the lower the sea ice internal stress. Resistance to deformation is a function of both sea ice concentration and floe size. The originality of the study by was to link the damage variable with wave-induced fragmentation, making the extension of the ice region behaving in free drift dependent on the wave field. This was done by linking the damage variable with wave-induced fragmentation. Using idealized simulations of waves compressing ice, they showed that the movement of the ice edge was not very sensitive to either wave fragmentation or the WRS. The investigations of were, however, limited to very idealized cases, and the EB rheology in neXtSIM has now been upgraded to the Maxwell elasto-brittle (MEB) rheology see, greatly improving sea ice deformation in pack ice . This upgrade could also affect the MIZ, as it led to the removal of an ice pressure term that was added to prevent damaged ice from piling up in EB but caused the modelled deformations to deteriorate too much with MEB rheology.

In this paper, we present the results obtained with a new coupled wave–sea-ice modelling system (WW3-neXtSIM). This modelling system benefits from recent wave–ice developments in WW3 and extends the work done by in neXtSIM. We again use the damage variable to link the sea ice dynamics and the fragmentation due to waves, allowing us to represent the link between the wave-induced fragmentation of sea ice and its mobility in MIZ areas. Our model also benefits from the advancements of FSD implementations in sea ice models done by and , and of the coupling of WW3 with the sea ice model LIM3 described by . We also propose a way to incorporate some floe size “memory” of previous fragmentation events due to waves by introducing two time-evolving FSDs in each grid cell. These developments provide a coupled wave–sea-ice framework able to provide year-round regional or pan-Arctic simulations. After describing the details of our implementation, we evaluate our new coupled framework to check if the produced wave attenuation, broken sea ice extent and refreezing timescales are reasonable. We then investigate the effects of wave-induced sea ice fragmentation using a regional case study. Finally, we discuss the different assumptions made in our study and suggest perspectives for future studies.

2 Implementation of the coupling between the wave and sea ice models

In this study, we make use of the spectral wave model WAVEWATCH III^® , building on the previous developments reported by , who included an FSD in WW3 as well as some representations of the different processes by which sea ice can affect the propagation and modulation of waves in the MIZ. These attenuation processes are scattering (which redistributes the wave energy without dissipation), friction under sea ice (with a viscous and a turbulent part depending on the wave Reynolds number) and inelastic flexure. All of these processes depend on sea ice thickness, concentration and floe size. Wave attenuation increases with thickness and concentration, and tends to decrease when the floe size is lower than the wavelength, as floes are not flexed anymore. This parameterization was chosen because it was shown to reproduce well wave attenuation in two different events in the Arctic: waves breaking a continuous sea ice cover near Svalbard as reported by and waves propagating in forming ice in the Beaufort Sea . As in the study by , we assume that deviations from the ice-free wave dispersion relationship induced by the presence of ice are small and can be neglected. This is likely to be the case once sea ice has been broken .

Figure 1

Summary of the exchanged variables in the neXtSIM–WW3 coupling framework.

[Figure omitted. See PDF]

The sea ice model we use for this coupling is neXtSIM , in which an FSD is first implemented as described in Sect. . The two models are coupled through the coupler OASIS-MCT . Figure  shows the variables that are exchanged. Briefly, WW3 determines if the waves will break the ice and calculates a representative wavelength ${\mathit{\lambda }}_{\text{break}}$ in the manner of . This is then used by neXtSIM to modify the FSD as described below in Sect. . WW3 also computes the WRS, which is used in the momentum equation of neXtSIM. neXtSIM gives WW3 the sea ice concentration and thickness, and the mean and maximum floe size, which are used by WW3 to determine the amount of attenuation. The mean floe size is used to determine the amount of scattering, while the maximum floe size is used to determine the amount of dissipation due to inelastic attenuation. The evolution of the FSD and the computation of the exchanged variables are described in more detail in the following paragraphs. Parameters introduced in this section are summarized in Table .

2.1 Wave radiation stress on sea ice

As in , the WRS is computed in WW3 and sent to the ice model. This computation provides an estimate of the WRS, which is likely to be an upper bound of its real value, as it assumes that all the momentum lost by attenuated waves is transferred to sea ice, therefore ignoring a potential partitioning of this momentum transfer between the ocean, sea ice and atmosphere. In neXtSIM, the WRS is added to the sea ice momentum equation in the same way as in . As discussed by , the estimation of the WRS and its distribution in the MIZ strongly depend on the parameterization chosen for the wave attenuation.

2.2 Floe size distribution modelling

As mentioned in the introduction, this study makes use of two different FSDs to represent the evolution of the floe size. The first FSD represents the evolution of the floe size considering floes as pieces of ice separated from each other by leads or cracks. This is the FSD that would be seen in satellite images or aerial photography. In this FSD, floe size growth is governed by mechanisms like surface refreezing and floe welding as in, e.g.. In freezing conditions, sea ice forming at the surface and floes welding together can therefore turn fragmented small floes into a continuous ice cover (i.e. not separated by leads) in a few hours to a few days. This FSD is particularly relevant for thermodynamical processes like lateral melting, which are likely to be unaffected by the mechanical properties of the ice cover. In neXtSIM, this FSD is represented by the variable we call “fast-growth” FSD, $g_{\text{fast}}(D_{\text{fast}},\mathit{x},t)$, where $\mathit{x}$ is the position vector, $t$ is time and $D_{\text{fast}}$ is defined as the mean calliper diameter of the floes separated from each other by leads or cracks, as introduced by .

The second FSD considers the floe size as the length scale associated with mechanically homogeneous pieces of ice, whether they are welded with other floes by thinner ice or not. In this second FSD, the floe size grows more slowly than for the first FSD, as it takes time for the ice joining the consolidated floes to thicken. The timescale associated with the consolidation is certainly more similar to the mechanical “healing” of sea ice in neXtSIM described in with values of $\simeq$ 10–30 d. In neXtSIM, this FSD is represented by the variable we call “slow-growth” FSD, $g_{\text{slow}}(D_{\text{slow}},\mathit{x},t)$ , where $D_{\text{slow}}$ is defined as the mean calliper diameter of the floes considered as mechanically homogeneous pieces of ice.

Figure 2

Aggregate of ice floes cemented together by thin ice in the marginal ice zone (Weddell Sea, Antarctica). The size of the largest floes is on the order of 10 m. Picture taken on board the RRS James Clark Ross in March 2014. Credit: Heather Regan.

[Figure omitted. See PDF]

The picture in Fig.  illustrates the different definitions of floe size in our study. On one hand, sea ice concentration is about 100 % with a continuous sea ice cover, represented by the fast-growth FSD with unbroken floes. Processes like lateral melting are unlikely to occur in these conditions. On the other hand, consolidated floes in Fig.  are easy to distinguish from the thin ice joining them. The slow-growth FSD represents the distribution of sizes of these consolidated floes. In this case, the slow-growth FSD is dominated by small floes (on the order of 10 m). This information can be useful for the study of mechanical processes. For example, it can represent the inhomogeneous nature of the ice cover, which is particularly relevant for wave attenuation processes like scattering and flexure-induced dissipation, for which the mechanical properties and ice thickness continuity of the ice cover are the quantities of interest.

The introduction of this second slow-growth FSD is motivated by the fact that sea ice cover, as a dynamical system, exhibits memory properties that can be illustrated by, e.g. scaling laws of deformation in both time and space . An ice cover such as the one illustrated in Fig.  is likely to have a very different mechanical strength/response under external stresses compared to a consolidated continuous ice cover due to the high likelihood of break-up of the thin ice between the consolidated floes. In the case of a new wave event, the fragmentation of the ice cover is likely to occur at its weakest points, i.e. at the thin ice joints. There have been very few reports of sea ice break-up events in the literature , but our intuition is supported by the recent Antarctic observations made by , who report waves breaking ice preferentially at refrozen cracks and pressure ridges. Introducing the slow-growth FSD in our model is therefore a way to keep a memory of the floe size resulting from the last wave-induced fragmentation event in the model.

These two FSDs are implemented in neXtSIM as areal FSDs as done by, e.g., which are defined as 1$$\underset{D}{\overset{D+\text{d}D}{\int }}g(D,\mathit{x},t)\text{d}D={\displaystyle \frac{\mathrm{1}}{A}}{a}_D(D,D+\text{d}D)$$ and 2$$\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}g(D,\mathit{x},t)\text{d}D=\mathrm{1}.$$

In these definitions, $g\left(D\right)$ represents an FSD where $D$ is its associated floe size, $A$ is the total area considered around the position $\mathit{x}$ at a time $t$ and $a_D$ are the areas within $A$ covered by sea ice with floes having diameters between $D$ and $D+\text{d}D$. The value $D=\mathrm{0}$ corresponds to open water and Eq. () is equivalent to ${\int }_{\mathrm{0}+}^{\mathrm{\infty }}g(D,\mathit{x},t)\text{d}D=c$, where $c$ is the sea ice concentration. In practice, we have a number $N$ of FSD bins of constant width $\mathrm{\Delta }D$ with edges between a minimum and a maximum floe size, respectively, $D_{\mathrm{0}}$ and $D_N$. Thus, floes with sizes in $\left[D_{\mathrm{0}}{D}_{\mathrm{1}}\right]$ cannot be broken into smaller pieces, and we refer to floes with sizes in $\left[D_{N-\mathrm{1}}{D}_N\right]$ as unbroken floes. Using fixed-width bins may bias our ability to represent or examine scale-invariant behaviour , but it has the advantage of being simple and the study of the FSD evolution and its impact on sea ice is outside the scope of this study.

Both FSDs evolve as in or : 3$${\displaystyle \frac{\partial g\left(D\right)}{\partial t}}=-\mathrm{\nabla }\cdot \left(\mathit{u}g\right(D\left)\right)+{\mathrm{\Phi }}_{\text{th}}+{\mathrm{\Phi }}_{\text{m}},$$ in which $\mathit{u}$ corresponds to the sea ice velocity vector, ${\mathrm{\Phi }}_{\text{th}}$ is a redistribution function of floe size due to thermodynamic processes (i.e. lateral growth/melt) and ${\mathrm{\Phi }}_{\text{m}}$ is a mechanical redistribution function associated with processes like fragmentation, lead opening, ridging and rafting. In our implementation, we assume that the only mechanical process modifying the shape of the FSD is the wave-induced sea ice fragmentation, and ${\mathrm{\Phi }}_{\text{m}}$ is therefore the redistribution term associated with this process. The advection terms for both FSDs are identical and similar to what is done for other conservative sea ice properties in neXtSIM. The terms ${\mathrm{\Phi }}_{\text{th}}$ and ${\mathrm{\Phi }}_{\text{m}}$ differ between $g_{\text{fast}}$ and $g_{\text{slow}}$ and are described below.

2.2.1 Lateral sea ice melt/growth

In this section, we describe the implementation of the terms ${\mathrm{\Phi }}_{\text{th,slow}}$ and ${\mathrm{\Phi }}_{\text{th,fast}}$ that represent the thermodynamical redistribution of floes associated with lateral melt/growth in each FSD. The evolution of the FSD due to ice growth and melt processes is first performed in the fast-growth FSD, and is quite similar to the implementation described by :

4$$\begin{array}{rl}& {\mathrm{\Phi }}_{\text{th,fast}}=\\ & -\mathrm{2}{G}_{\text{r}}\left(-{\displaystyle \frac{\partial g_{\text{fast}}}{\partial D}}+{\displaystyle \frac{\mathrm{2}}{D}}{g}_{\text{fast}}\right)+\mathit{\delta }(D-D_N){\dot{c}}_{\text{new}}+{\mathit{\beta }}_{\text{weld}},\end{array}$$ where $G_{\text{r}}$ is the lateral melt rate of floes, ${\dot{c}}_{\text{new}}$ is the rate of formation of new ice and ${\mathit{\beta }}_{\text{weld}}$ is the FSD redistribution term associated with welding of floes using the Smoluchowski equation as implemented by .

Lateral melting is implemented following and . Here, we neglect lateral melt for the largest floe size category as floes with size $O\left(\mathrm{100}\right)$ m and more are not resolved in this study and are expected to have very little contribution to lateral melt. We also do not make any distinction between what they call the “lead region” and the “open water fraction” of each grid cell, which means the factor called ${\mathit{\varphi }}_{\text{lead}}$ in is taken to be 1. Note that lateral melt is included in the model in this study but is not discussed here, as we focus on the impact of waves on sea ice dynamics during a time period dominated by freezing.

In contrast to , sea ice is assumed to be unbroken when initialized in our model, and there is therefore no need for an explicit thermodynamical lateral growth due to the agglomeration of frazil ice at the edge of existing floes. If, after a wave-induced fragmentation event the sea ice concentration reaches 1 in freezing conditions, it is assumed that the newly formed sea ice is filling all leads, creating joints between the floes. The fast-growth FSD is therefore redistributed so that all ice is considered to be made of unbroken floes.

The growth of small floes resulting from wave-induced fragmentation in our model is also ensured by welding, which was shown by to generate a lateral growth rate one order of magnitude higher than that arising from the lateral accumulation of frazil ice. We, however, found the algorithm they used to be very dependent on the choice of the FSD categories. After some discussion with the authors of , we decided to carry on with this formulation but with an appropriate tuning of the coefficient that called $\mathit{\kappa }$, which represents the rate at which the number of floes decreases due to welding per surface area. We tune $\mathit{\kappa }$ so that the timescale at which a delta-function FSD made of the smallest floes allowed in the model evolves into a delta-function FSD made of the biggest possible floes is similar in our model and in the model by . To give an idea of the timescales involved, the model of , starting from a delta-function FSD only made of floes with an average size of 20 m, ends up with half the ice cover being made of floes larger than 200 m in about 5 d within compact sea ice ($c=\mathrm{0.95}$). In the simulations presented in this study, setting $\mathit{\kappa }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{8}}$ m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ reproduces a similar evolution of the FSD for our choice of FSD categories.

As mentioned earlier, the redistribution of the slow-growth FSD due to lateral growth in ${\mathrm{\Phi }}_{\text{th,slow}}$ is expected to happen on longer timescales that are related to the time needed by the fractures to heal. This healing phenomenon is related to the thickening and the consolidation of the joints, the formation of which is described by ${\mathrm{\Phi }}_{\text{th,fast}}$. It is very similar to the “damage healing” already included in neXtSIM see, and associated with a timescale ${\mathit{\tau }}_{\text{heal}}$, that we reuse in the computation of ${\mathrm{\Phi }}_{\text{th,slow}}$: 5$${\mathrm{\Phi }}_{\text{th,slow}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\text{heal}}}}\left(g_{\text{fast}}\left(D_{\text{fast}}\right)-g_{\text{slow}}\left(D_{\text{slow}}\right)\right)+{\mathrm{\Phi }}_{\text{conservation}},$$ where ${\mathrm{\Phi }}_{\text{conservation}}$ is an ad hoc term ensuring the conservation of sea ice area (i.e. Eq. ) for the slow-growth FSD when sea ice is created or melted. The slow-growth FSD therefore relaxes to the fast-growth FSD over a time ${\mathit{\tau }}_{\text{heal}}$, representing the (slow) strengthening of the joints between the floes. Note that this healing only occurs if the sea ice is exposed to freezing conditions. If new sea ice is formed, all new ice is added to the largest floe size category, as in the fast-growth FSD. If sea ice is melting, we assume that lateral melting has little effect on the shape of the FSD. In this case, ${\mathrm{\Phi }}_{\text{conservation}}$ is used to scale the slow-growth FSD without changing its shape.

2.2.2 Wave-induced sea ice fragmentation

In this section, we describe the implementation of the terms ${\mathrm{\Phi }}_{\text{m,slow}}$ and ${\mathrm{\Phi }}_{\text{m,fast}}$ that represent the mechanical redistribution of floes associated with the fragmentation of sea ice by waves in each FSD.

Similar to the work by , the occurrence of sea ice fragmentation in our coupled system is decided in the wave model. In WW3, sea ice breaks up if the wave curvature induces a stress that exceeds the sea ice flexural strength. The shortest wavelength for which the wave-induced stress exceeds the critical stress for flexural failure, which we call ${\mathit{\lambda }}_{\text{break}}$, is passed to neXtSIM (see Fig. ) and used in the mechanical redistribution scheme of the FSD. The determination of the value of ${\mathit{\lambda }}_{\text{break}}$ in WW3 was explained in detail in Sect. 2.3 of (where it is called ${\mathit{\lambda }}_{i,\text{break}}$). If no fragmentation has occurred in WW3, neXtSIM receives ${\mathit{\lambda }}_{\text{break}}=\mathrm{1000}$ m, which corresponds to the default unbroken value, and no fragmentation occurs in neXtSIM (resulting in ${\mathrm{\Phi }}_{\text{m,slow}}=\mathrm{0}$ and ${\mathrm{\Phi }}_{\text{m,fast}}=\mathrm{0}$). If neXtSIM receives a value of ${\mathit{\lambda }}_{\text{break}}<\mathrm{1000}$ m, then FSD redistribution occurs in neXtSIM. Fragmentation occurrence is then determined every coupling time step.

If fragmentation occurs, we assume that the thin ice joining aggregated floes is very likely to break, as reported by . This is where the memory of previous cracks stored in the slow-growth FSD plays a role. In the model, the quick failure of the cementing thin ice is represented by relaxing the fast-growth FSD, $g_{\text{fast}}$, to the slow-growth FSD, $g_{\text{slow}}$. In practice, we set ${\mathrm{\Phi }}_{\text{m,thermo}}\mathrm{\Delta }{t}_{\text{ice}}=g_{\text{fast}}-g_{\text{slow}}$, where $\mathrm{\Delta }{t}_{\text{ice}}$ is the ice model time step, therefore assuming that this relaxation is almost instantaneous. We justify this short relaxation time by the fact that (i) waves can fragment a consolidated sea ice cover in a few tens of minutes only and (ii) the fast-growth FSD $g_{\text{fast}}$ is only used for thermodynamical processes associated with timescales of at least a few hours and is therefore relatively unaffected by the choice of a relaxation time value one order of magnitude lower.

Once fragmentation has occurred, the relaxation of the fast-growth FSD to the slow-growth FSD gives $g_{\text{fast}}=g_{\text{slow}}$. This equality represents the fact that $D_{\text{fast}}=D_{\text{slow}}$ when fragmentation occurs and before thin ice starts cementing the floes. In the model, the shape of the two FSDs after a fragmentation event is then controlled by the mechanical redistribution occurring in the slow-growth FSD represented by the term ${\mathrm{\Phi }}_{\text{m,slow}}$. This term can be written in the same form as in :

6$$\begin{array}{rl}& {\mathrm{\Phi }}_{\text{m,slow}}=\\ & -Q\left(D\right)g_{\text{slow}}\left(D\right)+\underset{{\mathrm{0}}^{+}}{\overset{\mathrm{\infty }}{\int }}Q\left(D^{\prime }\right)\mathit{\beta }(D^{\prime },D)g_{\text{slow}}\left(D^{\prime }\right)\text{d}{D}^{\prime },\end{array}$$ where $Q\left(D\right)$ is a redistribution probability function characterizing which proportion of floes of a given size $D$ is broken (with $D$ representing indistinctively $D_{\text{fast}}$ and $D_{\text{slow}}$ during fragmentation) and $\mathit{\beta }(D^{\prime },D)$ is a redistribution factor quantifying the fraction of sea ice concentration transferred from floe size $D^{\prime }$ to $D$ as fragmentation occurs. The choices of $Q\left(D\right)$ and $\mathit{\beta }(D^{\prime },D)$ therefore shape the FSD resulting from wave-induced fragmentation. This shape is important as it strongly impacts processes involved in wave attenuation and lateral melting , but the evolution of the FSD during wave-induced fragmentation is still not well understood.

have suggested from field observations that the shape of the FSDs could be interpreted as two truncated power laws separated by a cut-off floe size and hypothesize that this cut-off corresponds to the critical floe size under which floes cannot be broken by waves. A cut-off floe size also seems to be visible in the FSDs that obtain from a laboratory experiment. The existence of a cut-off floe size is, however, contested. The use of cumulative distribution functions to interpret FSDs may give a false impression of scale invariance and the apparent change of regime could originate from finite measurement windows . The division of the FSD into two truncated power laws suggested by has nevertheless been used to redistribute FSDs in the wave-in-ice models by and , which have been reused for the wave-in-ice attenuation parameterization implemented in WW3 by and in many other studies using wave-in-ice models see, e.g..

In this study, the FSD is mostly used to provide WW3 with information on the floe size to estimate the wave attenuation. The FSD does not impact the amount of new ice formed and we focus on periods during which lateral melting can be neglected. In these conditions it is advantageous to stay close to what has been done already for the FSD in WW3 in order to ensure that wave attenuation is not too different from the one evaluated by and . We therefore build on the work by and to suggest a new parameterization for $Q\left(D\right)$ and $\mathit{\beta }(D^{\prime },D)$ that redistributes the FSD in a similar way to the assumptions made in wave-in-ice models derived from . However, there are two main differences from the FSDs assumed in : we allow the exponent of the power law corresponding to the “small-floe” regime to vary and we smooth the transition between the small-floe regime and the “large-floe” regime, as can be seen in the data from .

$Q\left(D\right)$ represents the amount of ice in each floe size category that will be broken. and a number of further studies using their approach used a step function with $Q\left(D\right)=\mathrm{1}$ if $D\ge \mathrm{0.5}{\mathit{\lambda }}_{\text{break}}$ and $D\ge D_{\text{FS}}$, where $D_{\text{FS}}$ is the minimum floe size for which flexural failure can occur see, and $Q\left(D\right)=\mathrm{0}$ if not. The transition of $Q\left(D\right)$ from 0 to 1 occurs at a critical floe size under which floes do not break, corresponding to what interpreted as the cut-off floe size. In the model of , the transition therefore occurs at $max(D_{\text{FS}},\mathrm{0.5}{\mathit{\lambda }}_{\text{break}})$. $D_{\text{FS}}$ depends on sea ice properties and is computed as 7$$D_{\text{FS}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{{\mathit{\pi }}^{\mathrm{4}}Y{h}^{\mathrm{3}}}{\mathrm{48}\mathit{\rho }g(\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}})}}\right)}^{\mathrm{1}/\mathrm{4}},$$ where $g$ is gravity, $Y$ is the Young modulus of sea ice, $\mathit{\nu }$ is Poisson's ratio, $h$ is the mean sea ice thickness and $\mathit{\rho }$ is the density of sea water. For ice thinner than 1 m, which is often the case in the MIZ, $D_{\text{FS}}$ is lower than 15 m and the cut-off floe size is likely to be determined by the value of ${\mathit{\lambda }}_{\text{break}}$. To define $Q\left(D\right)$, we take an approach similar to and set 8$${\displaystyle }Q\left(D\right){\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\text{WF}}}}{p}_{\text{FS}}(D,D_{\text{FS}})p_{\mathit{\lambda }}(D,{\mathit{\lambda }}_{\text{break}}),$$ in which ${\mathit{\tau }}_{\text{WF}}$ is a relaxation time associated with wave-induced fragmentation events, and $p_{\text{FS}}$ and $p_{\mathit{\lambda }}$ are probabilities that floes break depending on their size. We introduced ${\mathit{\tau }}_{\text{WF}}$ to avoid dependency of the FSD redistribution to the coupling time step. It represents the timescale needed for the FSD of a fragmenting sea ice cover to reach a new equilibrium under a constant sea state. We set it to 30 min, as this corresponds to the timescale of the fragmentation event described in . The probability functions $p_{\text{FS}}$ and $p_{\mathit{\lambda }}$ express the idea that the smaller the floes are, the less chance they have to break. The function $p_{\mathit{\lambda }}$ compares the floe size $D$ with the value of $D_{\text{FS}}$ and $p_{\mathit{\lambda }}$ compares the floe size $D$ with ${\mathit{\lambda }}_{\text{break}}$, introducing a dependency of $Q\left(D\right)$ on the wave field. The difference with the model of is that instead of step functions we use hyperbolic tangents to get a continuous transition of $Q\left(D\right)$ between 0 and 1: $$\begin{array}{}\text{9a}& {\displaystyle }& {\displaystyle }{p}_{\text{FS}}\left(D\right)=max\left(\mathrm{0},\tanh \left({\displaystyle \frac{D-c_{\mathrm{1},\text{FS}}{D}_{\text{FS}}}{c_{\mathrm{2},\text{FS}}{D}_{\text{FS}}}}\right)\right),\text{9b}& {\displaystyle }& {\displaystyle }{p}_{\mathit{\lambda }}(D,{\mathit{\lambda }}_{\text{break}})=max\left(\mathrm{0},\tanh \left({\displaystyle \frac{D-c_{\mathrm{1},\mathit{\lambda }}{\mathit{\lambda }}_{\text{break}}}{c_{\mathrm{2},\mathit{\lambda }}{\mathit{\lambda }}_{\text{break}}}}\right)\right),\end{array}$$ in which $c_{\mathrm{1},\text{FS}}$, $c_{\mathrm{2},\text{FS}}$, $c_{\mathrm{1},\mathit{\lambda }}$ and $c_{\mathrm{2},\mathit{\lambda }}$ are parameters of the FSD that control the range of floe size that will be broken or not. The use of a continuous $Q\left(D\right)$ instead of a step function aims to relax the constraint on the FSD shape imposed by . With a step function, the probability of having floes larger than the cut-off floe size is 0, i.e. above the cut-off floe size the FSD is suddenly infinitely steep. This approach is particularly problematic. First, the FSDs reported in show a gradual steepening rather than a sharp transition. Second, the steepening of the FSD slope that led to the identification of the two floe size regimes by could actually be due to windowing issues . Here, instead of having a single cut-off floe size, we have a transition occurring in the floe size range for which $\mathrm{0}<Q\left(D\right)<\mathrm{1}$. The width of this range is controlled by $c_{\mathrm{2},\text{FS}}$ and $c_{\mathrm{2},\mathit{\lambda }}$. $Q\left(D\right)$ tends towards a step function when $c_{\mathrm{2},\text{FS}}$ and $c_{\mathrm{2},\mathit{\lambda }}$ tend towards 0. We found that setting $c_{\mathrm{2},\text{FS}}=\mathrm{2}$ and $c_{\mathrm{2},\mathit{\lambda }}=\mathrm{2}$ leads to FSDs with a gradual steepening that look very similar to what is reported by or in the lab experiments by , for instance. The location of the floe size range at which the transition occurs is controlled by $c_{\mathrm{1},\text{FS}}$ and $c_{\mathrm{1},\mathit{\lambda }}$. Like , we set $c_{\mathrm{1},\text{FS}}=\mathrm{1}$. Instead of using $c_{\mathrm{1},\mathit{\lambda }}=\mathrm{0.5}$ as in , we set $c_{\mathrm{1},\mathit{\lambda }}=\mathrm{0.3}$. It is consistent with the hypothesis made by that floes smaller than $\mathrm{0.3}{\mathit{\lambda }}_{\text{break}}$ are tilted by waves but do not bend and therefore have no chance of suffering flexural failure. Reducing the value of $c_{\mathrm{1},\mathit{\lambda }}$ reduces the value of the cut-off floe size we should get, which generally leads to floes smaller than assumed in . It is, however, compensated for by using a continuous $Q\left(D\right)$, which gives more weight to large floes than the FSD assumed in the model of . FSDs generated by this redistribution function are presented and discussed in Sect. .

The redistribution factor $\mathit{\beta }$ controls the shape of the FSD after it has been redistributed. Similar to the work by and , the redistribution factor $\mathit{\beta }$ follows the form 10$$\left\{\begin{array}{rl}{\displaystyle }& {\displaystyle }\mathit{\beta }(D_n,D)={\displaystyle \frac{D_n^q-D_{n-\mathrm{1}}^q}{D^q-D_{\mathrm{0}}^q}}\text{if}{D}_n\ge D\\ {\displaystyle }& {\displaystyle }\mathit{\beta }(D_n,D)=\mathrm{0}\text{otherwise }\end{array}\right.,$$ where $n$ corresponds to the index of the $n$th FSD category and $q$ is an exponent that controls the shape of the redistributed FSD. used $q=\mathrm{1}$ arguing that since the fragmentation of floes is a stochastic process, any floe size lower than the initial unbroken floe can be generated with a similar probability. It leads to power-law FSDs after a succession of fragmentation events. noted that $q$ is related to the (negative) exponent $\mathit{\alpha }$ of the power-law FSD associated with the smaller floe sizes resulting from the redistribution as assumed by and later by with $q=\mathrm{2}+\mathit{\alpha }$. In their study, related the exponent $\mathit{\alpha }$ to a quantity they call fragility, which represents the probability of fragmentation of floes. assumed that fragility is constant and equal to 0.9, giving $\mathit{\alpha }\simeq -\mathrm{1.85}$. set the value $q$ to prescribe power-law FSDs with this same value of $\mathit{\alpha }$ but remain consistent with what was done in wave models before. This is, however, a big constraint on the shape of the FSD, and it contradicts the variations of the exponent of the power law fitted to the small-floe regime reported by . Here, we have already expressed the probability of fragmentation of floes, hence the fragility, as the product of $p_{\text{FS}}$ and $p_{\mathit{\lambda }}$. Using this relationship, we get 11$$q=-{\log }_{\mathrm{2}}\left(p_{\text{FS}}\left(D\right)p_{\mathit{\lambda }}(D,{\mathit{\lambda }}_{\text{break}})\right).$$ This definition of $q$ still generates FSDs that tend to a power law for floe sizes lower than the cut-off floe size, but it does not prescribe a fixed value to the exponent $\mathit{\alpha }$ of this power law. Instead, $\mathit{\alpha }$ decreases as the ice thickens and wavelength increases as is shown and discussed in Sect. .

The processes we use for the wave-in-ice attenuation computation in WW3 require the estimation of two floe size parameters: the average floe size $⟨D⟩$ and the maximum floe size $D_{max}$ . When WW3 is run in stand-alone mode, $D_{max}$ is taken to be ${\mathit{\lambda }}_{\text{break}}/\mathrm{2}$ and an assumption is made on the shape of the FSD after fragmentation allows us to estimate the value of $⟨D⟩$. When WW3 is coupled with neXtSIM, the FSD is free to evolve to the sea ice model and it is necessary to estimate the value of $D_{max}$ and $⟨D⟩$ in neXtSIM so it can be sent to WW3.

The average floe size $⟨D⟩$ can be simply defined as 12$$⟨D⟩=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}Dg\left(D\right)\text{d}D.$$ Here we use the slow-growth FSD, $g_{\text{slow }}\left(D_{\text{slow}}\right)$, assuming that wave scattering, which is the wave attenuation process depending on $⟨D⟩$, is more affected by consolidated floes than by the thin ice joining them. The maximum floe size, $D_{max}$, definition is less straightforward, as it was originally designed in the FSD parameterization of to represent the largest floe size of a fragmented sea ice cover, and if no fragmentation had occurred it was set to a large default value. Here, this definition needs to be extended to a coupled system with an FSD free to evolve under the effects of both mechanical and thermodynamical processes and able to represent a mix of fragmented floes and large ice plates. We suggest a definition based on the percentage of the ice cover area occupied by large floes, computing $D_{max}$ as the $\mathrm{90}$th percentile of the areal FSD. Besides, the flexure dissipation mechanisms included in WW3 by require us to discriminate between sea ice cover made of large floes with size on the order of $O\left(\mathrm{100}\right)$ m and an unbroken sea ice cover for which the default $D_{max}$ in WW3 is set to 1000 m. This is because flexure only occurs if the wavelength is shorter or of the same order as the floe size. Knowing that long swells can have wavelengths on the order of $O\left(\mathrm{100}\right)$m, they are only fully attenuated by inelastic dissipation if floe size is on the order of $O\left(\mathrm{1000}\right)$ m, which can be larger than the floe size range covered by the FSD defined in neXtSIM. In the case where $D_N<\mathrm{1000}$ m, to make sure that swells are still attenuated in an unbroken sea ice cover by WW3, we linearly increase the value of $D_{max}$ sent to WW3 from $D_{max}=D_N$ to $D_{max}=\mathrm{1000}$ m with the proportion of sea ice in the largest floe size category ${\int }_{D_{N-\mathrm{1}}}^{D_N}{g}_{\text{slow}}\left(D\right)\text{d}D/c$.

2.3 Link between wave-induced sea ice fragmentation and damage

As mentioned in the Introduction, it is expected that sea ice fragmentation by waves results in lowering the ice internal stress. The lowering of sea ice resistance to deformation due to a high density of cracks is already represented in neXtSIM by the variable called damage. This variable takes a value between 0 and 1, with 0 corresponding to undamaged sea ice, and 1 to a highly damaged sea ice cover, i.e. presenting a high density of cracks. Sea ice damaging in neXtSIM is usually due to the wind. In our study, we want it to have an additional dependence on wave-induced sea ice fragmentation. In our implementation, it is possible to quantify the sea ice cover area that is susceptible to being broken by waves if WW3 provides a value of ${\mathit{\lambda }}_{\text{break}}<\mathrm{1000}$ m as $c_{\mathrm{broken}}=c\left(\mathrm{1}-{\mathrm{e}}^{-\mathrm{\Delta }{t}_{\text{cpl}}/{\mathit{\tau }}_{\text{WF}}}\right)$ and $c_{\mathrm{broken}}=\mathrm{0}$ if ${\mathit{\lambda }}_{\text{break}}\ge \mathrm{1000}$ m. As the fragmentation of sea ice by waves can break ice plates into floes with sizes up to a few hundred metres and the horizontal resolution of the model mesh is in general at least a few kilometres, we hypothesize that areas of the sea ice cover fragmented by waves are associated with high values of damage, i.e. close to 1. We thus suggest computing the new damage value associating a value of $d_w=\mathrm{0.99}$ to $c_{\mathrm{broken}}$, which gives the following evolution of the damage $d$:

13$$d=min\left(\mathrm{1},d(\mathrm{1}-c_{\mathrm{broken}})+d_w{c}_{\mathrm{broken}}\right).$$ This process is repeated every time fragmentation occurs in the sea ice model. Note that because wave events generally last for a few hours, this damaging process is generally repeated enough times to result in little sensitivity of the model to values of $d_w$ between 0.1 and 1. Note also that floe size and damage are not explicitly linked by this relationship, but the relaxation time associated with the healing of damage and the slow-growth FSD are the same, making their evolution parallel in the regions of broken ice.

3 Model set-up

3.1 General description of the pan-Arctic configuration

Similarly to , the coupled framework is run on a regional 0.25${}^{\circ }$ grid (CREG025), which covers the Arctic Ocean at an approximate resolution of 12 km as well as some of the North Atlantic. As neXtSIM is a finite-element sea ice model using a moving Lagrangian mesh, it is not run on a grid. Its initial mesh is, however, based on a triangulation of the CREG025 grid, giving a prescribed mean resolution (i.e. mean length of the edges of the triangular elements) of 12 km. In coupled mode, neXtSIM interpolates the fields to be exchanged onto the fixed grid so that the coupler, OASIS, is only required to send and receive different fields (i.e. OASIS does not need to do any interpolation).

neXtSIM is run with a time step of 20 s, and WW3 with a time step of 800 s. Fields between the two models are exchanged every 2400 s. Atmospheric forcings are provided by 6-hourly fields from the CFSv2 atmospheric re-analysis . In addition, neXtSIM is also forced by ocean fields from the TOPAZ4 re-analysis . For more details on the forcings used by neXtSIM, see . Wave–current interactions in WW3 are not considered in this study.

For the FSD in neXtSIM, we use 20 categories with a width $\mathrm{\Delta }D=\mathrm{10}$ m. The lower bound $D_{\mathrm{0}}$ is set to 10 m, which is about the size of the smallest floes susceptible to undergo flexural failure . The upper bound of the FSD, $D_N$, is therefore equal to 210 m, which is on the order of magnitude of the largest floes resulting from wave-induced sea ice fragmentation generated by WW3. The healing relaxation time ${\mathit{\tau }}_{\text{heal}}$ used by neXtSIM is set to its default value of 25 d .

3.2 Evaluation of the FSD implementation: model set-up

In Sect. , we use sea state and sea ice observations made in the Beaufort Sea in the framework of the Arctic Sea State and Boundary Layer Physics Program to evaluate the wave attenuation and broken sea ice extent in our coupled simulations, which is similar to what did with stand-alone WW3 simulations. To do so, we ran five simulations from 10–13 October 2015, a period covering the storm event investigated in a study by . This storm generated $\simeq$ 4 m waves fragmenting the sea ice edge in the Beaufort Sea from 11–12 October 2015.

The first of these simulations was a WW3 uncoupled simulation hereafter labelled ARD18, as it used the exact same parameterization as the one labelled REF2 in . The only difference is that here it was run on the CREG025 grid used for all our simulations. This parameterization of the wave model was chosen as it showed the best match with observations for both wave height and broken sea ice extent in the study by . We also used the same sea ice concentration data as to force the uncoupled wave model. They were obtained from a re-analysis of the 3 km resolution sea ice concentration dataset derived from the AMRS2 radiometer using the ASI algorithm

AMSR2 data are available at https://seaice.uni-bremen.de/data/amsr2/asi_daygrid_swath/n3125/2015/oct/Arctic3125/ (last access: November 2020).

. This re-analysis produced 12-hourly maps that gathered all the AMSR2 passes acquired between 00:00 and 13:59 UTC and 10:00 and 23:59 UTC for the morning (AM) and evening (PM) fields, respectively. Similarly to , sea ice thickness was set constant to 15 cm. Sea ice concentration was also kept constant to make the comparison with a coupled simulation easier. Sea ice concentration for this 3 d run corresponded to the conditions on the evening of 12 October provided by the AMSR-2 sea ice concentration re-analysis at the same time as illustrated in study. Initial wave conditions were provided by an initial 10 d run of this simulation, from 1–10 October 2015, in which sea ice concentration is updated every 12 h.

Secondly, we ran a coupled neXtSIM–WW3 simulation (hereafter labelled NXM/WW3). Initial conditions were the same as in ARD18. Sea ice dynamics and thermodynamics were switched off in neXtSIM so that we can compare the two simulations with a similar constant sea ice cover (thickness and concentration). The only difference between ARD18 and NXM/WW3 was the way the evolution of floe size was treated, which is what we wanted to evaluate.

Then, to illustrate the sensitivity of our results to the sea ice thickness value, we re-ran ARD18 and NXM/WW3 while this time setting the sea ice thickness to 30 cm. We named these two additional simulations ARD18_H30 and NXM/WW3_H30, respectively.

Finally, to investigate the sensitivity of the FSD evolution to the floe size categories used for the FSD, we ran a simulation similar to NXM/WW3 but with a refined FSD that we call NXM/WW3_refine. For this simulation, the number $N$ of categories was set to 41 instead of 20 and we set $\mathrm{\Delta }D=\mathrm{5}$ m and $D_{\mathrm{0}}=\mathrm{5}$ m. We also evaluated the evolution of the two FSDs with refreezing/healing using the CPL_DMG simulation described in Sect. .

3.3 Estimation of the impact of wave-induced fragmentation on sea ice dynamics: model set-up

In Sect. , we compare the results of three simulations in order to investigate the wave impact on sea ice dynamics in the MIZ. The first one (called REF) is a stand-alone simulation of neXtSIM. The second one (CPL_WRS) includes all the features presented in Sect.  but uses the relationship between wave-induced sea ice fragmentation and damage presented in Sect. . The third simulation (CPL_DMG) is similar to CPL_WRS except that it also includes a link between the damage variable $d$ and wave-induced sea ice fragmentation as described in Sect. . These simulations were run from 15 September to 1 November 2015. This period was selected because refreezing occurs in the MIZ, meaning that the differences between REF and the two coupled simulations were not due to the change in lateral melting parameterization. This period of the year is also characterized by the combination of a low sea ice extent (thus a large available fetch) and regular occurrence of storms in the Arctic, which increases the opportunities to evaluate the impact of waves on sea ice with fragmentation events over wide areas. The level of damage in the ice cover was initially set to zero where sea ice is present. Initial sea ice concentration and thickness were set from the TOPAZ4 re-analysis and sea ice is unbroken. The wave field in WW3 was initially at rest. The wave-in-ice attenuation parameterization in WW3 in CPL_WRS and CPL_DMG was the same as in ARD18 i.e. REF2 in. We investigated the results of these simulations from 1 October, thus allowing for 16 d of spin-up, which is enough for the wave and damage fields to develop.

Figure 3

Spatial distributions of maximum floe size (a) and significant wave height (b) in the Beaufort Sea taken from the NXM/WW3 simulation on 12 October at 17:00:00 GMT. Black arrows indicate the wave mean direction. The difference of the maximum floe size distribution with the ARD18 simulation is shown in panel (c). Evolution of the significant wave height for different simulations along the transect depicted in cyan in panels (a)–(c) is presented in panel (d), along with significant wave height estimated from Sentinel-1a SAR images seefor details and measured by an AWAC buoy. The AWAC position is depicted by a red square in panels (a)–(c). The green and magenta crosses indicate the position at which the FSDs are shown in Fig. . Solid and dashed black lines represent contours of sea ice concentration equal to 0.8 and 0.15, respectively. Point Barrow is indicated by a black triangle.

[Figure omitted. See PDF]

4 Results

4.1 Evaluation of wave–sea-ice interactions in the coupled framework

We first evaluate the representation of wave–ice interactions in our coupled framework. As our goal is to investigate the potential impact of waves on sea ice dynamics, we must ensure that the coupled framework produces a consistent wave-in-ice attenuation, as it is directly proportional to the WRS, as well as reasonable extents of broken sea ice and timescales of the ice recovery from fragmentation. In the following, we consider the wave attenuation and extent of fragmented sea ice on one hand, and the evolution of the FSDs after fragmentation events on the other hand.

4.1.1 Evaluation of wave attenuation and extent of fragmented sea ice

To evaluate the capacity of our coupled framework to produce reasonable wave-in-ice attenuation and extent of broken sea ice, we focused on the same event used to evaluate the WW3 parameterization of in the study by . Here, we used the ARD18 simulation as a reference, as it was shown to match well with observations for both the extent of broken ice and the wave attenuation in this particular case. The comparison was done on 12 October at 17:00 GMT (Fig. ). We also compared our model results with estimated wave height from SAR images and buoy measurements AWAC, see along a transect in Fig. d.

We evaluated the wave attenuation by looking at the evolution of the maximum floe size $D_{max}$ and significant wave height in the ice. Overall, the spatial distribution of these two quantities in ARD18 and NXM/WW3 was very similar and also very similar to the results of (see Fig. 9 of their study). Figure d shows the wave height evolution along a transect following the footprint of Sentinel 1-a; again, we see almost no difference between ARD18 and NXM/WW3. Both simulations showed reasonable agreement with the wave heights estimated from SAR and from the AWAC buoy. Similar to the results of , the model, however, seemed to slightly overestimate the wave height within the ice cover. This overestimation could have resulted from the assumption of constant thickness and low value (15 cm). This is visible in Fig. d where most observations actually show higher significant wave height values than the one yielded by ARD18_H30 and NXM/WW3_H30, which use a constant thickness of 30 cm.

Sea ice break-up occurrence depending on wave properties, comparable wave attenuation in between ARD18 and NXM/WW3 resulted in little difference in the extent of broken sea ice between the two simulations (Fig. a, c). Although the extent of broken ice was slightly smaller in the coupled run, the difference did not exceed two grid cells and therefore represented a distance of about 25 km, which is acceptable given the large uncertainties associated with wave attenuation in ice see, for instance,. Moreover, as in , the broken sea ice region extended up to about 15 km in-ice beyond the AWAC buoy (red square), which matched well with the SAR images observed for that same day.

4.1.2 Evaluation of the evolution of the FSDs

Our coupled framework introduces two FSDs to represent the evolution of the floe size from two different points of view. It also introduces a new redistribution scheme used when wave-induced sea ice fragmentation occurs. This section provides a brief evaluation of these new features.

We first look at the FSD resulting from wave-induced fragmentation in neXtSIM by plotting the cumulative distribution of floes CDF; see, e.g. for three different locations (Fig. ) for the NXM/WW3 and NXM/WW3_refine simulations. Note that because thermodynamical and dynamical processes are unactivated in NXM/WW3, the fast- and slow-growth FSDs are identical. The CDFs look very similar to the one reported by, e.g. , with the curve gradually steepening as floe size increases. Following the method of , two lines can be fitted to the FSD for small and large floes, defining two regimes intersecting at a cut-off floe size: (i) a small-floe regime that follows a power law and (ii) a large-floe regime with a much steeper slope of the CDF. This interpretation may result from an artefact arising from the use of the CDFs and finite-size windows . The use of CDFs can give a false impression of scale-invariance but has nonetheless been used in a number of wave–ice interaction studies, in particular the one used to calibrate the wave attenuation model we use here and in all the models using the work of . Here we use the CDFs to discuss the differences introduced by our parameterization of the FSD redistribution compared to the assumptions made by and discuss how they could impact wave attenuation.

Figure 4

Cumulative distribution of floes taken at three different locations from the NXM/WW3 simulation (a) and from the NXM/WW3_refine simulation (b). The three locations share the same longitude (150${}^{\circ }$ W), and their position is indicated by a symbol of the same colour in Fig. a and b (red colour corresponds to the AWAC position). The dashed lines correspond to the linear regression over the smallest floe size categories for each location. Values of the slope are given in the legend. The vertical dotted lines represent the $D_{\text{max}}$ values for each case.

[Figure omitted. See PDF]

As in the study by , the cut-off floe size and the (negative) power-law exponent of the small-floe regime increase with the distance from the ice edge. For the two simulations, all but one value of the (negative) power-law exponents related to the small-floe regime are greater than $-\mathrm{2}$, as expected for a CDF that depends on a two-dimensional fragmentation process . The one value of the exponent that does not lie in this range is obtained for the NXM/WW3 run close to the ice edge, where the cut-off floe size is too close to the value of $D_{\mathrm{0}}$ for the small-floe size regime to be resolved. In the model of , the exponent of the power law is equal to $\simeq -\mathrm{1.85}$. Close to the ice edge, our FSDs and the one assumed by models derived by are therefore likely to be very similar. However, further away in the ice, our parameterization gives less weight to small floes in the FSDs. For wave attenuation, this means that there is less scattering occurring as waves propagate towards pack ice. Little impact on wave attenuation is expected, as scattering is mostly efficient for short waves, which do no propagate far into the ice. $D_{max}$ in the model of corresponds to the “cut-off” floe size as interpreted by . Here, the values of $D_{max}$ for each location lie in the floe size range corresponding to the transition between the two regimes, agreeing with the definition used by . $D_{max}$ shows little sensitivity to the FSD definition, with a maximum difference between the two simulations presented in Fig.  not exceeding the value of $\mathrm{\Delta }D$ in the refined FSD simulation (5 m). Considering the large uncertainties due to the limited knowledge of wave–ice interactions, choosing $\mathrm{\Delta }D=\mathrm{10}$ m instead of more refined FSDs in our coupled framework therefore has little impact on the wave attenuation computed in WW3.

Figure 5

Pan-Arctic distribution of the area covered by floes with diameters larger than 200 m over the total sea ice cover area according to the fast-growth (a, b) FSD, slow-growth FSD (c, d), and distribution of the maximum floe size (e, f). Each column corresponds to a different time: 2 October 2015 00:00:00 (a, c, e) and 1 November 2015 00:00:00 (b, d, f). The red cross (a, c, e) indicates the location at which the FSDs shown in Fig.  and the floe size parameters in Fig.  are taken.

[Figure omitted. See PDF]

Figure 6

Areal (a, c) and cumulative distribution of floes (b, d) for the fast-growth (a, b) and slow-growth (c, d) FSDs taken at 45${}^{\circ }$ E and 83.5${}^{\circ }$ N. Each colour refers to a different time: 2 October 2015 00:00:00 (red) and 3 October 2015 00:00:00 (blue). In panel (a, c), the bar at 200$+$ m represents unbroken ice. The dashed lines correspond to the linear regression over the smallest floe size categories for each location. Values of the slope $\mathit{\alpha }$ are given in the legend. The vertical dotted lines represents the $D_{\text{max}}$ values for each date.

[Figure omitted. See PDF]

As the simulations in this study focused on the autumn period, when the sea ice cover expands due to freezing, we also check that the timescales associated with freezing and sea ice healing are reasonable. Figures  and  illustrate the evolution of the FSDs in the CPL_DMG simulation, in which both mechanical and thermodynamical processes are active. Figure a and c show the proportion of “unbroken” ice (the proportion of sea ice associated with the $N$th category of the FSD) for the fast- and slow-growth FSDs, respectively, 17 d after the beginning of the simulation, leaving enough time for the waves and sea ice to spin-up. The regions of broken ice are relatively similar for both FSDs with the exception of the Barents Sea area. They actually closely follow the contour of 1 m thick ice (not shown) after which the waves are too attenuated to fragment the sea ice. Floe size generally increases with distance from the ice edge as the ice gets thicker and shorter waves are quickly attenuated (Fig. e, f). The Barents Sea area shows a wide area of broken thin ice in the slow-growth FSD with only parts of this region being broken in the fast-growth FSD. This wide broken area is related to a strong wave event in this region that occurred between 30 September and 1 October 2015. This event is associated with wave height up to 9 m and waves with periods above 12 s propagating far into an ice cover made of relatively thin ice (less than 1 m). About 24 h after the event (not shown), the fast-growth FSD in pack ice has mostly “recovered” due to welding and freezing in leads. In pack ice, where floes are larger than at the ice edge, the speed of the floe size growth in the fast-growth FSD is mostly controlled by welding and therefore depends on the value chosen for the rate at which the number of floes decreases, $\mathit{\kappa }$. This is because, like , we use a constant value for $\mathit{\kappa }$, meaning that the fewer floes there are (i.e. the larger the floe size), the higher the proportion of floes that merge during a given time period. The growth of large floes in the slow-growth FSD takes much longer, with a timescale set by the value of ${\mathit{\tau }}_{\text{heal}}$ (25 d in CPL_DMG), and the end of the mechanical healing of the Barents Sea area was still visible on 1 November 2015 (Fig. d).

This difference in timescales is also visible in Fig. , which illustrates the evolution of the fast- and slow-growth FSDs (panels a, c) and CDFs (panels b, d) for the location indicated by a cross in Fig. a and c at the time shown on the snapshot (2 October 2015 at 00:00:00 GMT) and 24 h later. On 2 October 2015 the proportion of sea ice that was unbroken was almost 0 in the slow-growth FSD, while welding and refreezing in the fast-growth FSD allowed the re-formation of unbroken sea ice over more than 10 % of the ice-covered part of the mesh element. The action of welding and refreezing in the fast-growth FSD results in a steepening of the slope of the CDF associated with the small-floe regime and flattening of the slope of the large-floe regime. The slow-growth FSD shows no sign of any healing, the last fragmentation event being too recent, and its associated CDF clearly shows a small- and large-floe regime resulting from the fragmentation by waves (Fig. ). More than 95 % of the fast-growth FSD consists of unbroken sea ice 24 h later, while the slow-growth FSD is still very similar to what it was on the previous day, illustrating the memory effect of the slow-growth FSD.

Figure 7

Temporal evolution of the two floe size parameters computed in the ice model and sent to the wave model taken at 45${}^{\circ }$ E and 83.5${}^{\circ }$ N (see red cross in Fig. a, c, e).

[Figure omitted. See PDF]

The evolution of the slow-growth FSD involves the use of an ad hoc term to ensure the conservation of sea ice area as well as a relaxation towards the fast FSD (see Eq. ). These terms might affect the moments of the distribution in a non-trivial way, as was shown in Sect. 3.1 of . This might affect the estimations $D_{max}$ and $⟨D⟩$, which are both computed from moments of the slow-growth FSD. In the absence of observations that could help to constrain our model, we represent in Fig.  the evolution of $D_{max}$ and $⟨D⟩$ for the location indicated in Fig. a, c and e. This evolution is in line with what we could expect from the floe size evolution with sudden drops due to fragmentation events and a slower growth over a timescale of about 20 d.

In summary, once the wave activity has decreased, refreezing and welding allow for the re-generation of a completely unbroken sea ice cover on timescales of a few hours to a few days in the fast-growth FSD, depending on the initial level of fragmentation. The timescale over which the slow-growth FSD re-generates large floes is associated with the value of ${\mathit{\tau }}_{\text{heal}}$. As an indication of time, in the case illustrated here with freezing conditions and compact ice, it takes about 4 d for the value of $D_{max}$ to grow over 200 m.

Figure 8

Distributions of sea ice thickness (a), maximum floe size (b) and sea ice concentration (c) from the CPL_DMG simulation in the Barents Sea on 15 October 2015 at 00:00:00. The magenta line corresponds to the contour $D_{max}=\mathrm{200}$ m in CPL_DMG. The black thick lines delimit the domain used to analyse the results from the three simulations (REF, CPL_WRS, CPL_DMG) in Figs. , and .

[Figure omitted. See PDF]

4.2 Impact of sea ice fragmentation on sea ice dynamics in the MIZ

4.2.1 Case study: a fragmentation event in the Barents Sea (15–25 October 2015)

To better understand the impact of (i) the WRS and (ii) fragmentation-induced damage on sea ice dynamics, we compare the results given by the REF, CPL_WRS and CPL_DMG simulations in the Barents Sea. Focusing only on this region simplifies the analysis, as this area is exposed to wave and sea ice conditions that experience little variation over the investigated period. The available fetch, in particular, remains relatively constant and is large enough to allow for storm waves to penetrate far into the ice. The sea ice edge also remains oriented mostly east–west over this period. In our analysis, we can therefore consider southward winds to be mostly off-ice and conversely northward winds to be directed on-ice. The domain we define to perform our analysis is limited south and north by the 69 and 84${}^{\circ }$ N parallels, respectively, and west and east by the 16 and 60${}^{\circ }$ E meridians (see, for instance, Fig. ).

To highlight the various responses of sea ice to fragmentation depending on wind and waves conditions, we select a particular fragmentation event occurring on 15 October 2015 (see Fig.  for the initial sea ice conditions), which results in a growth of the surface occupied by broken ice (Fig. c). The 10 d following this event include both on-ice and off-ice conditions, allowing us to explore the impact of wave-induced sea ice fragmentation on sea ice dynamics in both cases.

Figure 9

Temporal evolution of (a) the wind speed (black solid line) and its meridional component (dashed grey line) averaged over the entire ice-covered part of the domain (see Fig.  for the domain definition), (b) the significant wave height averaged over the entire domain (black solid line) and ice-covered part of the domain only (dashed grey line) and (c) the ratio of the surface area of regions covered by recently broken sea ice (defined as $D_{max}\le \mathrm{200}$ m, black solid line) and compact sea ice that has been recently broken (defined as $D_{max}\le \mathrm{200}$ m and $c\ge \mathrm{0.8}$, grey dashed line) over the total sea ice-covered surface area. The time period shown covers 15–25 October 2015, for which initial conditions are given in Fig. . The sea ice-covered part of the domain is the area for which the sea ice concentration $c$ is greater than 0. The two orange vertical lines indicate the dates of the snapshots shown in Figs. , and .

[Figure omitted. See PDF]

The fragmentation event occurrs during on-ice wind conditions (see the positive meridional component in Fig. a) with high waves (up to 5 m at the ice edge) propagating far into the ice cover (see Figs. b and b). It results in an increase of the surface area made of recently broken sea ice in the domain (see, for instance, the evolution of the magenta contour between Figs. a and a) until it represents nearly 60 % of the sea ice-covered part of the domain (Fig. c). Here we define the sea ice-covered part of the domain as the area for which sea ice concentration $c>\mathrm{0}$. Recently broken ice is defined as the region for which $D_{max}\le \mathrm{200}$ m, thus corresponding to fragmented sea ice for which refreezing has not yet had time to heal a significant proportion of unbroken ice (at least 10 %). In the domain, the recently broken sea ice area is mostly made of compact sea ice (that we define as the area of the domain for which $c>\mathrm{0.8}$). At the end of the fragmentation event, 80 % of the recently broken sea ice is made of compact ice and represents nearly 40 % of the ice-covered part of the domain (Fig. c). Compact sea ice being broken is important in the scope of our study, as low-concentration sea ice experiences little resistance to deformation due to its low effective stiffness and is therefore largely unaffected by damage.

In the days following the fragmentation event, wind speed decreases and changes direction to become mostly southward (Fig. a). The wind speed then increases again, with a second maximum on 20 October, corresponding this time to off-ice winds. The generated waves tend to propagate away from the sea ice, resulting in low wave heights inside the ice (Fig. b). The increase in the wind speed coincides with a decrease in recently broken sea ice surface area, which stops when the wind turns again to become parallel to the sea ice edge (Fig. c). This decrease is mostly due to sea ice recovering in the absence of waves and can be seen by comparing the distributions of sea ice thickness and damage on 16 October (Figs. a and a) and 21 October (Figs. c and c), in which the limit between broken and unbroken ice tends to get closer to the ice edge, while damage in pack ice visibly decreases. The band of recently broken ice remains, however, much larger than it was initially (Fig. b), as fragmented sea ice produces lower wave attenuation, thus allowing for sea ice fragmentation even in low wave height conditions.

Figure 10

Snapshots of the distributions of sea ice thickness (a, c) and of the ratio between the WRS and the wind stress over sea ice (b, d) from the CPL_DMG simulation in the Barents Sea taken on 16 October 2015 at 18:00:00 GMT (a, b) and on 21 October 2015 at 09:00:00 GMT (c, d). The wind and WRS directions for each date are given by the red and blue arrows, respectively, in panels (b) and (d). The magenta line corresponds to the contour $D_{max}=\mathrm{200}$ m in CPL_DMG. The black thick lines delimit the domain used to analyse the results from the three simulations (REF, CPL_WRS and CPL_DMG) in Figs. , and .

[Figure omitted. See PDF]

4.2.2 Effects of linking damage and sea ice fragmentation on sea ice dynamics in the MIZ

The conditions of the described study period, the impact of adding the WRS, and a relationship between wave-induced fragmentation and damage can be investigated. We first proceed by comparing (in Fig. ) the ice drift velocity averaged over the ice-covered domain for the three simulations: REF, CPL_WRS and CPL_DMG. Overall, we observe no differences in the trends; however, the magnitudes of the ice drift velocities show intermittent differences between these three simulations. These differences have two maxima: the first on 16 October at 18:00:00 GMT and the second on 21 October at 09:00:00 GMT. To understand these differences, we compare the CPL_DMG and REF simulations on these two dates.

Figure 11

Snapshots of the distributions of sea ice damage in the CPL_DMG simulation (a, c) and of the differences in damage between the CPL_DMG and REF simulations (b, d) in the Barents Sea taken on 16 October 2015 at 18:00:00 GMT (a, b) and on 21 October 2015 at 09:00:00 GMT (c, d). The magenta line corresponds to the contour $D_{max}=\mathrm{200}$ m in CPL_DMG. The black thick lines delimit the domain used to analyse the results from the three simulations (REF, CPL_WRS and CPL_DMG) in Figs. , and .

[Figure omitted. See PDF]

On 16 October, the wind and waves are directed on-ice, thus compacting the sea ice (Fig. b). At the time of the snapshot, sea ice has been recently broken over a wide surface area (within the magenta contour in Figs. a, b, a, b and a, b). This results in the damage value being maximum everywhere in this recently broken sea ice area in the CPL_DMG simulation (Fig. a). Comparing damage between the CPL_DMG and REF simulations (Fig. b), we note that the increase in damage related to wave-induced fragmentation is strong at the immediate proximity of the ice edge but lower, although still reasonable ($\simeq$ 0.05), closer to the limit between broken and unbroken ice. Comparing Fig. b with Fig. b, it is interesting to note that the highest difference in the magnitude of ice drift velocities between CPL_DMG and REF is mostly located in this area, where the increase in damage due to waves is limited (Fig. b). This result, the biggest impact on the sea ice drift occurring where the impact on damage is the weakest, is at first counter-intuitive, but it is due to the nature of sea ice at the limit between broken and unbroken ice, which is thicker and more compact than at the proximity of the ice edge. As mentioned before, thin, loose sea ice does not provide much resistance to deformation. For such sea ice, the level of damage has therefore little impact on its behaviour. Conversely, thick compact ice is usually associated with high ice strength values, and its level of damage significantly impacts its resistance to deformation. The additional damage of thick compact ice due to wave-induced fragmentation, despite being small, allows for more sea ice convergence than in the REF and CPL_DMG simulations. Similarly, on 21 October, the difference in ice drift velocity between CPL_DMG and REF is mostly due to an acceleration of compact sea ice that was recently broken (Fig. d). This acceleration follows the wind direction, creating additional convergence north of Svalbard and divergence at the centre of the domain (Fig. c), thus increasing the export of sea ice this time.

To better quantify the impact of this additional damage on the dynamics of compact sea ice, we compare the ice drift velocity between the CPL_DMG and CPL_WRS simulations averaged over the area covered by compact and recently broken sea ice only in Fig. b. This area includes all ice within the domain delimited by the magenta and green solid lines in Fig. b and d. For this particular part of the sea ice cover, the differences in ice drift velocity magnitude are very significant, with an increase of more than 20 % on 16 October, and exceeding 40 % on 21 October. Over the whole 10 d following the fragmentation event, the ice drift velocity for recently broken and compact sea ice increases on average by 7 % in CPL_DMG compared to CPL_WRS.

We also note that for both events the maximum ice drift velocity difference between CPL_DMG and the two other simulations (Fig. a, b) does not happen when the wind speed, ice drift velocity and wave height reach their maxima, but rather a few hours or days later (Figs.a, b and a). A possible explanation is that for strong winds and waves, the magnitude of the sum of the external stresses applied to the ice is high enough to overcome the ice internal stress (in all three simulations). However, for lower wind speeds and wave heights, the magnitude of the external stress decreases. In these conditions, only sea ice with a high enough damage value continues to be deformed. The effect of additional damage is therefore more likely to be maximum in the wake of a storm than during the storm itself.

Waves also impact the sea ice dynamics through the WRS. Figure b and d show the relative importance of the WRS compared to the wind stress as well as the direction in which both apply. On 21 October (Fig. d), the WRS exceeds the wind stress over the first kilometres of the sea ice cover, where the sea ice is rather thin and not compact. As described in , the WRS direction tends to be aligned with the wind at the sea ice edge, but further into the ice cover it aligns with the gradient of ice concentration or thickness. This is because the mean wave direction in the ice cover corresponds to the least attenuated waves, which tend to be the waves that have travelled the shortest distance in-ice. As a consequence, the WRS is most often directed on-ice and is thus a source of sea ice convergence . When the WRS is taken into account, the external stress applied to the sea ice under on-ice (respectively, off-ice) wind conditions is enhanced (respectively, reduced).

Figure 12

Temporal evolution of (a) the ice drift velocity averaged over the entire ice-covered part of the domain (see Fig.  for the domain definition) for the three simulations (REF in green, CPL_WRS in red and CPL_DMG in blue), and (b) the difference of ice drift velocity in the region covered by compact sea ice that has been recently broken (defined as $D_{max}\le \mathrm{200}$ m and $c\ge \mathrm{0.8}$) between the CPL_WRS and CPL_DMG simulations. The time period shown covers 15–25 October 2015, for which initial conditions are given in Fig. . The sea ice-covered part of the domain corresponds to the area for which the sea ice concentration, $c$, is greater than 0. The two orange vertical lines indicates the dates of the snapshots shown in Figs. , and . In the panel legends, $\mathit{\mu }$ indicates the temporal mean associated with each curve and ${\mathit{\sigma }}_{\mathrm{rel}}$ is the standard deviation divided by the mean.

[Figure omitted. See PDF]

With our coupled model, the impact of the WRS on the sea ice dynamics in the MIZ strongly depends on the activation of the link between wave-induced sea ice fragmentation and damage. In particular, these interactions between the WRS and sea ice damage contribute to some of the differences in ice drift velocity between the simulations we see in Fig. a. For the 21 October case, we see a larger difference between CPL_DMG and CPL_WRS than between CPL_DMG and REF. In REF, the off-ice wind stress reduces the sea ice concentration, lowering its effective stiffness, so that sea ice drifts freely with the wind. In CPL_WRS, the WRS cancels some of the wind stress, as well as compacting the sea ice, which allows it to better resist the off-ice external stress. In CPL_DMG, the previous damage to the ice from the waves lowers its stiffness and it again drifts freely in the off-ice direction. The WRS thus limits the rate of sea ice export in off-ice wind events to an extent determined by the rheology of the broken ice. Conversely, on 16 October (Fig. b), the wind stress is directed on-ice and therefore aligned with the WRS. Moreover, the relative importance of the WRS exceeds that of the wind stress over a wide part of the recently broken sea ice area. In these conditions, the WRS contributes to the acceleration of the sea ice convergence, which results in a larger difference of average ice drift velocity over the domain between CPL_DMG and REF than between CPL_DMG and CPL_WRS (Fig. a).

Figure 13

Snapshots of the distributions of sea ice drift velocity in the CPL_DMG simulation (a, c) and of the differences in sea ice drift velocity between the CPL_DMG and REF simulations (b, d) in the Barents Sea taken on 16 October 2015 at 18:00:00 GMT (a, b) and on 21 October 2015 at 09:00:00 GMT (c, d). The green line delimits the area with compact sea ice (defined as $c\ge \mathrm{0.8}$) in CPL_DMG. The magenta line corresponds to the contour $D_{max}=\mathrm{200}$ m in CPL_DMG. The black thick lines delimit the domain used to analyse the results from the three simulations (REF, CPL_WRS, CPL_DMG) in Figs. , and .

[Figure omitted. See PDF]

Figure 14

Snapshots of the distributions of difference in sea ice thickness (a, c) and sea ice concentration (b, d) between the CPL_DMG and the REF simulations (a, b) and between the CPL_WRS and REF simulations (c, d). All snapshots were taken on 21 October 2015 at 09:00:00 GMT. The magenta line corresponds to the contour $D_{max}=\mathrm{200}$ m in CPL_DMG (a, b) and in CPL_WRS (c, d). The black thick lines delimit the domain used to analyse the results from the three simulations (REF, CPL_WRS, CPL_DMG) in Figs. , and .

[Figure omitted. See PDF]

These interactions between sea ice damaging related to wave-induced fragmentation and WRS also have an impact on sea ice properties in the MIZ. Figure a and b display the differences in the sea ice thickness and concentration fields between CPL_DMG and REF. Compaction of sea ice in the MIZ is clearly visible, with thicker, more compact sea ice in the area between the magenta and green solid lines, which corresponds to compact sea ice that has recently been broken. In contrast, between CPL_WRS and REF there are almost no differences in the sea ice thickness fields (Fig. c). Sea ice concentration is only impacted in the first few kilometres of the ice cover, with compaction of the sea ice visible on the western part of the domain (i.e. lower sea ice concentration than in REF at the ice edge, but higher concentration slightly further in the ice cover, Fig. d). From this, we can conclude that when sea ice is damaged by wave-induced fragmentation, it enables sea ice convergence over the entire broken sea ice area. As a consequence, on-ice wind events lead to thicker, more compact sea ice in the MIZ and this phenomenon is enhanced by the WRS. When wave-induced sea ice fragmentation has no impact on sea ice rheology, compacted sea ice, if it is not damaged, resists convergence, preventing sea ice from thickening in the MIZ in CPL_WRS.

Note finally that the relative thickening of sea ice in CPL_DMG increases the wave attenuation, leading to a lower broken sea ice extent in CPL_DMG compared to CPL_WRS (visible in Fig. b, d, for instance). As an example, throughout October in the domain we defined in the Barents Sea, sea ice thickness is on average 2.9 % thicker in CPL_DMG than in CPL_WRS, while the ratio of recently broken sea ice surface area over the total sea ice surface area decreases by 2.4 % between CPL_DMG and CPL_WRS.

5 Discussion

In our case study, the damage added by wave-induced sea ice fragmentation did not significantly enhance sea ice deformation during wave-induced fragmentation events, but after them, when the sea state relaxed. This is because these fragmentation events coincided with high wind speeds; wind stress dominated the internal stress of sea ice in all simulations regardless of the level of damage. Once the wind speed lowers, the internal stress of sea ice dominated the wind stress in places where the sea ice is compact and not damaged, and limited deformation. However, in regions that were previously damaged by wave-induced fragmentation, the level of damage remained high in the first few days following the storm and sea ice could still deform relatively freely. This high level of damage significantly enhanced sea ice mobility in the MIZ in the CPL_DMG simulation compared to the CPL_WRS simulation. This behaviour of the MIZ, with fragmentation events followed by calm periods during which sea ice mobility is enhanced, is not limited to the particular event we describe here. In the Barents Sea, for instance, maxima in the difference between ice drift velocities in the CPL_WRS and CPL_DMG simulations during October 2015 occur after maxima in the ice drift velocity magnitude (Fig. ). We also noted a similar behaviour in the Greenland Sea (not shown).

The impact of wave-induced fragmentation on sea ice dynamics is expected to vary spatially and temporally with waves and sea ice conditions in the Arctic. From our results, the magnitude of this impact depends mostly on the distance over which waves break the ice. At the beginning of October, the low sea ice extent combined with a high frequency of storms favour the occurrence of fragmentation events over a large area, and waves therefore have a strong impact on sea ice dynamics in the MIZ. This impact seems to decrease towards the end of October, as visible in Fig.  where the difference in drift speed between CPL_DMG and CPL_WRS does not show any peaks after the one on 21 October. It coincides with refreezing in the Barents Sea. The generation of thin sea ice over wide areas reduces the fetch and damps the waves over long distances before they can break more compact sea ice. It is therefore likely that, in autumn and winter, fragmentation has little effect on sea ice dynamics as waves are attenuated by the thin and loose forming sea ice. In contrast, in spring and summer, with the sea ice melt increasing the available fetch and exposing compact ice to the waves, we expect wave-induced sea ice fragmentation to have a significant impact on the MIZ. Note also that in melting conditions, sea ice healing does not occur, lengthening the effects of fragmentation over time.

In addition to the effects of interactions between sea ice fragmentation and the WRS on sea ice dynamics, our study allows us to isolate the effect of the WRS in neXtSIM. The effects of the WRS in a coupled wave–sea-ice model system have already been discussed by (coupling WW3 with the sea ice model LIM3) and (coupling neXtSIM and a simplified wave-in-ice model), so we therefore only comment here on the main differences and similarities between these two studies. When the effect of wave-induced sea ice fragmentation on the dynamics is not included in the ice model (in CPL_WRS), the WRS pushes the sea ice edge towards pack ice and increases the sea ice concentration gradient over the first kilometres of the MIZ. This compaction is, however, limited to regions where sea ice poses little resistance to deformation, either because it has a low concentration (hence ice strength) or because it has been previously damaged by the wind. This is similar to the results of the study by in which they note that sea ice drift is only impacted by the WRS in low-concentration sea ice areas. However, also found an acceleration by nearly 10 % of the sea ice drift velocity in areas where sea ice has been broken compared to an uncoupled reference simulation. This is not the case in our study. The acceleration observed by is attributed to the effect of the WRS on sea ice drift in regions with low sea ice concentration where waves can be generated in ice. In our simulations, low-concentration regions allowing for in-ice wave generation are very limited, and this effect is therefore not present. These differences in sea ice concentration distributions might be related to the differences in the sea ice models and external forcings but also to the different time period investigated. It is also interesting to note that reported a decrease of the sea ice melt in their coupled ocean–sea-ice–wave framework, explained by the on-ice drift force associated with the WRS pushing sea ice away from the open water. In our case, we have described how sea ice that has been damaged by wave-induced fragmentation can be exported by an off-ice wind despite a resultant stress reduced by the on-ice push of the WRS. It would therefore be interesting to add an ocean component to our wave–sea-ice coupled framework in order to compare the effects of wave–sea-ice interactions on sea surface properties with the results found by .

In the study by , the WRS was shown to have little effect on the ice edge location, even when sea ice fragmentation was lowering the sea ice internal stress. This difference with our study is likely to be due to the removal of the pressure term used with the EB rheology that gave some resistance to compression, which was especially needed once the ice became damaged . Fragmented ice in our model is thus easier to pile up when convergence occurs. In reality, the collisions and subsequent transfer of momentum between the floes should generate some internal stress that should resist sea ice convergence, lowering the thickening at the ice edge we report here. The acceleration of broken, compact sea ice we give here is therefore likely to be an upper bound of the effect of wave fragmentation on sea ice deformation.

Our results show how waves can modulate the extent of the region of sea ice drifting freely in the MIZ. In most sea ice models, this extent is generally a function of sea ice concentration and thickness, and ignores wave activity. have recently used ICESAT-2 observations to show how the extent of an MIZ defined on a sea ice concentration criterion can differ from an MIZ defined on a wave activity criterion. This is similar to what we show here and our model simply illustrates the potential effects of waves on sea ice dynamics in the MIZ. We show that to impact sea ice dynamics waves need to propagate far enough in the ice cover to fragment compact ice. One of the main uncertainties of our results therefore lies in the estimation of the ice cover area affected by waves. In our model, the extent of this area could be affected by two factors: the relaxation time associated with healing of damage (${\mathit{\tau }}_{\text{heal}}$) and the extent of broken ice.

The relaxation time associated with ${\mathit{\tau }}_{\text{heal}}$ controls the speed at which the ice “forgets” the impact of previous fragmentation events. Its impact was investigated by re-running our experiments using ${\mathit{\tau }}_{\text{heal}}=\mathrm{15}$ d instead of 25 d, the default value in neXtSIM, where 15 d corresponds to the lower limit of the range of ${\mathit{\tau }}_{\text{heal}}$ values for which neXtSIM reproduces well the multi-scaling of sea ice deformation , while 25 d is close to the upper limit of this range. We found that changing the value of ${\mathit{\tau }}_{\text{heal}}$ had very little effect on our results. The sensitivity to ${\mathit{\tau }}_{\text{heal}}$ is possibly low because, once sea ice is broken, waves are able to propagate more easily in the sea ice cover and maintain quite a high level of damage until the wave activity strongly decreases or the sea ice cover extends.

Figure 15

Temporal evolution of (a) the ice drift velocity averaged over all of the ice-covered parts of the domain (see Fig.  for the domain definition) for the CPL_DMG simulation throughout October 2015 and (b) the difference of ice drift velocity in the region covered by compact sea ice that has been recently broken (defined as $D_{max}\le \mathrm{200}$ m and $c\ge \mathrm{0.8}$) between the CPL_WRS and CPL_DMG simulations. The sea ice-covered part of the domain corresponds to the area for which the sea ice concentration $c$ is greater than 0. The two orange vertical lines indicates the dates of the snapshots shown in Figs. , and . In each panel legend, $\mathit{\mu }$ indicates the temporal mean of the plotted quantity.

[Figure omitted. See PDF]

The extent of broken ice depends strongly on the wave attenuation computed in WW3 as discussed in and . The stronger the attenuation, the narrower the extent of broken ice, and the less likely it is that waves can break compact ice and impact sea ice dynamics. We have extended here the parameterization evaluated by to a pan-Arctic simulation that can include very different sea ice conditions. This parameterization, like most wave-in-ice attenuation models, is very sensitive to the values of sea ice concentration and thickness , and the extent of broken sea ice we show here is therefore subject to very large uncertainties. The fact that storm waves can propagate and break the ice over tens to hundreds of kilometres in the Barents Sea is, however, not surprising see the event reported in and consistent with the results of , who found that the Barents Sea is the most consistent region of high wave–ice activity in the Arctic. Comparing our wave attenuation model with other parameterizations used in wave–sea interactions studies is not straightforward, as its attenuation varies strongly with sea ice thickness and floe size. It was shown by to generally yield much more attenuation than the scattering model by from which the attenuation model used by is derived. used the empirical wave attenuation formula by , which is also implemented in WW3. This empirical parameterization was shown by to produce rather strong wave attenuation compared to other empirical formulations. The formula used by does not account for floe size and thickness, but the sensors used in the experiment from which it was derived were located on floes with freeboards of 10 cm or more , and it is therefore calibrated for a total thickness of about 1 m or more. If we compare the formula with our parameterization in the Beaufort Sea case used for the evaluation (not shown), our parameterization leads to a faster decay of the wave height below 50 cm the level at which it stops breaking the ice in when sea ice is thicker than 50 cm (not shown) and slower when sea ice is thinner. In our case study in the Barents Sea, sea ice thickness in the area broken by waves varies between $\simeq$ 25 cm and 1 m, the wave attenuation we estimate should therefore be in line with the one given by the formula.

The estimation of the extent of broken ice is also likely to depend on the model used to determine the occurrence of sea ice break-up. We use a break-up model identical to most studies interested in wave–ice interactions e.g.. It remains extremely simplified and assumes that break-up only occurs in the case of flexural failure in one dimension. However, recent results from laboratory experiments tend to show that there is not such a clear relationship between the wave forcing and the floe size resulting from fragmentation. This is because a complete break-up model should include effects that are currently missing (e.g. from the floe shape and size, two-dimensional flexure modes, floe–floe collisions, rafting). We also point out that while our model includes some memory of previous fragmentation events, we do not account for the fatigue of the ice when determining if break-up occurs or not. The slow-growth FSD is used to keep a memory of the distribution of consolidated floes. It is associated with large-scale mechanical properties of the ice cover, while fatigue is related to the micro-structure of the ice. Accounting for fatigue could significantly lower the ice resistance to flexural failure in some events .

Although our main results are related to the effects of waves on sea ice dynamics and only depend on the FSD's capacity to provide a good estimate for the maximum and average floe size used in the wave model, this work also takes the opportunity to introduce a coupled wave–sea-ice framework using two FSDs to distinguish between two definitions of the floe size: one is more relevant to thermodynamical processes, which grow fast when refreezing occurs, and one more relevant to dynamical and mechanical processes, which are associated with slower growth of the floe size of interest. Theoretically, the processes included allow us to represent the evolution of the FSD year-round even though simulations over different and longer time periods require further evaluation. The introduction of these two FSDs highlights the importance of distinguishing the different timescales involved in the floe size evolution.

As in the studies by and , the main uncertainty related to the FSDs concerns the way sea ice is redistributed after fragmentation. In these early days of the implementation of FSDs in sea ice models, we have built on what was done in wave-in-ice models and used a redistribution scheme that yields FSDs relatively similar to the ones described by , although their methods and interpretations have been contested . Our redistribution scheme could easily be adapted if progress is made on the relative importance of the parameters affecting this cut-off floe size and the associated redistribution, either from observations thanks to new available FSD datasets or from discrete element modelling .

6 Conclusions

Using a coupled wave–sea-ice model, we have shown how waves may contribute to modifying the sea ice dynamics in the MIZ by modulating the extent of ice regions that pose little resistance to deformations. As noted by , this extent does not necessarily coincide with areas of low ice concentrations. With our model, we note a significant acceleration of compact sea ice in both convergent and divergent sea ice drift conditions once sea ice has been fragmented by waves. Even though some assumptions we make here require further evaluation, the results are of particular interest as they highlight missing physics in current modelling systems used for short- and long-term sea ice predictions, and concern key areas of the polar regions.

Reliable sea ice forecasts are essential to ensure the safety of human activities close to the MIZ. In this context, waves pose a hazard as they make sea ice more mobile. Our results therefore stress the need for the addition of wave effects in sea ice models used in forecasts. On longer timescales, the impact of waves on sea ice dynamics could affect the amount of sea ice that is exported to the ocean. Indeed, eddies and/or filaments are likely to play an important role in this export in cases where sea ice dispersion is possible . Moreover, the fragmentation of sea ice itself could also generate sub-mesoscale and mesoscale activity in the ocean . Future coupling with an ocean model could therefore bring new insight into the interactions between waves, sea ice and the ocean in the MIZ.

Appendix A Table A1

List of symbols used in this study and their associated quantities. Parameter default values, or references used to compute each quantity, are indicated in the last column.

Symbol	Quantity	Value/reference
$c$	Sea ice concentration	–
$h$	Sea ice thickness	–
$D$	Floe size (mean calliper diameter)	–
$D_{max}$	Maximum floe size	Initialized with 1000 m
$Y$	Young modulus	5.49 GPa
$\mathit{\nu }$	Poisson's ratio	0.3
$G_{\text{r}}$	Lateral melt rate	See
${\dot{c}}_{\text{new}}$	Rate of formation of new ice	See
${\mathit{\beta }}_{\text{weld}}$	FSD redistribution term associated with welding of floes	See
$\mathit{\kappa }$	Rate at which the number of floes decreases due to welding per surface area	$\mathrm{5}\times {\mathrm{10}}^{-\mathrm{8}}$ m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$
${\mathit{\tau }}_{\text{heal}}$	Relaxation time associated with damage (mechanical) healing	25 d (default)
$\mathrm{\Delta }{t}_{\text{ice}}$	Ice model time step	20 s
$\mathrm{\Delta }{t}_{\text{cpl}}$	Coupling time step	2400 s
${\mathit{\tau }}_{\text{WF}}$	Relaxation time associated with fragmentation	1800 s
$Q$	Redistribution probability function of floe size associated with fragmentation	See
$\mathit{\beta }$	Redistribution factor of floe size associated with fragmentation	See
${\mathit{\lambda }}_{\text{break}}$	Shortest wavelength triggering flexural failure of sea ice	See , Sect. 2.3
$D_{\text{FS}}$	Minimum floe size for flexural failure	See
$p_{\text{FS}}$	Probability that ice breaks depending on $D/D_{\text{FS}}$	Eq. ()
$p_{\mathit{\lambda }}$	Probability that ice breaks depending on $D/{\mathit{\lambda }}_{\text{break}}$	Eq. ()
$c_{\mathrm{1},\text{FS}}$	Value of $D/D_{\text{FS}}$ under which ice cannot be broken	1.
$c_{\mathrm{2},\text{FS}}$	Parameter controlling the range of $D/D_{\text{FS}}$ over which $p_{\text{FS}}$ goes from 0 to 1	2.
$c_{\mathrm{1},\mathit{\lambda }}$	Value of $D/{\mathit{\lambda }}_{\text{break}}$ under which ice cannot be broken	0.3
$c_{\mathrm{2},\mathit{\lambda }}$	Param. controlling the range of $D/{\mathit{\lambda }}_{\text{break}}$ over which $p_{\text{FS}}$ goes from 0 to 1	2.
$c_{\text{broken}}$	Concentration of broken sea ice	Eq. ()
$d_w$	Damage value associated with broken ice	0.99
$\mathit{\alpha }$	Exponent of the small-floe regime power-law FSD	See ; note
		that the sign is reversed
$q$	Exponent used in the redistribution factor	Eq. ()

Code and data availability

The wave–sea-ice coupling routines developed in WW3 will be included in the future release of WW3 and made publicly available at https://github.com/NOAA-EMC/WW3 (last access: November 2020). Model outputs are available upon request from [email protected].

Author contributions

GB and TW developed the theory and set up the coupling with the support of EO. GB led the analysis and wrote the manuscript in close collaboration with TW and with contributions from all co-authors.

Competing interests

The authors declare no competing interests.

Acknowledgements

This work was carried out as part of the Copernicus Marine Environment Monitoring Service (CMEMS) WIzARd project. CMEMS is implemented by Mercator Ocean in the framework of a delegation agreement with the European Union. It was also funded by institutional support from the Nansen Center. We thank the editor, the two anonymous reviewers and Christopher Horvat for their suggestions and their support for the publication of the manuscript. Finally, we thank Fabrice Ardhuin for providing the wave height data from SAR.

Financial support

This research was supported by the CMEMS (WIzARd) and the NERSC (internal funding).

Review statement

This paper was edited by Daniel Feltham and reviewed by Christopher Horvat and two anonymous referees.