2.1 Idealized configuration

In order to test and illustrate the effect of the coupling (Sect. 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 by 100 grid cells, with a resolution of 0.03^∘ in both directions (corresponding roughly to 3 km). Both models are run on the same grid with the same time step set to 300 s. The coupling time step is also set to 300 s. The sea ice is forced solely by waves, without prescribing any wind or ocean currents. 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 is 0.84^∘ E). Ice thickness also increases from west to east following $h_i=\mathrm{2}(\mathrm{0.1}+{\mathrm{e}}^{-N_x/\mathrm{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^∘ latitude, and propagate to the east. The wave spectrum used as forcing at the boundary is extracted at a point south of Svalbard from an Arctic hindcast performed with WW3 described by Stopa et al. (2016). It covers the period of 2 and 3 May 2010, during which a storm occurred in this particular area (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 30 April, at 02:00 GMT, and the attenuation processes (scattering, bottom friction and inelastic flexion dissipation) use the same parameterization as in the reference simulation described by Ardhuin et al. (2018).

2.2 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 LIM3 and WW3, the configuration here also includes the ocean component of NEMO 3.6 (Madec, 2008). CREG025 encompasses the Arctic and parts of the North Atlantic down to 26^∘N and has 75 vertical levels and a nominal horizontal resolution of 0.25^∘ (≃ 12 km in the Arctic). Both NEMO-LIM3 and WW3 are run on the same grid. Initial conditions for the ocean are taken from the World Ocean Atlas 2009 climatology for temperature and salinity. The initial sea ice thickness and concentration are taken from a long ORCA025 simulation performed by the Drakkar consortium. Along the lateral open boundaries, monthly climatological conditions (comprising sea surface height, 3-D velocities, temperature and salinity) 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 ice flexural strength has however been increased from 0.27 to 0.6 MPa, which is the highest value used in Ardhuin et al. (2018). This choice makes sea ice harder to break and has been made to compensate for the fact that no lateral growth of 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 1 January 2002 to the end of 2010, in which the already-existing lateral-melt parameterization in LIM3 is activated. The first years of the simulations are allowing for the adjustment of the ocean and sea ice conditions, and we only analyse results from August and September 2010. During that period, the sea ice extent reaches its annual minimum, providing some fetch for the generation of waves, 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 divergence in this region in this particular year. This specific period has also been selected as it includes a few storm events, so extreme wave conditions can be investigated. Another simulation (CPL) is initialized from NOT_CPL on 1 August 2010 and run until 9 September 2010. After this date, sea ice extent starts to increase again, and, as our FSD does not allow for the refreezing of sea ice floes, we cannot realistically represent the processes at play during this period.

Finally, we run a simulation over the same period, based solely on WW3 (referred to as 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 exclude the first 3 d from the analysis. In the following, all the results are for the 37 d period between 4 August and 9 September 2010.

Figure 2Schematic summary of the exchanged information between the sea ice model LIM3 and the wave model WAVEWATCH III^® in our coupled framework. The two boxes correspond to the processes accounted for in a given model, while the variables exchanged between the models are listed in the bubbles.

Download

3.1 Wave radiation stress

Waves transport momentum, and, when they are attenuated by either dissipation or reflection, this momentum is transferred to what has caused 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; and $S_{\mathrm{ice}}(\mathit{x};\mathit{\omega },\mathit{\theta })$ is the source term (negative by convention in WW3) 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):

in which m is the total mass of ice and snow per unit of area; u is the ice velocity vector; σ is the internal stress tensor; f is the Coriolis parameter; k is a unit vector pointing upwards; η is the sea surface elevation; c is the sea ice concentration; and τ_a and τ_o are the atmospheric and oceanic stresses, respectively. In contrast to τ_a and τ_o, τ_w,i does not need to be multiplied by c, as the wave attenuation estimation in WW3 (and hence the WRS) is already scaled by the sea ice concentration to account for the partial sea ice cover.

Figure 1 illustrates the effect of the implementation of the WRS in our simple model. Here, the sea ice thermodynamics are switched off so that we only simulate the effect of waves pushing sea ice. Under the action of waves, the sea ice edge shifts eastward, resulting in an increase in 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, d). When comparing panels (e) and (f), it is clear that wave attenuation also responds to this change in the sea ice properties: waves tend to penetrate further eastward when the sea ice edge retreats to the east but are then attenuated faster in the compacted sea ice.

3.2 Floe size distribution and sea ice fragmentation

As mentioned earlier, waves can break sea ice and thus impact the sea ice floe size. It is therefore necessary to exchange parameters defining the FSD between the two models. An 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), 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(10 m), while D_max depends on the local wave properties and is used to estimate the level of sea ice fragmentation. As D_min is assumed to be constant in WW3, the FSD in the wave model is thus only a function of D_max. Ideally, WW3 would therefore send the value of D_max to the sea ice model, where it would be advected and updated due to the effects of thermodynamical and mechanical processes and then sent back to WW3. Yet, D_max is not an area-conserved quantity and therefore cannot be advected as a tracer (Williams et al., 2017). Instead, we thus choose to define an FSD in LIM3, from which a maximum floe size can be estimated and then sent to the wave model.

There is no FSD included in the standard version of LIM3. However, recent work by Zhang et al. (2015) and Horvat and Tziperman (2015) has proposed ways to implement an 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). In their study, Horvat and Tziperman (2015) use a joint thickness and floe size distribution in order to represent the evolution of sea ice floes affected by a great variety of processes not necessarily related to waves (e.g. welding, refreezing, ridging). The approach by Zhang et al. (2015) is simpler and computationally cheaper, as it assumes that all floes of a given size have the same ice thickness distribution, allowing the FSD to be treated independently of the sea ice thickness distribution. To do so, they hypothesize that the FSD mostly results from the fragmentation of large unbroken floes randomly yielding floes of any size smaller than the original ones. In this study, we choose to follow the simpler approach of Zhang et al. (2015), as we only consider the effects of wave-induced sea ice fragmentation and lateral melt on the FSD evolution, and our formulation of lateral melt does not depend on sea ice thickness (see Sect. 3.3). Therefore, we can consider the distribution of sea ice thickness and floes independently, defined, respectively, as

and

h being the sea ice thickness and D being the caliper diameter of the floes as defined by Rothrock and Thorndike (1984). In these definitions, A is the total area considered, and a_h and a_D are the areas within A covered by sea ice with thickness between h and h+dh and floes with diameters between D and D+dD, respectively. The evolution of the FSD depends on sea ice advection, thermodynamics and mechanical processes and is given by

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 and melt); and Φ_m is a mechanical redistribution function associated with processes like fragmentation, lead opening, ridging and rafting. In our implementation, the only processes affecting the FSD are lateral melt and sea ice fragmentation. Other processes driving the evolution of the sea ice concentration do not modify the shape of the FSD. From a technical point of view, the FSD in LIM3 is implemented as an areal distribution divided into floe size categories. It is advected in the same way as other sea ice tracers like sea ice concentration or thickness.

In their sea ice model neXtSIM, Williams et al. (2017) have also implemented an FSD but in a different way than Zhang et al. (2015). Indeed, they assume that the FSD follows a truncated power law between a minimum and maximum floe size, similar to the assumption made in WW3. Thus, in their model, the FSD always obeys a power law with a constant exponent and only evolves when sea ice fragmentation results in a reduction of the maximum floe size. Here we combine the approaches of Zhang et al. (2015) and Williams et al. (2017) and implement an FSD in LIM3 that evolves following Eq. (7) but that also undergoes a redistribution after each fragmentation event that makes it tend towards the power law assumed in WW3. Thus, in contrast to Williams et al. (2017), we do not make any assumption about the shape of the FSD in general, but, just as they do, we assume that after being fragmented by waves, the FSD follows a power law with a maximum floe size that depends on the local sea state. Our implementation only differs from the one carried out by Zhang et al. (2015) in the way we have implemented the redistribution: instead of assuming that the broken sea ice is redistributed uniformly among the smaller floes after each fragmentation event (that actually quickly tends towards a power-law distribution with a varying exponent depending on the wave state in their study), we assume that wave-induced fragmentation results in a truncated power law with a constant exponent defined for floe sizes up to the value of D_max received from WW3, under the constraint of conservation of area. As a result, the FSD resulting from fragmentation in LIM3 is the same as the one assumed in WW3, ensuring coherence between the two models. The exponent of the truncated power law corresponds to a value of about −1.85 for the cumulative-number-of-floes distribution, which is the distribution generally used to represent the FSD in the scientific literature (e.g. Toyota et al., 2011; Herman et al., 2018). This value originates from the parameterization of wave-induced fragmentation by Dumont et al. (2011) and is the one used in WW3 by Boutin et al. (2018). This redistribution is instantaneous in LIM3, which we justify by the fact that fragmentation is a violent process that can completely change an FSD on timescales of a few minutes to a few hours (see Collins et al., 2015, for a description of such an event). Note that the redistribution of the FSD due to fragmentation transfers sea ice from large floes to smaller floes. Combined with lateral melt, a process that also reduces the floe size, the action of waves in our model always reduces the floe size. The details of the mechanical redistribution function Φ_m are mostly following what has been proposed by Zhang et al. (2015) and are given in Appendix A.

Now that both models include an FSD, the two models can be coupled in order to represent the effect of the wave-induced sea ice fragmentation, the occurrence of which in LIM3 is determined depending on information provided by WW3. As mentioned earlier, sea ice fragmentation in WW3 is determined by local wave properties, and fragmentation events result in an update of the maximum floe size D_max which controls the FSD. It is thus logical to define a similar parameter D_max,LIM3 from the FSD of LIM3 that ideally equals the value of D_max in WW3. However, estimating D_max,LIM3 is not straightforward. Indeed, our FSD implementation requires not only 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 FSD in LIM3 from this power law under the effects of sea ice advection and thermodynamics. Calling $g_{D,\mathrm{P}.\mathrm{L}}$ the distribution corresponding to an FSD following the assumed power law, we thus define D_max,LIM3 as the greatest value of D for which the following condition applies:

in which $k_{D_{\max }}$ is an ad hoc parameter allowing the value of D_max,LIM3 to remain unchanged when the FSD slightly deviates from the assumed power law (after lateral melt or advection for instance). Setting $k_{D_{\max }}$ to equal 1 results in an immediate reduction of D_max if there is a negligible reduction in the proportion of large floes in the FSD after a fragmentation event (when lateral melt occurs for instance). This can be problematic, as reducing the value of D_max during a fragmentation event allows for the wave to propagate further, generating more fragmentation. Values of $k_{D_{\max }}$ between 0.5 and 0.8 allow for little variation in the FSD while keeping the same value of D_max, thus giving more weight to fragmentation than to other processes in the estimation of D_max. Overall the choice of $k_{D_{\max }}$ in this range does not significantly affect our results. In the following, $k_{D_{\max }}$ 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 consider a discrete number N of floe size categories, the Nth 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 sea ice in a grid cell as unbroken only if most of its floes belong to this Nth category, so $D_{\max ,\mathrm{LIM}\mathrm{3}}=D_N$ only if g_N>0.5c, g_N being the value of the FSD function associated with the Nth category and c being the total sea ice concentration.

In all of our simulations, sea ice is initialized as unbroken everywhere so that g_N=c, and $D_{\max ,\mathrm{LIM}\mathrm{3}}=D_N$. As soon as wave-induced fragmentation 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 with the simplified domain were also performed with different numbers and widths of floe categories to investigate the sensitivity of the results to those parameters. This sensitivity remains very small as long as the widths of the categories are smaller than 10 m and the categories cover a range of floe sizes larger than 300 m. In the following, we use N=60 floe size categories that we define as follows:

• A first category corresponds to the sea ice floes that are already broken but cannot be broken anymore (D₀=8 m, D₁=13 m). D₀ represents the smallest floe size possible in the model and is set to 8 m in order to agree with the minimum floe size used in LIM3 to estimate lateral melt from the formula by Lüpkes et al. (2012). D₀ is also on the same order as the size of the smallest floes that can be generated by wave-induced fragmentation (Mellor, 1986) and therefore is an acceptable value for the lower limit D_min that the truncated power law is assumed to follow after wave-induced fragmentation.

• There are 58 categories for which $D_n-D_{n-\mathrm{1}}=\mathrm{5}$ m, with $\mathrm{1}\le n\le N-\mathrm{1}$.

• A last category represents unbroken floes ($D_{N-\mathrm{1}}=\mathrm{298}$ m, D_N=1000 m). This value of D_N=1000 m is set as it is 1 order of magnitude higher than the floe size generated by waves (Toyota et al., 2011).

Figure 3Snapshots of D_max from the uncoupled WW3 (a) and the WW3–LIM3 (b) simulations after 72 h and the difference between the two (c). Panel (d) shows the FSD from the WW3–LIM3 run at two locations indicated with crosses on (a–c). The dashed black line in (d) corresponds to the theoretical power-law FSD assumed in WW3.

Download

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 looking at the same two simulations in our idealized configuration as presented in Fig. 1 and comparing an uncoupled WW3 simulation with a coupled WW3–LIM3 simulation (Fig. 3). Sea ice thermodynamics are switched off in the WW3–LIM3 simulation. In the uncoupled WW3 simulation, D_max evolves depending on the sea state, but sea ice thickness and concentration are constant. In the WW3–LIM3 coupled simulation, sea ice properties are all evolving as sea ice is pushed by the WRS, and the FSD is advected in LIM3. The comparison between D_max estimated from the WW3 simulation and D_max,LIM3 from the coupled simulation is shown in Fig. 3a–c. The pattern of broken sea ice is broadly similar in the two simulations (Fig. 3a, b), despite the sea ice retreat due to the WRS in the coupled case. Differences in D_max (Fig. 3c) follow the wave height differences already commented on in Fig. 1e and f. Indeed, the retreat of the ice edge due to the WRS allows for waves to propagate further with less attenuation, thus involving more sea ice fragmentation and a lower maximum floe size close to the open ocean in the coupled simulation. Further east in the MIZ, the sea ice compacted by the WRS effect generates stronger wave attenuation and thus less sea ice fragmentation and a larger maximum floe size when compared to the uncoupled simulation. Both effects partly compensate for each other, so the shift in the ice edge position has little effect on the spatial extent of broken ice. Figure 3d shows the FSD at two locations in the domain. At both locations, the distribution of ice-covered area within the different categories agrees very well between LIM3 and the truncated power law assumed in WW3. The area covered by floes of the smallest possible size in LIM3 is nevertheless greater than it would be if the FSD was exactly following the truncated power law. This is because floes that have been broken down into the smallest possible size do not contribute to the redistribution (see Appendix A) and accumulate in this category since no 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 a coupled and an uncoupled simulation (not shown).

3.3 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 the lateral-melt rate, which depends on the difference between sea ice and sea surface 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). By default, 〈D〉, which represents the average floe size, is a function of the sea ice concentration obtained empirically from observational data by Lüpkes et al. (2012):

where β=1 and c^* is introduced to avoid a singularity at c=1 and is defined as

This relationship finds a value of 〈D〉 that increases very little from its minimum value (set to D₀) as long as the sea ice concentration remains lower than ≃ 0.6 (see Fig. 3 from Lüpkes et al., 2012). In the following, we refer to this lateral-melt parameterization as the parameterization of Lüpkes et al. (2012), although we acknowledge that the work of Lüpkes et al. (2012) only provides a relationship between the average floe size and the sea ice concentration.

In the case of our coupled model, we estimate an FSD, and thus it makes sense to implement a parameterization of the lateral melt that depends explicitly on the FSD rather than on sea ice concentration. Following the work by Horvat and Tziperman (2015) and Roach et al. (2018b), we estimate the lateral melt as

where Φ_th is the change in the FSD due to lateral melt (see Eq. 7). Note that lateral melt for floes in the unbroken category is computed assuming that all the floes have a size D of 1000 m. Note also that Horvat and Tziperman (2015) and Roach et al. (2018b) represent the evolution of the floe radius and therefore include in their equation a term they call G_r, representing $\frac{\partial r}{\partial t}$, the rate at which the floe radius changes due to lateral growth or melt. In our case, we represent the evolution of the floe diameter, and the rate at which it changes due to lateral growth or melt is thus equal to $\mathrm{2}\frac{\partial r}{\partial t}$ (hence the addition of a factor of 2 in Eq. 12).

Figure 4Lateral-melt rates (estimated as percentage of sea ice concentration lost per day) estimated by the coupled model after 72 h of simulation using the parameterization of Lüpkes et al. (2012) (a) or the parameterization developed in this study accounting for wave-induced sea ice fragmentation (b). The dashed black contour in panel (b) indicates $D_{max}=\mathrm{500}$ m and thus represents the limit between broken and unbroken sea ice.

Download

Figure 5Temporal evolution of the sea ice volume loss due to lateral melt integrated over the whole domain for simulations similar to the one presented in Fig. 4 but for different values of D_min. The two colours correspond to the two lateral-melt parameterizations used in this study.

Download

We run two simulations in which the lateral melt is estimated either from the formulation of Lüpkes et al. (2012) or by our new formulation, which accounts for the FSD that is determined by both the sea ice and the wave models (Fig. 4). Here we only activate the lateral melt and turn off the basal and surface melt. The sea surface temperature is set constant to T=0.3 ^∘C. Floe size categories are the same as in Sect. 3.2. In the case of the Lüpkes et al. (2012) parameterization (Fig. 4a), lateral melt only depends on the sea ice concentration and thus follows its distribution. In the second case (Fig. 4b), lateral melt is highly constrained by both the distribution of the sea state and ice properties and is only significant where the sea ice is broken. Melt rates are overall higher when estimated from the Lüpkes et al. (2012) parameterization, mostly due to the fact that the average floe size in the uncoupled run is very close to D₀ for a wide range of concentrations. Unlike the parameterization that we propose here, the parameterization of Lüpkes et al. (2012) results in a significant lateral melt far from the ice edge, where sea ice is mostly compact and unbroken, which is likely unphysical.

We find that the results are also sensitive to the choice of D_min, regardless of the parameterization used. In the case of our FSD, the sensitivity arises from the use of D_min to compute the average diameter of the smallest floe size category, which is the category most affected by lateral melt. In order to quantify the sensitivity to the choice of D_min, we run similar experiments to the one presented in Fig. 4, varying this time the value of D_min (Fig. 5). When using the parameterization of Lüpkes et al. (2012) (Eqs. 10 and 11), the volume of sea ice melted laterally roughly doubles when D_min is divided by 2 compared to the reference simulation (using D_min=8 m). In the case of our FSD, the sensitivity is still large but greatly reduced, with an increase of only 20 % in response to the same change in D_min. Figure 5 also illustrates the strong differences in melted volume between the two parameterizations, with much less lateral melt when computed using the FSD developed here.

Figure 6Significant wave height (a) and maximum floe size (c) in the CPL simulation averaged over the period 3 August to 9 September 2010 and the differences with the WAVE simulations (b, d). The black and grey contours delimit the MIZ in the CPL simulation, defined here as $\mathrm{0}<\langle D_{max}\rangle <\mathrm{700}$ m. Note that the sea ice conditions from the NOT_CPL run are used as forcing for the WAVE run and are thus similar in the WAVE and the NOT_CPL runs.

4.1 Effect of the coupling at the pan-Arctic scale

4.1.2 Impact on the sea ice and sea surface properties

We now focus on the effect of adding a wave component on the sea ice properties by comparing results from the CPL and NOT_CPL simulations. Figure 7 shows the pan-Arctic distribution of the sea ice thickness and concentration averaged over the 37 d 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 exhibit different signs depending on the location. Positive and negative anomalies tend to compensate for each other, resulting in weak overall differences 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. 8b). 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 grid-cell-average thickness.

There are also differences in sea surface properties between the two simulations (Fig. 9), with average increases in sea surface temperature (SST) and salinity (SSS) in the MIZ as high as 0.5 ^∘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.

Figure 7Sea ice concentration and thickness in the CPL simulation (a, c) and the differences with NOT_CPL (b, d) averaged over the period 4 August 2010 to 9 September 2010. The black contours delimit the MIZ in the CPL simulation, defined here as $\mathrm{0}<\langle D_{max}\rangle <\mathrm{700}$ m.

Figure 8(a) Sea ice volume melted (in 10² km³) integrated over the MIZ and over the period between 4 August and 9 September 2010 in the CPL and the NOT_CPL simulations. Here the MIZ is defined as the region where $\mathrm{0}<\langle D_{max}\rangle <\mathrm{700}$ m in the CPL run. The contributions from lateral melt and bottom melt to the total melt for both simulations are also represented. (b) Mean sea ice volume (in 10² km³), area (in 10⁵ km²) and drift speed (in m s⁻¹) in the MIZ over the same period. Values for each simulation are found above their associated bar.

Download

4.1.3 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, and we now investigate this further. Figure 10a and b show the total sea ice volume loss from lateral melt during the study period in the CPL run as well as the difference between this and the same quantity from the NOT_CPL run. The sea ice volume melted by lateral melt shows very similar spatial patterns in the two simulations, although it is estimated from two very different parameterizations (Eqs. 9 and 12) and although lateral melt estimated by the parameterization from Lüpkes et al. (2012) tends to be lower in CPL. The difference is substantial, the sea ice volume melted in the MIZ in CPL being 30 % lower than in NOT_CPL (Fig. 8a). Another signal is found in the central Arctic, while no lateral melt occurs in this region in CPL. This signal actually arises from the combination of the drop in sea ice concentration that happens in the region in August 2010 (Zhao et al., 2018) and the use of the parameterization by Lüpkes et al. (2012) to estimate floe size and resulting lateral melt in NOT_CPL. Indeed, lower sea ice concentration values correspond to estimated average floe sizes below 100 m when estimated from the parameterization by Lüpkes et al. (2012), and thus some lateral melt is triggered. In contrast, the absence of waves in the middle of the sea ice pack in CPL leaves sea ice unbroken in this region, preventing lateral melt from occurring. An average floe size of ≃ 100 m in the middle of the pack seems somewhat unrealistic and highlights the limitations of estimating the floe size using the Lüpkes et al. (2012) parameterization in pan-Arctic configurations. Absence of lateral melt in the central Arctic in CPL explains the excess in sea ice concentration when compared to the uncoupled simulation (Fig. 7b). Moreover, the fact that differences in lateral melt between the two simulations are mostly negative means that it cannot explain the regional patterns found in the distribution of sea ice property differences.

Figure 9SST (a) and SSS (c) in the CPL run for the period between 4 August and 9 September 2010 and the difference with NOT_CPL (b, d). The black contours delimit the MIZ in the CPL simulation, defined here as $\mathrm{0}<\langle D_{max}\rangle <\mathrm{700}$ m.

Figure 10c and d show the differences in bottom and total ice melt between the CPL and NOT_CPL simulations. The spatial distributions of the differences in bottom and total ice melt are very similar, meaning that the variations in bottom melt dominate the differences in sea ice melt between CPL and NOT_CPL, although the bottom melt is computed the same way in the two simulations. This result is confirmed by rerunning a coupled and an uncoupled simulation of NEMO-LIM3 while deactivating lateral melt (not shown), which yields differences in total melt distribution almost identical to the ones presented in Fig. 10c and d.

The total melted sea ice volume, 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. 8a). In parallel, bottom melt slightly decreases by ≃ 1 % between these two simulations. This result masks the fact that the regional differences of total melt are dominated by bottom melt. An explanation is that bottom and lateral melt both depend on the available heat in the surface layer, either directly for bottom melt or indirectly through lateral melt that depends on the SST. If lateral melt occurs, it removes heat from the surface layer, therefore reducing the bottom melt capacity. Conversely, if this heat is not used for lateral melt, it remains available for bottom melt. The overall decrease in bottom melt in the MIZ between CPL and NOT_CPL visible in Fig. 8a therefore mostly results from the compensation for the increase in lateral melt due to the change in parameterization, as can be seen in Fig. 10b and c. This compensation mechanism has been reported before by Tsamados et al. (2015) and Roach et al. (2018a) and also by Roach et al. (2018b) and Bateson et al. (2019), who compare two runs of the sea ice model CICE, one with the standard lateral-melt parameterization using a constant floe size of 300 m and one using an FSD (allowing floes smaller than 300 m). In our case, this compensation is strong enough, and completely deactivating lateral melt in both runs has a negligible effect on the quantity of melted ice in our simulations.

4.1.4 Dynamical effect of the coupling

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 in Figs. 7 and 9. We thus investigate the impact of the WRS on the sea ice conditions and melt. Figure 10e and f show the mean directions of the wind stress and the WRS in the CPL simulation and the ratio of WRS magnitude to 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, where the wind stress and the WRS can be of comparable magnitude (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. 7). 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 on the order of a few percent (Fig. 10f). The opposite situation is visible in the Barents Sea, where the high relative influence of the WRS does not result in a significant increase in the sea ice melt when the effect of the waves is included. Therefore, there is no direct relationship between the intensity of the WRS and the differences in sea ice and sea surface properties between the coupled and uncoupled simulations. In the Southern Ocean, Stopa et al. (2018a) found that the orientation of the WRS, which tends to be orthogonal to the sea ice edge, might explain why WRS is as important as the wind (that tends to vary in direction much more over time) in determining 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. This 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 travelled the shortest distance in sea ice. This is particularly visible in some parts of the Greenland and Kara seas, where wind and wave stresses have opposite directions 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 directions of the WRS and the wind stress roughly align and thus work together in setting the position of the sea ice edge in the CPL run. However, at the pan-Arctic scale, there is no clear relationship between the WRS direction and the differences in sea ice melt induced by the WRS in the CPL simulation.

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 average sea ice drift velocity increases by ≃ 9 % between the CPL and the NOT_CPL runs (Fig. 8b). This overall increase in the sea ice velocity can be explained by the fact that both WRS and sea ice drift have a dependency on wind direction. As 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 detail the wave–sea ice interactions in two regions during storms. Indeed, although the differences between the CPL and NOT_CPL run at the pan-Arctic scale remain 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 the waves, wind, sea ice and ocean surface.

Figure 10Volume of sea ice melted by lateral melt in the CPL simulation over the period between 4 August and 9 September 2010 (a). Differences between the CPL and the NOT_CPL runs of lateral melt (b), bottom melt (c) and total melt (d). (e) Wind stress (red) and WRS (blue) averaged over the same period between 4 August and 9 September 2010 in the CPL simulation. Note that the WRS has been multiplied by a factor of 10 in order to improve readability. (f) Distribution of the relative magnitude of WRS over the wind stress. The grey contours represent the position of the ice edge (c=0.15) on the first day (solid line) and last day (dashed line) of the period considered in the CPL simulation.

Figure 11Significant wave height and wave mean direction of propagation (a, d), wind speed (b, e) and WRS (c, f) simulated by the CPL run during the storms that occurred in the Beaufort Sea (a–c) and in the Barents Sea (d–f) on 16–17 August 2010. The grey contours indicate the position of the sea ice edge determined from the averaged sea ice concentration (c=0.15).

4.2 Regional impacts of wave–sea ice interactions during storm events

4.2.1 Case 1 – storm in the Beaufort Sea (16–17 August 2010)

We first focus on a storm event that occurred near the MIZ in the Beaufort Sea on 16–17 August 2010 (Figs. 11a–c and 12a, e). During the storm, waves and winds are oriented towards the north-west on the west side of the domain but towards the west on the east side. Wave height and wind speed reach up to 3 and 12 m s⁻¹ (Fig. 11a, b), respectively, while they do not exceed 1 and 7 m s⁻¹, respectively, during the 3 d 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 contour) is highly irregular, with the presence of an ice tongue centred around 72^∘ N and 155^∘ W that is exposed upwind on its eastern side but downwind on its western side during the storm. This sea ice tongue is composed of relatively thick ice (≥ 1 m). 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. 11b, c), accelerating the drift that is directed north-west (Fig. 12a, c), as 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. 11c). As a consequence, the sea ice drift is substantially accelerated (Fig. 12c). The effect of the waves results in large changes in 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. 12g). 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. 12e) and the relatively small effect of the WRS compared to the wind stress (Fig. 11c).

The differences in sea ice properties around the sea ice tongue between the two runs also result in changes in SST and SSS, with increases of around 1 ^∘C and 1 psu, respectively, on the eastern side of the sea ice tongue and a decrease of roughly the same magnitude on the western side (Fig. 13c, g). These differences arise from changes in sea ice melt, as differences in the total heat flux at the sea surface (Fig. 14a) are largely determined by bottom melt (Fig. 14b), the lateral-melt contribution being 1 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 a 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, resulting in higher 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, 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).

In our model, bottom melt arises from heat fluxes determined by two distinct processes: (i) a conductive heat flux, the intensity of which 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 WRS could in principle affect the turbulent heat flux 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 the combination of a colder SST and sea ice reduction.

Figure 12Mean sea ice drift (a, b) and sea ice thickness (e, f) simulated by the CPL run during the storms that occurred in the Beaufort Sea (a, e) and in the Barents Sea (b, f) on 16–17 August 2010. 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). The dashed black contour delimits the border between broken and unbroken ice ($D_{max}=\mathrm{500}$ m).

4.2.2 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 as a second and stronger storm in the Barents Sea, with wave heights of up to 5 m and south-westerly winds reaching ≃ 15 m s⁻¹ on average over the 2 d (Fig. 11d, e). In the CPL run, waves fragment sea ice over a very large area (Fig. 12f). 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. 11d). The transition is however much smoother here than in the Beaufort Sea as the gradient is much weaker (Fig. 12f). In the CPL run, sea ice is drifting southward (Fig. 12b), with a slight deviation from the wind direction and speeds twice as high as in the Beaufort Sea, due to stronger winds and thinner and less concentrated sea ice.

Figure 13SST (a, b) and SSS (e, f) simulated by the CPL run during the storms that occurred in the Beaufort Sea (a, e) and in the Barents Sea (b, f) on 16–17 August 2010. 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).

In contrast to the effects of the storm in the Beaufort Sea, the WRS in the CPL run reaches large values (Fig. 11f). Indeed, the strong storm generates high waves, inducing a WRS as large as the wind stress close to the sea ice edge where most of the attenuation takes place, 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 over 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. 12d).

The differences in sea ice drift between the CPL and the NOT_CPL runs also result in differences in bottom melt (Fig. 14d) and more specifically in the part associated with the turbulent heat flux (not shown). This increase in 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 in the conductive heat flux. However, in the Barents Sea, the sea ice thickness and concentrations are lower than in the Beaufort Sea, while the sea ice temperature is higher overall (not shown). This results in only moderate differences of the conductive heat flux between the CPL and NOT_CPL runs.

The differences in SST and SSS exhibit similar patterns to the differences in heat flux (Figs. 13d, h and 14c), but the magnitude of the differences is much weaker than in the Beaufort Sea, not exceeding a few tenths of a degree Celsius and practical salinity units for SST and 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 wave and wind conditions: 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 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.

Figure 14Averaged differences between the CPL and NOT_CPL simulations of (a, c) the heat flux into the ocean and (b, d) the contribution to the heat flux into the ocean coming from the sea ice bottom melt during the storms in the Beaufort Sea (a, b) and in the Barents Sea (c, d) which occurred on 16–17 August 2010. Grey contours indicate the position of the ice edge determined from the averaged sea ice concentration (c=0.15).