Toward a coupled model to investigate wave-sea ice interactions in the Arctic marginal ice zone

. The Arctic Marginal Ice Zone (MIZ), where strong interactions between sea ice, ocean and atmosphere are taking place, is expanding as the result of the on-going sea ice retreat. Yet, state-of-art models are not capturing the complexity of the varied processes occurring in the MIZ, and in particular the processes involved in the ocean-sea ice interactions. In the present study, a coupled sea ice - wave model is developed, in order to improve our understanding and model representation of those interactions. The coupling allows us to account for the wave radiative stress resulting from the wave attenuation by 5 sea ice, and the sea ice lateral melt resulting from the wave-induced sea ice break-up. We found that, locally in the MIZ, the waves can affect the sea ice drift and melt, resulting in signiﬁcant changes in sea ice concentration and thickness as well as sea surface temperature and salinity. Our results highlight the need to include the wave-sea ice processes in models aiming at forecasting sea ice conditions on short time scale, although the coupling between waves and sea ice would probably required to be investigated in a more complex system, allowing for interactions with the ocean and the atmosphere. 10


Introduction
Numerical models exhibit large biases in their representation of the Arctic sea ice concentration and thickness, regardless of their complexity or resolution (Stroeve et al., 2014;Chevallier et al., 2017;Wang et al., 2016;Lique et al., 2016). Comparing 10 reanalyses based on state-of-the-art ocean-sea ice models against observations, Uotila et al. (2018) found that the model biases are the largest in the Marginal Ice Zone (MIZ). Indeed, the MIZ is characterized by a wide variety of processes resulting 15 from the highly non linear interactions between the atmosphere, the ocean and the sea ice (Lee et al., 2012), and many of these processes are only crudely (if at all) taken into account in models. Some of these processes result from the interactions between surface wave and sea ice, and are thought to be key for the dynamics and evolution of the MIZ . These interactions are the focus of the current paper.
Summer sea ice in the Arctic has been drastically receding over the past few decades (Comiso et al., 2017), resulting in an 20 expansion of the MIZ which is expected to be intensified in the future (Aksenov et al., 2017). This offers an expanding fetch for waves to grow and propagate (Thomson and Rogers, 2014), suggesting an overall increase of wave heights in the Arctic (Stopa et al., 2016b). Once generated, waves can then propagate into sea ice, impacting strongly the dynamical and thermody-sea ice edge. The implementation of FSDs in different sea ice models as done by Zhang et al. (2015) or Horvat and Tziperman (2015) has also allowed the assessment of the potential enhancement of lateral melt by wave-induced ice fragmentation (Zhang et al., 2016;Bennetts et al., 2017), but the representation of waves remains too crude to simulate the full effect of waves on the evolution of sea ice. 5 In the present study, we go beyond the simple inclusion of the forcing of wave by sea ice or of sea ice by wave in models, by proposing a full coupling between a spectral wave model and a state-of-art sea ice model. The coupled framework allows us to investigate the interactions between waves and sea ice in the Arctic, and the impact that including these effects in a model has on the representation of the wave, ocean, and sea ice properties in the Arctic MIZ. The remainder of this paper is organized as follows. The different models and configurations used in this study are described in Section 2. Section 3 is devoted 10 to the theoretical and practical implementation of the coupling between the two models. In section 4, we compare pan Arctic simulations in which the coupling between wave and sea ice is implemented or not, in order to quantify the dynamical and thermodynamical impacts on the coupling on the sea ice and ocean surface properties. A summary and conclusions are given in Section 5.

15
In this study we make use of the spectral wave model WAVEWATCH III ® (hereafter WW3; The WAVEWATCH III ® Development Group, 2016), building on the previous developments performed by Boutin et al. (2018) who included a FSD in WW3 as well as a representation of the different processes by which sea ice can affect the propagation and modulation of waves in the MIZ. We also use the sea ice model LIM3 (Vancoppenolle et al., 2009;Rousset et al., 2015), in which a FSD is first implemented as described in Section 3.2. The two models are coupled through the coupler OASIS-MCT (Craig et al., 2017). 20 Two configurations of different complexities are used in the following and briefly described in the remaining of this section.

Idealized configuration
In order to test and illustrate the effect of the coupling (Section 3), we make use of a simple idealized configuration (see Fig. 1), in which LIM3 is used in a stand alone mode (without any ocean component). The configuration is a squared domain with 100 × 100 grid cell, with a resolution of 0.03 o in both directions (corresponding roughly to 3 km). Both models are run 25 on the same grid, and with the same time step set to 300s. The coupling time step is set to 300s too. The sea ice is only forced by waves, without any wind or ocean current. Following Boutin et al. (2018), the simulation starts at rest, with distributions of sea ice concentration (Fig. 1a) and thickness (Fig. 1c) set to represent roughly the conditions that can be found in the MIZ.
Starting from the western border, the domain is free of ice over the first 10 km, after which the ice concentration c increases linearly from 0.4 to 1 about 90 km further eastward (longitude=0.84 o E). Ice thickness also increases from west to east follow- 30 ing h i = 2(0.1 + e −Nx/20 ), where N x is the number of grid cells in the x direction starting from the western border of the ice covered domain. Waves radiate from part of the western border of the domain, between 1.2 and 1.8 o of latitude, and propagate to the east. The wave spectrum at the boundary is extracted from an Arctic realistic simulation using WW3 described by Stopa et al. (2016a) for May 2nd and 3rd, 2010, during which a storm happened south of Svalbard (Collins et al., 2015). Here we rotate the spectrum so that the direction with the largest density of wave energy is lined up with our x-axis. Simulations start on April 30th, at 02:00 a.m., and the attenuation processes (scattering, bottom friction and inelastic dissipation) use the same parameterization as in the reference simulation described by Ardhuin et al. (2018).

Pan-Arctic configuration
We also make use of the CREG025 configuration (Dupont et al., 2015;Lemieux et al., 2018), which is a regional extraction of the global ORCA025 configuration developed by the Drakkar consortium (Barnier et al., 2006). Although the coupling is solely between the LIM3 and WW3, the configuration here also includes the ocean component of NEMO 3.6 (Madec, 2008). are taken from the same ORCA025 simulation. Regarding the atmospheric forcing, we use the latest version of the Drakkar Forcing Set (DFS 5.2, which is an updated version of the forcing set described in Brodeau et al., 2010). The choices regarding the parameterization of the wave-ice attenuation are following the ones made in the REF simulation by Ardhuin et al. (2018).
The value of the ice flexural strength has however been increased from 0.27 MPa to 0.6 MPa, which is the highest value used in Ardhuin et al. (2018). This choice makes sea ice harder to break, and aim at compensating the fact that no lateral growth of 20 sea ice is included in our coupled framework.
Three simulations are performed. First we run a simulation based solely on NEMO-LIM3 (referred to as NOT_CPL), covering the period from January 1st, 2002 to the end of 2010, in which the lateral melt parameterization from Lüpkes et al. (2012) is activated. The first years of the simulations are allowing for the adjustment of the ocean and sea ice conditions and we only 25 analyze results from August and September 2010. During that period, the sea ice extent reaches its annual minimum, providing some fetch for the generation of sea ice, and in particular in the Beaufort Sea. The model sea ice extent during the summer of 2010, and more generally the distribution of the sea ice concentration, compares reasonably well with satellite observations (not shown. Note that this period includes a drop in sea ice concentration in the Central Arctic, found both in model results and in satellite observations, that has already been documented by Zhao et al. (2018) and attributed to an enhancement of ice 30 divergence in this region this particular year. This specific period has also been chosen as some storms occurred during it, so that extreme waves conditions can be investigated. Another simulation (CPL) is initialized from NOT_CPL on August 1st 2010 and run until September 9th 2010. After that date, sea ice extent starts to increase again, and as our FSD distribution does not allow for the refreezing of the sea ice floes, we cannot represent realistically the processes at play during that period. Finally, we run a simulation over the same period, based solely on WW3 (referred to WAVE), in which the wave model is forced by sea ice conditions from the NOT_CPL simulation. In order to allow for some spin up for the waves to develop and break the ice, we remove the first 3 days. In the following, all the results are for the 37-days period between August 4th and September 9th 2010.
3 Implementation of the coupling between the wave and the sea ice models. 5 The objective of this section is to present the theoretical background and the practical implementation of the coupling between LIM3 and WW3. Fig. 2 shows the principle of the coupling and the variables that are exchanged between the two models.
Briefly, LIM3 provides sea ice floe size, thickness and concentration to WW3 in order to estimate the wave attenuation and wave-induced sea ice break-up. Note that LIM3 being a multi-category sea ice model, the actual state variable is a sea ice thickness distribution g h , from which the mean sea ice thickness can be defined either by doing a grid-cell average or by 10 doing an ice-cover average. Here we choose to exchange the ice-cover average sea ice thickness, although this choice does not affect significantly our results. WW3 then returns the WRS to LIM3, as well as the updated floe size. LIM3 takes into account the WRS in its ice transport equation, and advects the sea ice and its information on floe size. If break-up has occurred in the wave model, floe size is actualized to match the FSD assumed in WW3. The floe size is then used to estimate lateral melt.

15
In the following, we describe in more details the modifications that have been done in LIM3 and WW3 in order to couple them, and how variables are exchanged between the two models. The coupling allows a new formulation for the sea ice lateral melt in LIM3 (section 3.3).

Wave Radiative Stress
Waves transport momentum, and when they are attenuated either by dissipation or reflections, this momentum is transferred 20 to the cause of this attenuation (Longuet-Higgins, 1977). In the case of sea ice, this momentum loss thus acts as a stress that pushes sea ice in the direction of attenuated waves. Following the study of Williams et al. (2017), in which a WRS was implemented in neXtSIM, the WRS τ w,i is computed as: where ρ w is the water density, g is gravity, ω, θ and k are respectively the radial frequency, direction and wavenumber of waves 25 and S ice (x; ω, θ) is the source term corresponding to wave attenuation by sea ice at a given position.
Once estimated by WW3, the WRS is then sent to the sea ice model and added as an additional term in the momentum equation of LIM3 (Rousset et al., 2015): Fig. 1 illustrates the effect of the implementation of the WRS in our simple model. Here, the sea ice thermodynamics is switched off, so that we only simulated the effect of waves pushing sea ice. Under the action of waves, the sea ice edge shifts eastward, resulting in an increase of the sea ice concentration (panel b). As the sea ice near the sea ice edge is compacted, it creates a sharp gradient in sea ice concentration and thickness (panels b,e). When comparing panels (e) and (f), it is clear that wave attenuation also responds to this change of the sea ice properties: waves tend to penetrate further eastward when the sea 10 ice edge retreats to the east, but are then attenuated faster in the compacted sea ice.

Floe size distribution and sea ice break-up
As mentioned earlier, waves can break sea ice and thus impact the sea ice floe size. It is thus required to exchange a FSD between the two models. A FSD has been previously implemented in WW3 by Boutin et al. (2018), and is used to estimate the wave attenuation due to inelastic flexure and scattering. Following the work by Toyota et al. (2011) and Dumont et al. (2011), 15 we assume that the FSD in WW3 follows a truncated power law between a minimum floe size, D min and a maximum floe size, D max . D min corresponds to the minimum floe size that can be generated by waves and is of the order of O(10m) , while D max depends on the local waves properties and is used to estimate the level of sea ice fragmentation.
There is no FSD included in the standard version of LIM3. However, recent work by Zhang et al. (2015) and Horvat and 20 Tziperman (2015) have proposed ways to implement a FSD in sea ice models, following what is done for the sea ice thickness distribution (which is a state variable of any multi category sea ice model). Here we start by following the approach of Zhang et al. (2015), distributing ice concentration into bins corresponding to different floe sizes by defining a FSD function g D . The evolution of the FSD depends on sea ice advection, thermodynamics and mechanical processes, and is given by: 25 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 their sea ice model neXtSIM, Williams et al. (2017) have implemented a FSD that enables the floes to be advected once they have been broken by waves, making the assumption that the FSD follows a truncated power-law between a minimum and a maximum floe size, similarly to the assumption made in WW3. Here we take a similar 30 approach and implement a FSD in LIM3 that evolves following eq.(3). Yet, in contrast to Williams et al. (2017), we do not make any assumption on its shape in general, but the FSD is forced to follow the power-law assumed in WW3 as soon as wave-induced sea ice break-up occurs. This ensure coherence between the FSDs in LIM3 and WW3. We acknowledge that this assumption on the FSD is strong, and as discussed in Roach et al. (2018a), it is not a suitable way to proceed when studying the sea ice evolution, since the FSD should evolve freely and observation have regularly shown that power-law distributions are not always followed (e.g. Inoue et al., 2004). However, understanding the details of the FSD evolution is beyond the scope of this study, and assuming a power-law FSD is coherent with a distribution caused by a succession of break-up event (Toyota et al., 2011;Dumont et al., 2011). The details of the mechanical redistribution function Φ m are mostly following what has been 5 proposed by Zhang et al. (2015) and are given in appendix A.
Now that both models include a FSD, the coupling between the two models can be done in order to represent the effect of the wave-induced sea ice break-up, whose occurrence in LIM3 is determined depending on information provided by WW3. As mentioned earlier, sea ice break-up in WW3 is controlled by local wave properties and break-up events result in an update of the 10 maximum floe size D max . It is thus logical to define a similar parameter D max,LIM3 from the LIM3's FSD, that would ideally equals D max . Yet estimating D max,LIM3 is not straightforward. Indeed, our FSD implementation requires that D max,LIM3 corresponds to the upper limit of the power law followed by the FSD in both WW3 and LIM3, but also that D max,LIM3 can evolve with the deviations of the LIM3's FSD from this power-law under the effects of sea ice advection and thermodynamics.
Calling g D,P.L the distribution corresponding to a FSD following the assumed power-law, we thus define D max,LIM3 as the 15 greatest value of D for which the following condition applies: in which k Dmax is an ad hoc parameter allowing the value of D max,LIM3 to remain unchanged when the FSD slightly deviates away from the assumed power-law (after lateral melt or advection for instance). Setting k Dmax =1 is a too strong constraint, and results in noisy D max distributions, since the smallest change in the FSD after a break-up event results in a change of 20 D max,LIM3 . Values between 0.5 and 0.8 lead to smoother FSDs, but overall the choice of k Dmax does not significantly affect our results. In the following, k Dmax is set to 0.5.
Floes that have never been broken by waves have no physical reason to follow this truncated power-law. In practice, if we 25 consider a discrete number N of floe size categories, the N th category should represent these unbroken floes, with a different condition to set the value of D max,LIM3 to D N (the upper size limit of this category). We thus consider that sea ice in a grid cell can be qualified as unbroken only if most of its floes belongs to this N th category, so that D max,LIM3 = D N only if g N > 0.5c, g N being the value of the FSD function associated to the N th category and c the total sea ice concentration.

30
In all our simulations, sea ice is initialized as unbroken everywhere, so that g N = c, and D max,LIM3 = D N . As soon as wave-induced break-up occurs, D max,LIM3 is updated. To do so, the received value of D max is rounded up to the upper limit of the category it lies in. D max,LIM3 is therefore slightly greater than the value received from WW3, with an error that depends on the width of its associated floe size category.
Tests on the simplified domain were also performed to investigate the sensitivity of the results to the number and width of floes categories. This sensitivity remains really small as long as the widths of the categories are smaller than 10 m and that the categories cover a range of floe sizes larger than 300 m. In the following, we used N =60 floe size categories, that we define ) and therefore an acceptable value for the lower limit D min that the truncated power-law is assumed to follow after wave-induced break-up.
-A last category representing unbroken floes [D N −1 = 298 m, D N = 1000 m]. This value of 1000 m was set as it is one order of magnitude higher that the floe size generated by waves (Toyota et al., 2011). 15 We evaluate the effect of this part of the coupling between WW3 and LIM3, as well as the robustness of the implementation of the FSD in LIM3, by performing 2 simulations in our idealized configuration, based on WW3 only or the coupled WW3- greater than it would be if the FSD was exactly following the truncated power-law. This is because floes that have been broken to the smallest possible size do not contribute to the redistribution (see section A) and accumulate in this category since no 30 lateral growth occurs. Note that a coupled simulation in which advection had been deactivated was also run to ensure that, in a case with unaffected initial sea ice properties, no significant discrepancies were noticeable for both significant wave height and maximum floe size between an coupled and a not-coupled simulation (not shown).

Lateral melt
A parameterization to account for the sea ice lateral melt is already implemented in LIM3. Its formulation follows Steele (1992): where c is the sea ice concentration, w lat is a lateral melt rate, which depends on the difference between sea ice and sea surface 5 temperatures taken from Maykut and Perovich (1987), and α is a coefficient which varies with the floe geometry. By default, α = 0.66, which is the average value of the non-circularity of floes obtained by Rothrock and Thorndike (1984).

15
In the case of our coupled model, we estimate a FSD, and it thus makes sense to implement a parameterization of the lateral melt that depends explicitly on the FSD rather than the sea ice concentration. Following the work by Horvat and Tziperman (2015) and Roach et al. (2018a), we estimate the lateral melt as: where Φ th is the change in area covered by floes of a size D due to lateral melt (see Eq.3). Note that lateral melt for floes in the unbroken category is computed assuming that all the floes have a size D of 1000 m.
We run two simulations, in which the lateral melt is either estimated from the formulation of Lüpkes et al. (2012), or by our new formulation, which accounts for the actual FSD that is determined by both the sea ice and the wave models (Fig. 4). 25 Here we only activate the lateral melt, and turn off the basal and surface melt. The sea surface temperature is set constant to where sea ice is mostly compact and unbroken, which is likely unphysical.
4 Importance of wave-sea ice interactions.
In this section we compare the three simulations performed with the CREG025 configuration described in Section 2.2, in order to quantify the impact of the including the wave-sea ice interactions for the waves, sea ice and ocean surface properties.

5
Remember that although the coupling is only between the wave and the sea ice components, our coupled model includes an ocean component, which is only interacting with the sea ice model but not the wave model. This means that we only consider here the impact that waves may have on the ocean through their impact on the sea ice conditions.
To evaluate the impact of waves in the MIZ, we first need to define the MIZ in our model. Various criteria, relying either on 10 sea ice concentration, floe size or the region where waves impact the sea ice floe size, have been previously used to delimit the MIZ (see for instance Dumont et al., 2011;Strong and Rigor, 2013;Sutherland and Dumont, 2018). Here we take the following definition based on the maximum floe size: 0 < D max < 700 m, where < D max is the average of the maximum floe size over the studied period. Physically, it roughly corresponds to the region where sea ice has been broken during a time period that is long enough for the averaged maximum floe size to be under 1000 m (which is the limit between the broken and unbroken indicating that the wave-induced break-up is similar in the two simulations (Fig. 5f). Locally, in the Barents and Greenland seas for instance, the differences of D max can be significant, due to the specific ice drift conditions in these regions. Indeed, the overall southward drift of sea ice tends to bring unbroken sea ice from the central Arctic to regions where sea ice is broken up, increasing D max in the CPL simulation. The signs of the differences in H s and D max vary regionally. This might be due to the differences in sea ice concentration and thickness, as the wave attenuation in sea ice is very sensitive to sea ice properties (see 25 for instance Ardhuin et al., 2018). Indeed, the pattern of the differences in H s between the CPL and WAVE runs is consistent with the differences in sea ice concentration and thickness between the CPL and the NOT_CPL simulations (Fig. 6). One should keep in mind that the sea ice conditions from the NOT_CPL run are used as forcing for the WAVE run), with higher waves found in regions where ice is less concentrated and thinner.

Impact on the sea ice and sea surface properties
We now focus on the effect of adding a wave component for the sea ice properties, by comparing results from the CPL and NOT_CPL simulations. Fig. 6 shows the Pan-Arctic distribution of the sea ice thickness and concentration averaged over the 37 days considered in the CPL simulation, as well as the differences with the NOT_CPL simulation. These differences are concentrated in the vicinity of the ice edge and exhibits different signs depending on the location. Positive and negative anoma-5 lies tend to compensate, resulting in weak overall difference in sea ice extent and volume when averaging over the full Arctic Basin. If we only consider the MIZ, the sea ice volume and area decrease by about 3% and 2%, respectively, between CPL and NOT_CPL ( Fig. 7b). Locally, however, these variations can be much larger. In the MIZ of the Beaufort Sea for instance, the relative changes can be as high as 10% for mean sea ice thickness. 10 There are also difference in sea surface properties between the two simulations ( Fig. 8), with an average increase in sea surface temperature (SST) and salinity (SSS) in the MIZ of the order as high as 0.5 o C and 0.8 psu locally, respectively. It is worth noting that, in contrast to the sea ice properties, the sign of the differences in SST and SSS tends to be positive, i.e.
warmer and saltier in the CPL experiment compared to the NOT_CPL one.

Thermodynamical effect of the coupling
Given that there is no coupling between the ocean and the wave components, the difference in sea surface properties must arise from variations in sea ice conditions, and in particular the sea ice melt, that we investigate further. Fig. 9(a,b) shows the total sea ice volume melted laterally during the studied period in the CPL run as well as its difference with the same quantity from the NOT_CPL run. The sea ice volume melted by lateral melt shows very similar spatial patterns between the two sim-20 ulations, although it is estimated from two very different parameterizations (Eq. 5 and Eq. 6), although lateral melt estimated by the parameterization from Lüpkes et al. (2012) tends to be larger in NOT_CPL. The difference is substantial, the sea ice volume melted in the MIZ in NOT_CPL being 30% larger than in CPL (7a). Another signal is found in the central Arctic, where the value of lateral melt in the NOT_CPL run are small but positive. This is due to the drop in sea ice concentration that happens in the region in August 2010 (Zhao et al., 2018), resulting in a reduction of the average floe size below 100 m when ice concentration decrease and lateral melt in the NOT_CPL simulation therefore explains the deficit in sea ice concentration reported in the central Arctic when compared to the coupled simulation (Fig. 6b). Moreover, the differences in lateral melt between the two simulations being mostly negative, it cannot explain the regional patterns found in the distribution of sea ice properties differences. simulations. This result is confirmed by rerunning a coupled and an uncoupled simulation of NEMO-LIM3 while de-activating lateral melt (not shown), which yields differences in total melt distribution almost identical to the ones presented on Fig. 9(c,d).
The total sea ice volume melted once integrated over the MIZ increases by 3% between CPL and NOT_CPL, mainly due to the larger volume of sea ice melted laterally in NOT_CPL (Fig. 7a). In parallel, bottom melt slightly decreases by 1% between  20 The differences in lateral melt between the CPL and the NOT_CPL runs cannot explain the differences in sea ice and sea surface properties seen on Figs. 6 and 8. We thus investigate the impact of the WRS on the sea ice conditions and melt. Fig. 9(e,f)

Dynamical effect of the coupling
show the mean directions of the wind stress and the WRS in the CPL simulation and the ratio of WRS magnitude on wind stress respectively. This ratio is generally low, not exceeding 15% of the wind stress in the eastern Barents Sea, where the WRS reaches its highest magnitude. This is much smaller than the values retrieved from satellite observations in the Southern Ocean, 25 where the wind stress and the WRS can be of comparable magnitudes (Stopa et al., 2018a). It is also worth noting that the regions where this relative importance of the WRS compared to the wind is large do not always coincide with regions where differences in sea ice properties are significant (Fig. 6). In the Beaufort Sea for instance, there is substantially less sea ice melt in the CPL simulation than in the NOT_CPL one, although the ratios of WRS over the wind stress are only of the order of a few percents (Fig. 9f). The opposite situation is visible in the Barents Sea, where the high relative influence of the WRS does not 30 result in a significant increase of the sea ice melt when the effect of the waves is included. Therefore, the amplitude of the WRS alone does not allow to conclude on the mechanism through which the WRS impact sea ice melt. In the Southern Ocean, Stopa et al. (2018a) found that the orientation of the WRS, that tends to be orthogonal to the sea ice edge, might explain why WRS might be as important as the wind (that tends to vary much more its direction over time) to determine the position of the sea ice edge. Similarly, here, we found that the WRS is very often orientated orthogonally to the ice edge, towards packed ice. It is due to the fact that the longer waves encounter sea ice on their path, the more they are attenuated. The direction of propagating waves at a given point in sea ice is then generally imposed by the waves that have traveled the shortest distance in sea ice. This is particularly visible in some part of the Greenland and the Kara seas, where wind and wave stresses have opposite direction on average. In the Chukchi and the eastern Beaufort seas, the WRS is orthogonal to the wind stress. In contrast, in the Laptev sea, The primary effect of the WRS is to push sea ice, modifying the intensity and the direction of the sea ice drift. This impact is significant in the MIZ, where the averaged sea ice drift velocity increases by 9% between the CPL and the NOT_CPL 10 runs (Fig. 7b). This overall increase of the sea ice velocity can be explained by the fact that both WRS and sea ice drift have a dependency on wind direction. As it was the case for sea ice thickness and concentration, the distribution of the differences in sea ice drift velocity between the two simulations varies strongly depending on the region considered (not shown), but exhibits no clear relationship at the Pan-Arctic scale that could explain the differences in sea ice melt induced by the WRS.
In the following we investigate in further details the wave-sea ice interactions in two regions during storms. Indeed, although 15 the differences between the CPL and NOT_CPL run at the pan-Arctic scale remains small, it is clear that the way the waves can influence the sea ice and the ocean surface would depend on the local properties of wave, wind, sea ice and ocean surface.  (Fig. 10a,b), respectively, while they do not exceed 1 m and 7 m/s during the 3 days preceding the storm (not shown). Before the event, the south Beaufort Sea is ice-free, and the position of the sea ice edge (defined at the 15% sea ice concentration) is highly irregular, with the presence of an ice tongue centered around 72 o N and 155 o W, that is exposed upwind (and waves) on its eastern side but downwind on 25 its western side during the storm. This sea ice tongue is composed of relatively thick ice (≥ 1m). During the storm, sea ice breaks all over the ice tongue in the western part of the domain, but not further than 40 km after the sea ice edge. Both the waves and the wind stresses push the ice to the west (Fig. 10b,c), accelerating the drift that is directed north-west (Fig. 11a,c), as it was already the case before the storm (not shown). The wave action is particularly effective at the location of the sea ice tongue, where the WRS has an amplitude comparable to the wind stress over sea ice (Fig. 10c). As a consequence, the sea ice 30 drift is substantially accelerated (Fig. 11c). Considering the effect of the waves results in large changes of the sea ice thickness pattern (when comparing the CPL and NOT_CPL runs), with a decrease on the eastern part of the tongue but an increase on the western part (Fig. 11g). Outside of the sea ice tongue, the differences between the simulations are very small, likely because of the sharp sea ice thickness gradient opposing internal resistance to deformation (Fig. 11e), and the relative small effect of the WRS compared to the wind stress (Fig. 10c).
The differences of sea ice properties around the sea ice tongue between the two runs also result in changes in SST and SSS, with increase around 1 o C and 1psu, respectively, on the eastern side of the sea ice tongue and a decrease of roughly the same 5 magnitude on the western side (Fig. 12c,g). This differences arises from changes in sea ice melt, as differences of the total heat flux at the sea surface (Fig. 13a) are largely determined by bottom melt (Fig. 13b), the lateral melt contribution being one order of magnitude lower in this case. On the eastern side of the sea ice tongue, waves tend to push the sea ice away from the edge in the CPL run, and thus away from surface waters with warmer SST, resulting in a smaller amount of heat in the surface layer available for bottom melt. As the sea ice melt decreases, it also reduces the amount of freshwater received by the ocean surface, 10 resulting in larger SSS. On the western side of the south end of the ice tongue, where the sea ice is thicker in the CPL run than in the NOT_CPL one, the opposite effect happens, eventually explaining the lower SST and SSS values. One should note that the effects of this storm are particularly strong, due to the specific conditions before the storm, with warm waters brought very close to the sea ice edge during the storm (not shown) .

15
In our model, bottom melt arises from heat fluxes determined by two distinct processes: (i) a conductive heat flux, which intensity is controlled by the difference between sea ice temperature and SST, and (ii) a turbulent heat flux in the surface layer, which depends on both the SST and the shear between the sea ice and the sea surface currents . The inclusion of the effect of the waves and the WRS could in principle modified the total bottom melt through its effect on the sea ice drift, but it is not the case here, suggesting that the deficit of sea ice melt on the eastern side of the sea ice tongue in the CPL run is therefore due to 20 the combination of colder SST and the sea ice reduction.

Case 2: Storm in the Barents Sea (16-17 August 2010)
The storm that we just examined in the Beaufort Sea occurred on the same date than a second and stronger storm in the Barents Sea, with wave heights up to 5 m and south-westward winds reaching 15m/s on average over the two days (bottom panels of 25 Fig. 10d,e). In the CPL run, waves break-up sea ice over a very large area (Fig. 11f). Similarly to what we see in the Beaufort Sea, the mean direction of propagation of the waves aligns with the direction of the wind over the ice-free ocean, and is rotated orthogonally to the gradient in sea ice thickness once in the sea ice pack (Fig. 10d). The transition is however much smoother here than in the Beaufort Sea as the gradient is much weaker (Fig. 11f). In the CPL run, sea ice is drifting southward (Fig. 11b), with a slight deviation from the wind direction, and speeds twice larger than in the Beaufort Sea, due to stronger winds and 30 thinner and less concentrated sea ice.
In contrast to the effect of the storm in the Beaufort Sea, the WRS in the CPL run reaches large values (Fig. 10f). Indeed, the strong storm generates very high waves of which attenuation induces WRS close to the sea ice edge as large as the wind stress, although the WRS does not align with the direction of the wave propagation in ice. This is due to the low sea ice concentration in this region that allows for wave generation on a large region, even if partially ice-covered. The attenuation of these short in-ice generated waves dominates the WRS that is therefore aligned with the wind direction, thus accelerating the ice drift, especially close to the ice edge (Fig. 11d). 5 The differences in sea ice drift between the CPL and the NOT_CPL runs also result in differences in bottom melt (Fig. 13d), and more specifically of the part associated with the turbulent heat flux (not shown). This increase of the turbulent heat flux, which occurs in the Barents Sea but not in the Beaufort Sea, can be explained by the larger ice drift velocities driven by the WRS, which intensify the shear between the sea ice and the ocean, and therefore the turbulence in the surface mixed layer.
The differences in sea ice drift between the two runs also result in changes of the conductive heat flux. Yet, in the Barents Sea, 10 the sea ice thickness and concentrations are lower than in the Beaufort Sea while the sea ice temperature is overall higher (not shown). This results in only moderate differences of the conductive heat flux between the CPL and the NOT_CPL runs.
The differences in SST and SSS exhibit similar pattern than the differences in heat flux (Fig. 12d,h and Fig. 13c), but the magnitude of the differences are much weaker than in the Beaufort Sea, not exceeding a few tenths of o C and psu for SST and 15 SSS respectively. These small differences can be explained by two causes: (i) the small differences of sea ice properties between the two simulations result in small changes in melt, and (ii) the initial state before the storm is also different with higher SST and SSS in CPL (not shown). This difference in the initial state can be related to previous waves and wind conditions (not shown): low wind speeds are not sufficient to generate waves in the MIZ, implying that the WRS must be directed northward in the same direction as the propagating waves. It therefore compacts the sea ice edge, and thus reduces sea ice melt in the MIZ 20 in the CPL run. As seen in the Beaufort Sea case, this in turn leads to higher SST and SSS values in the vicinity of the ice edge.

What determines the impact of the waves?
From these two particular cases we suggest a generalization of the mechanisms by which the waves can impact the sea ice and ocean properties in the MIZ. It is based on a simple principle: if sea ice is moved towards warmer water, it tends to melt 25 more, and vice versa. The direction of the WRS compared to the orientation of the sea ice edge is thus fundamental if we are to understand the impact of the waves. In compact sea ice, waves are quickly attenuated and the direction of the WRS is generally towards the packed ice, thus impeding part of the sea ice melt and increasing the SST and SSS (Fig. 8). In regions where the sea ice is less concentrated and thinner, waves can be generated locally, so that the WRS aligns with the wind, whose direction determines the impact of the WRS (enhanced melt for off-ice wind and reduced melt for on-ice wind). Another key 30 factor determining the impact of the WRS onto sea ice is the internal stress of sea ice (a.k.a the rheology; see Eq.2). The impact of the WRS is larger in regions of the MIZ where the sea ice is thin and low concentrated, as the internal stress tends to be negligible (Hibler III, 1979), making the sea ice easier to deform and to drift freely. Close to the sea ice edge in the Barents Sea for instance, the WRS in storm-induced high waves conditions can be larger than the wind stress, strongly accelerating the sea ice drift towards the open ocean, which also result in an increase of the ice/ocean shear, enhancing the turbulent heat flux under sea ice and the sea ice melt.

Discussion and conclusion
The goal of this study was to examine the wave-sea ice interactions in the MIZ of the Arctic Ocean during the melt season, as these processes are thought to be important for determining the sea ice conditions but are not accounted for in the state-of- Among the wave-sea ice interaction processes considered in this study, we found that the dynamical effect of the waves (the WRS) has a larger impact than the thermodynamical one (through the additional lateral source melt). Our simulations were however limited to only a few weeks during the melting season and it is unclear if that result would hold if longer timescales 15 were considered. To make progress on this question, we would need to implement a parameterization that account for the refreezing of the floes, through lateral growth and welding. A first parameterization of that kind has been very recently developed by Roach et al. (2018a). We also anticipate that running simulation over longer time period would highlight new impacts of the WRS. Indeed, observations have revealed that heat stored during melt season below the mixed layer can significantly affect the sea ice growth the following year (Jackson et al., 2010;Timmermans, 2015). In regions where the WRS contributes to reduce 20 the ice melt, an excess of summer heat could likely accumulate under the mixed layer, possibly modulating the future evolution of the sea ice melt and growth. Recently, Smith et al. (2018) have for instance observed that a large amount of heat stored under the mixed layer could be released to melt sea ice during a storm. The significant changes of SST and SSS found locally over 37 days also highlight that wave-sea ice interactions should be considered when trying the forecast the Arctic sea ice conditions on short timescale (up to a few weeks), as these surface ocean changes can greatly affect melting and refreezing conditions. 25 The coupling developed in the present study marks a valuable new step toward an improved representation of waves and sea ice interactions in models, which might improve the representation of the dynamics of the MIZ. Yet, our coupling relies on a number of assumptions, which are most likely leading to an underestimation of the impact of the wave on the ocean and sea ice conditions. For instance, in our coupling, the sea ice rheology is unaffected by fragmentation, which is unlikely to be 30 the case (McPhee, 1980). Moreover, the sea ice model used here does not retain any memory of the past sea ice conditions, while wave would most likely affect differently sea ice that has been previously broken (Langhorne et al., 1998). Developing a similar coupling using a model that consider a state variable accounting for the previous sea ice conditions (such as the state variable 'damage' included in the sea ice model neXtSIM (Rampal et al., 2016;Williams et al., 2017)) would probably reveal new mechanisms via which waves can modulate the ocean and sea ice conditions in the MIZ. 35 Finally, the coupling we have developed here is also only considering the interactions between wave and sea ice, without any direct coupling with the ocean and the atmosphere. Yet, we know that wave dissipation would also likely impact the mixed layer, by enhancing turbulence (Couvelard et al., Submitted), and eventually modulate the rate of sea ice melt and formation (Martin and Kauffman, 1981;Rainville et al., 2011;Lee et al., 2012;Smith et al., 2018). Similarly, the effect of the waves is probably damped due to the lack of feedbacks with the atmosphere (Khon et al., 2014). Future coupling should include some 5 of these features in order to fully capture the complexity of the MIZ dynamics.
Code and data availability. Will be made available before final submission Appendix A: Floe size redistribution in the sea ice model LIM3 Here we provide the details of the calculation and implementation of the FSD, and in particular of the mechanical redistribution function Φ m that accounts for processes such as sea ice fragmentation, lead opening, ridging, and rafting. Following Zhang represents sea ice ridging and rafting, and Φ f represents the wave-induced floes fragmentation. Here we compute Φ o and Φ r in a similar way to Zhang et al. (2015), assuming that all the floes of different sizes have the same ice thickness distribution, so that changes in sea ice concentration due to open water creation or ridging affects all floes equally. As a result, the shape of the FSD and its evolution are independent from these two terms.

15
Assuming that, in a given grid cell, sea ice fragmentation does not induce any change of the sea ice concentration, Φ f can be written as (Zhang et al., 2015): where D is the floe size, Q(D) is a redistribution probability function characterizing which floes are going to be broken depending on their size, and β(D , D) is a redistribution factor quantifying the fraction of sea ice concentration transferred 20 from one floe size to another as break-up occurs. Φ f is thus used to transfer sea ice concentration from large floes to smaller floes. To ensure the conservation of sea ice area during fragmentation, β must respect (Zhang et al., 2015): In the absence of a wave model to simulate the sea state, Zhang et al. (2015) has defined β so that it redistributes uniformly the sea ice concentration of the large broken floes into the smaller floe sizes categories of the FSD. Their redistribution probability 25 function Q(D) thus assumes that a constant fraction of the sea ice cover is broken by waves during each break-up event. Their definition of Q(D) also ensures that larger floes contribute more to the redistribution than smaller floes. In our coupled model, sea ice break-up is initially computed by WW3 (for details see Boutin et al., 2018), and accounts for the sea state variability. In WW3, the FSD resulting from wave-induced break-up is assumed to follow a truncated power-law between a minimum (D min ) and a maximum (D max ) floe size. For consistency, the FSD in LIM3 after a given break-up event must follow the same power-law, defined for D taken in [D min D max ] as: where P (D > D * ) is the probability of having D > D * , and p(D) is the associated probability density. In WW3, a break-up event occurs if, firstly, waves with a wavelength λ applies a strain on sea ice greater than a given threshold, and secondly if λ/2 which is assumed to be the value of the new maximum floe size is lower than the current D max value in the wave model (Dumont et al., 2011). Therefore, a break-up event in WW3 corresponds to a decrease of D max .

10
As detailed in section 3.2, we define a maximum floe size in LIM3, D max,LIM3 , that is compared to the value of the maximum floe size received from WW3, D max,WW3 . Initially, ice is unbroken and D max,LIM3 = D max,WW3 . If break-up has occurred in WW3, then we have D max,WW3 < D max,LIM3 . In this case, D max,LIM3 must be updated to D max,WW3 value, and Φ f must be computed so that it forces the FSD in LIM3 to match the FSD assumed in WW3.
In practice in LIM3, we define a given number N of floe size categories, such that each floe size category n ∈ [0, N ] implemented in LIM3. This value is close to choices done in previous studies (see Williams et al., 2013;Bennetts et al., 2017). 20 If break-up occurs, the update of D max,LIM3 is done as follows: To force the FSD to follow this power-law during the computation of the mechanical redistribution term Φ f , in LIM3 we introduce changes in the computation of β and Q(D). When using N floe size categories, the redistribution equation (A1) becomes: Following Zhang et al. (2015), the redistribution factor β(m, n) must respect Eq.A2. β(m, n) should also allow to switch from completely unbroken ice to a truncated power-law distribution with lower limit With this choice of β(m, n), the FSD of each floe size category n < n * is equal to the distribution function derived from the power-law assumed in WW3 (g n,P.L ), given by: c being the sea ice concentration.
If sea ice in a given grid cell has already been broken, the FSD may have deviated from the truncated power-law distribution 5 (due to advection or melting). If break-up occurs again at a latter model time step, we force the FSD to be reset to the powerlaw assumed in WW3, by adjusting the fraction of each floe size category contributing to the redistribution through the value Q n . This ensures that the FSD in LIM3 and WW3 are identical. After a break-up event, D max,LIM3 is the new maximum floe size in LIM3. The sea ice contained in floe size categories associated with floes larger than D max,LIM3 is therefore entirely redistributed into smaller floe size categories by setting: The smallest floe size category (i.e D ∈ [D 0 , D 1 ]) does not contribute to the floe size redistribution, assuming that this category accounts for floes too small to be broken by waves (Toyota et al., 2011). It therefore forces Q 1 = 0. For a given floe size category n, we define ∆g th,n as the difference between the actual and theoretical values of the FSD for this floe size category (∆g th,n = g n − g n,P.L , and the theoretical value is given by the truncated power-law between D 0 and D max,LIM3 ). After the 15 redistribution of floes between categories, ∆g th,n needs to be zero, which is achieved through the adjustment of Q n in order to obtain Φ f,n = ∆g th,n . The following system thus needs to be solved: Φ f,2 = (−1 + β 2,2 )Q 2 g 2 + β 3,2 Q 3 g 3 + ... + β n * ,2 Q n * g n * + N n>n * β n≥n * ,2 g n Φ f,3 = (−1 + β 3,3 )Q 3 g 3 + β 4,3 Q 4 g 4 + ... + β n * ,n * Q n * g n * + N n>n * β n≥n * ,3 g n ...
This system consists in a triangular matrix in which all diagonal terms are non-zero. It is solved by doing: Q n * = max 0, ∆g th,n * − N n>n * β n≥n * ,n * g n g n * (β n * ,n * − 1) ... The constraint Q n > 0 ensures that the redistribution can only be done toward categories containing smaller floe size. This constraint thus implies that, in the case where ∆g th,n > 0, the FSD in LIM3 is reset to the truncated power-law only if there is enough sea ice in large floes categories to be redistributed into smaller floes categories. Besides, setting Q 1 = 0 means that the sea ice concentration associated with the smallest floe size category is never redistributed. In the absence of lateral growth, a succession of break-up events leads to an accumulation of floes in this category, deviating the FSD from the theoretical 5 power-law for floe sizes between D 0 and D 1 (see Fig. 3).        Panels (c, d, g, h) show the differences for these quantities between the CPL and NOT_CPL simulations.
Grey contours indicate the position of the ice edge determined from the averaged sea ice concentration (c = 0.15).