Wave–sea-ice interactions in a brittle rheological framework

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 (Meier, 2017) 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 (Thomson and Rogers, 2014) 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. Langhorne et al., 1998; Collins et al., 2015). 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 (Tietsche et al., 2014). On shorter timescales, large uncertainties remain in the forecasts of the position of the sea ice edge (Schweiger and Zhang, 2015; DeSilva and Yamaguchi, 2019); however, this information is essential for the safety of the increasing number of human activities in polar regions (Yumashev et al., 2017). 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, Longuet-Higgins, 1977), 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, Stopa et al. (2018b) found it to be as important as wind stress over the first 50 km of the MIZ in the Southern Ocean, whereas Alberello et al. (2020) 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. McPhee, 1980), having lower internal stress. This seems to be consistent with the deformation observations of Oikkonen et al. (2017) 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. Shen et al. (1986) 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. Feltham (2005) 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 Johannessen et al. (1983). 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 (Horvat et al., 2019). Vichi et al. (2019) 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 (Thomson et al., 2018; Squire, 2018). 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 (Wadhams et al., 1986; Montiel et al., 2016), 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 Squire et al. (1995) and Squire (2007, 2020) gave a more detailed history of this area of work. Liu and Mollo-Christensen (1988) and Collins et al. (2015) 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 (Dumont et al., 2011; Williams et al., 2013a, b). 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_max) depending on the wave field (Dumont et al., 2011). This assumption about the shape of the FSD is based on observations made by Toyota et al. (2011) of power-law FSD in the MIZ that they explain by successive fragmentation of floes by waves. More recently, Boutin et al. (2018) included a parameterization in the spectral wave model WAVEWATCH III (WW3: The WAVEWATCH III Development Group, 2019) 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 Collins et al. (2015). Ardhuin et al. (2018) 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 (Williams et al., 2017; Bennetts et al., 2017; Bateson et al., 2020; Roach et al., 2018) and more recently using more complex spectral wave models like WW3 (Boutin et al., 2020; Roach et al., 2019). 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 (Zhang et al., 2015; Horvat and Tziperman, 2015). 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 (Asplin et al., 2012).

The dynamical aspects of wave–ice interactions have received less attention. Boutin et al. (2020) 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, Rynders (2017) 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 (Hunke, 2010). Williams et al. (2017) 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, Rampal et al., 2019) in a stand-alone set-up. The elasto-brittle (EB) rheology (Girard et al., 2011) 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 (Bouillon and Rampal, 2015; Rampal et al., 2016). 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 Williams et al. (2017) 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 Williams et al. (2017) 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 Dansereau et al., 2016), greatly improving sea ice deformation in pack ice (Rampal et al., 2019). 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 (Ardhuin et al., 2016; Boutin et al., 2018; Ardhuin et al., 2018) and extends the work done by Williams et al. (2017) 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 Zhang et al. (2015) and Horvat and Tziperman (2015), and of the coupling of WW3 with the sea ice model LIM3 described by Boutin et al. (2020). 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^® (The WAVEWATCH III Development Group, 2019), building on the previous developments reported by Boutin et al. (2018), 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 Collins et al. (2015) and waves propagating in forming ice in the Beaufort Sea (Ardhuin et al., 2018). As in the study by Ardhuin et al. (2018), 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 (Sutherland and Rabault, 2016).

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

The sea ice model we use for this coupling is neXtSIM (Rampal et al., 2019), in which an FSD is first implemented as described in Sect. 2.2. The two models are coupled through the coupler OASIS-MCT (Craig et al., 2017). Figure 1 shows the variables that are exchanged. Briefly, WW3 determines if the waves will break the ice and calculates a representative wavelength λ_break in the manner of Boutin et al. (2018). This is then used by neXtSIM to modify the FSD as described below in Sect. 2.2.2. 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 A1.

2.1 Wave radiation stress on sea ice

As in Boutin et al. (2020), 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 Williams et al. (2017). As discussed by Williams et al. (2017), 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. Roach et al., 2018, 2019). 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 x is the position vector, t is time and D_fast is defined as the mean calliper diameter of the floes separated from each other by leads or cracks, as introduced by Rothrock and Thorndike (1984).

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 Rampal et al. (2016) with values of ≃ 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_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.

The picture in Fig. 2 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. 2 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 (Rampal et al., 2008; Marsan and Weiss, 2010). An ice cover such as the one illustrated in Fig. 2 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 (Liu and Mollo-Christensen, 1988; Collins et al., 2015; Kohout et al., 2016), but our intuition is supported by the recent Antarctic observations made by Kohout et al. (2016), 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. Zhang et al., 2015; Roach et al., 2018; Bateson et al., 2020; Boutin et al., 2020), which are defined as

$$\begin{array}{}\text{(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)\end{array}$$

and

$$\begin{array}{}\text{(2)}& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}g(D,\mathit{x},t)\text{d}D=\mathrm{1}.\end{array}$$

In these definitions, g(D) represents an FSD where D is its associated floe size, A is the total area considered around the position 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+dD. The value D=0 corresponds to open water and Eq. (2) 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 ΔD with edges between a minimum and a maximum floe size, respectively, D₀ and D_N. Thus, floes with sizes in [D₀ D₁] cannot be broken into smaller pieces, and we refer to floes with sizes in [D_N−1 D_N] as unbroken floes. Using fixed-width bins may bias our ability to represent or examine scale-invariant behaviour (Stern et al., 2018), 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 Roach et al. (2018) or Boutin et al. (2020):

$$\begin{array}{}\text{(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}},\end{array}$$

in which u corresponds to the sea ice velocity vector, Φ_th is a redistribution function of floe size due to thermodynamic processes (i.e. lateral growth/melt) and Φ_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 Φ_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 Φ_th and Φ_m differ between g_fast and g_slow and are described below.

2.2.1 Lateral sea ice melt/growth

In this section, we describe the implementation of the terms Φ_th,slow and Φ_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 Roach et al. (2018):

$$\begin{array}{}\text{(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}\end{array}$$

where G_r is the lateral melt rate of floes, ${\dot{c}}_{\text{new}}$ is the rate of formation of new ice and β_weld is the FSD redistribution term associated with welding of floes using the Smoluchowski equation as implemented by Roach et al. (2018).

Lateral melting is implemented following Horvat and Tziperman (2015) and Roach et al. (2018). Here, we neglect lateral melt for the largest floe size category as floes with size O(100) 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 ϕ_lead in Roach et al. (2018) 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 Roach et al. (2018), 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 Roach et al. (2018) 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 Roach et al. (2018), we decided to carry on with this formulation but with an appropriate tuning of the coefficient that Roach et al. (2018) called κ, which represents the rate at which the number of floes decreases due to welding per surface area. We tune κ 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 Roach et al. (2018). To give an idea of the timescales involved, the model of Roach et al. (2018), 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=0.95). In the simulations presented in this study, setting $\mathit{\kappa }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{8}}$ m⁻² s⁻¹ 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 Φ_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 Φ_th,fast. It is very similar to the “damage healing” already included in neXtSIM (see Rampal et al., 2016), and associated with a timescale τ_heal, that we reuse in the computation of Φ_th,slow:

$$\begin{array}{}\text{(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}},\end{array}$$

where Φ_conservation is an ad hoc term ensuring the conservation of sea ice area (i.e. Eq. 2) 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 τ_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, Φ_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 Φ_m,slow and Φ_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 Boutin et al. (2020), 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 λ_break, is passed to neXtSIM (see Fig. 1) and used in the mechanical redistribution scheme of the FSD. The determination of the value of λ_break in WW3 was explained in detail in Sect. 2.3 of Boutin et al. (2018) (where it is called λ_i,break). If no fragmentation has occurred in WW3, neXtSIM receives λ_break=1000 m, which corresponds to the default unbroken value, and no fragmentation occurs in neXtSIM (resulting in Φ_m,slow=0 and Φ_m,fast=0). If neXtSIM receives a value of λ_break<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 Kohout et al. (2016). 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_fast, to the slow-growth FSD, g_slow. In practice, we set ${\mathrm{\Phi }}_{\text{m,thermo}}\mathrm{\Delta }{t}_{\text{ice}}=g_{\text{fast}}-g_{\text{slow}}$, where Δt_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 (Collins et al., 2015) and (ii) the fast-growth FSD g_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_fast=g_slow. This equality represents the fact that D_fast=D_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 Φ_m,slow. This term can be written in the same form as in Zhang et al. (2015):

$$\begin{array}{}\text{(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}\end{array}$$

where Q(D) is a redistribution probability function characterizing which proportion of floes of a given size D is broken (with D representing indistinctively D_fast and D_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^′ to D as fragmentation occurs. The choices of Q(D) 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 (Boutin et al., 2018) and lateral melting (Bateson et al., 2020), but the evolution of the FSD during wave-induced fragmentation is still not well understood.

Toyota et al. (2011) 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 Herman et al. (2018) 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 (Burroughs and Tebbens, 2001; Stern et al., 2018). The division of the FSD into two truncated power laws suggested by Toyota et al. (2011) has nevertheless been used to redistribute FSDs in the wave-in-ice models by Dumont et al. (2011) and Williams et al. (2013b), which have been reused for the wave-in-ice attenuation parameterization implemented in WW3 by Boutin et al. (2018) and in many other studies using wave-in-ice models (see, e.g. Aksenov et al., 2017; Williams et al., 2017; Bennetts et al., 2017; Bateson et al., 2020; Boutin et al., 2020).

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 Boutin et al. (2018) and Ardhuin et al. (2018). We therefore build on the work by Zhang et al. (2015) and Boutin et al. (2020) to suggest a new parameterization for Q(D) 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 Williams et al. (2013b). However, there are two main differences from the FSDs assumed in Williams et al. (2013b): 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 Toyota et al. (2011).

Q(D) represents the amount of ice in each floe size category that will be broken. Williams et al. (2013b) and a number of further studies using their approach (Bennetts et al., 2017; Bateson et al., 2020; Boutin et al., 2020) used a step function with Q(D)=1 if D≥0.5λ_break and D≥D_FS, where D_FS is the minimum floe size for which flexural failure can occur (see Mellor, 1986), and Q(D)=0 if not. The transition of Q(D) from 0 to 1 occurs at a critical floe size under which floes do not break, corresponding to what Toyota et al. (2011) interpreted as the cut-off floe size. In the model of Williams et al. (2013b), the transition therefore occurs at $max(D_{\text{FS}},\mathrm{0.5}{\mathit{\lambda }}_{\text{break}})$. D_FS depends on sea ice properties and is computed as

$$\begin{array}{}\text{(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}},\end{array}$$

where g is gravity, Y is the Young modulus of sea ice, ν is Poisson's ratio, h is the mean sea ice thickness and ρ is the density of sea water. For ice thinner than 1 m, which is often the case in the MIZ, D_FS is lower than 15 m and the cut-off floe size is likely to be determined by the value of λ_break. To define Q(D), we take an approach similar to Williams et al. (2013b) and set

$$\begin{array}{}\text{(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}}),\end{array}$$

in which τ_WF is a relaxation time associated with wave-induced fragmentation events, and p_FS and p_λ are probabilities that floes break depending on their size. We introduced τ_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 Collins et al. (2015). The probability functions p_FS and p_λ express the idea that the smaller the floes are, the less chance they have to break. The function p_λ compares the floe size D with the value of D_FS and p_λ compares the floe size D with λ_break, introducing a dependency of Q(D) on the wave field. The difference with the model of Williams et al. (2013b) is that instead of step functions we use hyperbolic tangents to get a continuous transition of Q(D) 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_1,FS, c_2,FS, c_1,λ and c_2,λ are parameters of the FSD that control the range of floe size that will be broken or not. The use of a continuous Q(D) instead of a step function aims to relax the constraint on the FSD shape imposed by Williams et al. (2013b). 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 Toyota et al. (2011) 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 Toyota et al. (2011) could actually be due to windowing issues (Stern et al., 2018). 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_2,FS and c_2,λ. Q(D) tends towards a step function when c_2,FS and c_2,λ 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 Toyota et al. (2011) or in the lab experiments by Herman et al. (2018), for instance. The location of the floe size range at which the transition occurs is controlled by c_1,FS and c_1,λ. Like Williams et al. (2013b), we set $c_{\mathrm{1},\text{FS}}=\mathrm{1}$. Instead of using $c_{\mathrm{1},\mathit{\lambda }}=\mathrm{0.5}$ as in Williams et al. (2013b), we set $c_{\mathrm{1},\mathit{\lambda }}=\mathrm{0.3}$. It is consistent with the hypothesis made by Boutin et al. (2018) that floes smaller than 0.3λ_break are tilted by waves but do not bend and therefore have no chance of suffering flexural failure. Reducing the value of c_1,λ reduces the value of the cut-off floe size we should get, which generally leads to floes smaller than assumed in Williams et al. (2013b). It is, however, compensated for by using a continuous Q(D), which gives more weight to large floes than the FSD assumed in the model of Williams et al. (2013b). FSDs generated by this redistribution function are presented and discussed in Sect. 4.1.2.

The redistribution factor β controls the shape of the FSD after it has been redistributed. Similar to the work by Zhang et al. (2015) and Boutin et al. (2020), the redistribution factor β follows the form

$$\begin{array}{}\text{(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.,\end{array}$$

where n corresponds to the index of the nth FSD category and q is an exponent that controls the shape of the redistributed FSD. Zhang et al. (2015) used q=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. Boutin et al. (2020) noted that q is related to the (negative) exponent α of the power-law FSD associated with the smaller floe sizes resulting from the redistribution as assumed by Toyota et al. (2011) and later by Williams et al. (2013b) with $q=\mathrm{2}+\mathit{\alpha }$. In their study, Toyota et al. (2011) related the exponent α to a quantity they call fragility, which represents the probability of fragmentation of floes. Williams et al. (2013b) assumed that fragility is constant and equal to 0.9, giving $\mathit{\alpha }\simeq -\mathrm{1.85}$. Boutin et al. (2020) set the value q to prescribe power-law FSDs with this same value of α 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 Toyota et al. (2011). Here, we have already expressed the probability of fragmentation of floes, hence the fragility, as the product of p_FS and p_λ. Using this relationship, we get

$$\begin{array}{}\text{(11)}& q=-{\log }_{\mathrm{2}}\left(p_{\text{FS}}\left(D\right)p_{\mathit{\lambda }}(D,{\mathit{\lambda }}_{\text{break}})\right).\end{array}$$

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 α of this power law. Instead, α decreases as the ice thickens and wavelength increases as is shown and discussed in Sect. 4.1.2.

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 (Boutin et al., 2018). When WW3 is run in stand-alone mode, D_max is taken to be λ_break∕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

$$\begin{array}{}\text{(12)}& \langle D\rangle =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}Dg\left(D\right)\text{d}D.\end{array}$$

Here we use the slow-growth FSD, g_slow (D_slow), 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 Dumont et al. (2011) 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 90th percentile of the areal FSD. Besides, the flexure dissipation mechanisms included in WW3 by Boutin et al. (2018) require us to discriminate between sea ice cover made of large floes with size on the order of O(100) 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(100)m, they are only fully attenuated by inelastic dissipation if floe size is on the order of O(1000) m, which can be larger than the floe size range covered by the FSD defined in neXtSIM. In the case where D_N<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 λ_break<1000 m as $c_{\mathrm{broken}}=c\left(\mathrm{1}-{\mathrm{e}}^{-\mathrm{\Delta }{t}_{\text{cpl}}/{\mathit{\tau }}_{\text{WF}}}\right)$ and c_broken=0 if λ_break≥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=0.99 to c_broken, which gives the following evolution of the damage d:

$$\begin{array}{}\text{(13)}& d=min\left(\mathrm{1},d(\mathrm{1}-c_{\mathrm{broken}})+d_w{c}_{\mathrm{broken}}\right).\end{array}$$

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 Boutin et al. (2020), the coupled framework is run on a regional 0.25^∘ 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 (Saha et al., 2014). In addition, neXtSIM is also forced by ocean fields from the TOPAZ4 re-analysis (Sakov et al., 2012). For more details on the forcings used by neXtSIM, see Rampal et al. (2016). Wave–current interactions in WW3 are not considered in this study.

For the FSD in neXtSIM, we use 20 categories with a width ΔD=10 m. The lower bound D₀ is set to 10 m, which is about the size of the smallest floes susceptible to undergo flexural failure (Mellor, 1986). 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 τ_heal used by neXtSIM is set to its default value of 25 d (Rampal et al., 2016).

3.2 Evaluation of the FSD implementation: model set-up

In Sect. 4.1, 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 (Thomson et al., 2018) to evaluate the wave attenuation and broken sea ice extent in our coupled simulations, which is similar to what Ardhuin et al. (2018) 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 Ardhuin et al. (2018). This storm generated ≃ 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 Ardhuin et al. (2018). 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 Ardhuin et al. (2018). We also used the same sea ice concentration data as Ardhuin et al. (2018) 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 (Kaleschke et al., 2001; Spreen et al., 2008)¹. 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 Ardhuin et al. (2018), 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 Ardhuin et al. (2018) 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 ΔD=5 m and D₀=5 m. We also evaluated the evolution of the two FSDs with refreezing/healing using the CPL_DMG simulation described in Sect. 3.3.

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

In Sect. 4.2, 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. 2 but uses the relationship between wave-induced sea ice fragmentation and damage presented in Sect. 2.3. 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. 2.3. 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 (Sakov et al., 2012) 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 Ardhuin et al., 2018). 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 (see Stopa et al., 2018a; Ardhuin et al., 2018, for 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. 4. 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.

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 Boutin et al. (2018) in the study by Ardhuin et al. (2018). 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. 3). We also compared our model results with estimated wave height from SAR images (Stopa et al., 2018a) and buoy measurements (AWAC, see Thomson et al., 2018) along a transect in Fig. 3d.

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 Ardhuin et al. (2018) (see Fig. 9 of their study). Figure 3d 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 Ardhuin et al. (2018), 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. 3d 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. 3a, 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, Nose et al., 2020). Moreover, as in Ardhuin et al. (2018), 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. Toyota et al., 2011; Herman et al., 2018) for three different locations (Fig. 4) 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. Toyota et al. (2011), with the curve gradually steepening as floe size increases. Following the method of Toyota et al. (2011), 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 (Stern et al., 2018). 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 Williams et al. (2013b). Here we use the CDFs to discuss the differences introduced by our parameterization of the FSD redistribution compared to the assumptions made by Williams et al. (2013b) 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^∘ W), and their position is indicated by a symbol of the same colour in Fig. 3a 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_max values for each case.

As in the study by Toyota et al. (2011), 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 −2, as expected for a CDF that depends on a two-dimensional fragmentation process (Toyota et al., 2011). 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₀ for the small-floe size regime to be resolved. In the model of Williams et al. (2013b), 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 Williams et al. (2013b) 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 Williams et al. (2013b) corresponds to the “cut-off” floe size as interpreted by Toyota et al. (2011). 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 Williams et al. (2013b). D_max shows little sensitivity to the FSD definition, with a maximum difference between the two simulations presented in Fig. 4 not exceeding the value of ΔD in the refined FSD simulation (5 m). Considering the large uncertainties due to the limited knowledge of wave–ice interactions, choosing ΔD=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. 6 and the floe size parameters in Fig. 7 are taken.

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^∘ E and 83.5^∘ 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 α are given in the legend. The vertical dotted lines represents the D_max values for each date.

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 5 and 6 illustrate the evolution of the FSDs in the CPL_DMG simulation, in which both mechanical and thermodynamical processes are active. Figure 5a and c show the proportion of “unbroken” ice (the proportion of sea ice associated with the Nth 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. 5e, 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, κ. This is because, like Roach et al. (2018), we use a constant value for κ, 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 τ_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. 5d).

This difference in timescales is also visible in Fig. 6, 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. 5a 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. 5). 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^∘ E and 83.5^∘ N (see red cross in Fig. 5a, c, e).

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. 5). These terms might affect the moments of the distribution in a non-trivial way, as was shown in Sect. 3.1 of Horvat and Tziperman (2017). 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. 7 the evolution of D_max and 〈D〉 for the location indicated in Fig. 5a, 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 τ_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. 9, 12 and 15.

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^∘ N parallels, respectively, and west and east by the 16 and 60^∘ E meridians (see, for instance, Fig. 8).

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. 8 for the initial sea ice conditions), which results in a growth of the surface occupied by broken ice (Fig. 9c). 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. 8 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≥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. 8. 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. 10, 11 and 13.

The fragmentation event occurrs during on-ice wind conditions (see the positive meridional component in Fig. 9a) with high waves (up to 5 m at the ice edge) propagating far into the ice cover (see Figs. 9b and 10b). 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. 8a and 10a) until it represents nearly 60 % of the sea ice-covered part of the domain (Fig. 9c). Here we define the sea ice-covered part of the domain as the area for which sea ice concentration c>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>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. 9c). 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. 9a). 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. 9b). 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. 9c). 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. 10a and 11a) and 21 October (Figs. 10c and 11c), 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. 8b), 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. 9, 12 and 15.

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. 12) 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. 9, 12 and 15.

On 16 October, the wind and waves are directed on-ice, thus compacting the sea ice (Fig. 10b). At the time of the snapshot, sea ice has been recently broken over a wide surface area (within the magenta contour in Figs. 10a, b, 11a, b and 13a, b). This results in the damage value being maximum everywhere in this recently broken sea ice area in the CPL_DMG simulation (Fig. 11a). Comparing damage between the CPL_DMG and REF simulations (Fig. 11b), 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 (≃ 0.05), closer to the limit between broken and unbroken ice. Comparing Fig. 11b with Fig. 13b, 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. 13b). 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. 13d). This acceleration follows the wind direction, creating additional convergence north of Svalbard and divergence at the centre of the domain (Fig. 13c), 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. 12b. This area includes all ice within the domain delimited by the magenta and green solid lines in Fig. 13b 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. 12a, 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.9a, b and 12a). 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 10b 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. 10d), 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 Boutin et al. (2020), 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 (Stopa et al., 2018b; Sutherland and Dumont, 2018). 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. 8 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≥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. 8. 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. 10, 11 and 13. In the panel legends, μ indicates the temporal mean associated with each curve and σ_rel is the standard deviation divided by the mean.

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. 12a. 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. 10b), 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. 12a).

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≥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. 9, 12 and 15.

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. 9, 12 and 15.

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 14a 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. 14c). 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. 14d). 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. 14b, 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. 15). 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. 15 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 Boutin et al. (2020) (coupling WW3 with the sea ice model LIM3) and Williams et al. (2017) (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 Boutin et al. (2020) in which they note that sea ice drift is only impacted by the WRS in low-concentration sea ice areas. However, Boutin et al. (2020) 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 Boutin et al. (2020) 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 Boutin et al. (2020) 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 Boutin et al. (2020).

In the study by Williams et al. (2017), 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 (Rampal et al., 2016). 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 (Shen et al., 1986) 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. Horvat et al. (2020) 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 (τ_heal) and the extent of broken ice.

The relaxation time associated with τ_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 τ_heal=15 d instead of 25 d, the default value in neXtSIM, where 15 d corresponds to the lower limit of the range of τ_heal values for which neXtSIM reproduces well the multi-scaling of sea ice deformation (Rampal et al., 2016), while 25 d is close to the upper limit of this range. We found that changing the value of τ_heal had very little effect on our results. The sensitivity to τ_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. 8 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≥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. 10, 11 and 13. In each panel legend, μ indicates the temporal mean of the plotted quantity.

The extent of broken ice depends strongly on the wave attenuation computed in WW3 as discussed in Boutin et al. (2018) and Ardhuin et al. (2018). 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 Ardhuin et al. (2018) 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 (Doble and Bidlot, 2013; Nose et al., 2020), 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 Collins et al., 2015) and consistent with the results of Horvat et al. (2020), 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 Boutin et al. (2018) to generally yield much more attenuation than the scattering model by Kohout and Meylan (2008) from which the attenuation model used by Roach et al. (2018) is derived. Bennetts et al. (2017) used the empirical wave attenuation formula by Meylan et al. (2014), which is also implemented in WW3. This empirical parameterization was shown by Collins and Rogers (2017) to produce rather strong wave attenuation compared to other empirical formulations. The formula used by Meylan et al. (2014) 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 (Kohout et al., 2016), and it is therefore calibrated for a total thickness of about 1 m or more. If we compare the Meylan et al. (2014) 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 Ardhuin et al., 2018) 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 ≃ 25 cm and 1 m, the wave attenuation we estimate should therefore be in line with the one given by the Meylan et al. (2014) 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. Williams et al., 2017; Bateson et al., 2020; Roach et al., 2019). 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 (Herman et al., 2018; Dolatshah et al., 2018) 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 (Langhorne et al., 1998).

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 Zhang et al. (2015) and Boutin et al. (2020), 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 Toyota et al. (2011), although their methods and interpretations have been contested (Stern et al., 2018; Horvat et al., 2019). 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 (Hwang et al., 2017; Horvat et al., 2019) or from discrete element modelling (Herman, 2018).

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 Horvat et al. (2020), 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 (Manucharyan and Thompson, 2017). Moreover, the fragmentation of sea ice itself could also generate sub-mesoscale (Horvat et al., 2016) and mesoscale activity in the ocean (Dai et al., 2019). 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

Maykut and Perovich (1987)Rampal et al. (2016)Roach et al. (2018)Zhang et al. (2015)Zhang et al. (2015)Boutin et al. (2018)Mellor (1986)Toyota et al. (2011)

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.