Ocean surface wave radiation stress represents the flux of momentum due to the waves. When waves are dissipated or reflected by sea ice, that momentum is absorbed or reflected, resulting in a horizontal forcing that frequently compresses the ice. In this work, wave radiation stress is used to estimate the compressive force applied by waves to the marginal ice zone (MIZ). It is balanced by an ice internal compressive stress based on Mohr–Coulomb granular materials theory. The ice internal stress can be related to ice thickness, allowing this force balance to be used as a model for the estimation of MIZ ice thickness. The model was validated and tested using data collected during two field campaigns in the St. Lawrence estuary in 2016 and 2017. Modeled ice thickness was found to be consistent with the mean measured ice thickness over the conditions available. The range of validity of the model is discussed, and a definition of MIZ extent, based on the relative strength of wind and wave forcing, is proposed.
Interactions between surface waves and sea ice are becoming increasingly important in polar regions. In the Arctic, sea ice extent has reduced dramatically in recent years (e.g., Stroeve et al. 2014), and the associated increase of open water provides fetch for an emerging wave climate (Thomson and Rogers 2014). In the Antarctic, forcing by surface waves provides a stabilizing effect on the marginal ice zone (MIZ), both compacting the ice edge and resisting equatorward drift (Stopa et al. 2018).
Wave radiation stress is the flow of momentum due to surface waves (Longuet-Higgins and Stewart 1964). When waves are damped or reflected by a floating body, horizontal momentum is absorbed or reflected as well (Longuet-Higgins 1977), producing a horizontal force on the body and potentially causing movement and deformation (e.g., Christensen and Terrile 2009). Consequently, when waves are attenuated by sea ice, momentum is transferred to the ice, making wave radiation stress an important component of ice dynamics in the MIZ (e.g., Perrie and Hu 1996; Williams et al. 2017). Typically, the contribution of this stress results in the compaction of sea ice, for example, in the formation of ice bands (Wadhams 1983) or in the accumulation of grease ice (Martin and Kauffman 1981). Depending on the attenuation rate, compaction of the ice in some areas of the MIZ by wave radiation stress can be far higher (orders of magnitude) than wind or current forcing.
Horizontal compressive forces on sea ice are known to cause ice thickening. This thickening is observed across a full range of ice conditions, from frazil and shuga ice a few centimeters thick, up to the thickest ice observed north of the Canadian Arctic Archipelago. Deep within solid ice pack, this compression can cause ridging, but nearer to the ice edge, where floe sizes are typically small, compression of the ice pack causes rafting. Ice rafting, the vertical stacking of ice floes due to strong horizontal compression, can result in vertical growth rates much faster than thermal growth.
During the ice compaction process, compressive forces increase the concentration of the floes until the maximum single layer concentration is reached (typically approximately 80%; e.g., Kawaguchi and Mitsudera 2008). At this point the ice pack internal stress quickly grows, and any further compression requires floes to “ride up” upon one another, forming porous “rafts” of ice. The rafting process is assisted by wave motions that constantly rearrange the interfloe geometry. As rafts thicken, the resistance to compression increases, as does the capacity to attenuate waves. When the force required to compress the ice rafts further is equal to the force applied by external forcing, the ice stops thickening. The internal resistance to compression of rafting ice has frequently been described using Mohr–Coulomb (MC) granular materials theory (e.g., Uzuner and Kennedy 1976; Hopkins and Tuhkuri 1999; Dai et al. 2004).
In this work, ice internal stress is modeled using MC theory, and compressive stress is modeled using wave radiation stress. Wind stress and current stress are also easily included in the formulation presented. By balancing these forcing mechanisms, it is possible to estimate ice thickness in the MIZ.
Previously, Dai et al. (2004) used floe–floe collisions to estimate a pressure force due to wave motions within the MIZ. They then equated that force to MC internal stress to estimate an equilibrium ice raft thickness. In comparison, the wave radiation stress formulation for ice compaction presented in this work provides a more direct, clear, relation to wave attenuation and does not require an estimate of floe size, a value that may not be representative of a highly heterogeneous MIZ, characterized by a broad and nonnormal spectrum of sizes.
The structure of this article is as follows: The governing equations for wave radiation stress and Mohr–Coulomb internal stress are presented in section 2, and the equilibrium model is introduced in section 3. The experiment and measurements used to test the model are described in section 4. Section 5 discusses the tuning of the model, comparison of the model and the data, and the implications and limitations of the results. Section 6 is a concluding summary. Appendix A describes the dispersion relation (DR) of waves in the MIZ and the exploitation thereof for estimation of ice thickness. Appendix B details how the forcing from tidal currents has been estimated and provides a justification for neglecting this forcing in the analysis.
2. Governing equations
Here is the ice velocity, is the ice density, is the ice thickness, is the Coriolis parameter (with Earth’s angular rotation , latitude ϕ, and vertical unit vector ), g is gravitational acceleration, ξ is the large-scale sea surface height, and is the ice internal stress. The applied stress,
is the vector sum of the air-side (wind) stress , the water-side (current) stress , and the momentum flux from the wave field [defined later in Eq. (16)].
Assuming stationary equilibrium and neglecting sea surface slope, Eq. (1) reduces to a balance between applied stress and the divergence of ice internal stress:
While not universally true (e.g., in the case of wave groups or other rapidly changing conditions), this simplification allows the investigation of equilibrium dynamics. The formulations used to describe each of the terms in Eq. (3) are given in the following subsections.
a. Wind and current stress
Wind stress over sea ice can be approximated by
(e.g., McPhee 1978), where is the air density, is the air-side drag coefficient, and is the 10-m wind velocity. The air-side drag coefficient has been shown to be highly variable and to depend on ice roughness (Andreas 2011). Field measurements have returned values of approximately (e.g., Banke and Smith 1973). In this work, a value of is used.
A similar formulation can be written for drag on the ice pack by the upper ocean (e.g., Perrie and Hu 1997). Here the relevant velocity scale is the difference between the ice velocity and the mean water velocity directly below the ice . Water-side stress can then be written as
where is the water density and is the water-side drag coefficient.
b. Wave radiation stress
Wave radiation stress is defined as the excess flow of momentum due to the presence of waves (Longuet-Higgins and Stewart 1964). The wave radiation stress tensor for linear waves traveling in the x direction in water of finite depth can be written as follows (Longuet-Higgins and Stewart 1960):
Here is the group speed, and c is the phase speed. The E is the wave energy, defined as
where η is the surface elevation. The group speed and phase speed are written
The angular frequency ω and the wavenumber k are related by the dispersion relation, which in ice-free waters is written as
where h is the water depth ( appendix A discusses the dispersion relation in ice-covered water).
The terms in and in Eq. (7) can be considered to represent an isotropic pressure. In deep water these terms disappear, and the radiation stress simplifies to
in the direction of wave propagation. This is directly analogous to the standard expression for the energy flux within the wave field, (e.g., Lamb 1945, chapter IX). For waves traveling at some azimuth angle with respect to the coordinate system, the nonisotropic component of Eq. (7) is modified (Longuet-Higgins and Stewart 1964) to produce a stress tensor with the form
Equation (13) can be rewritten in spectral form to describe the radiation stress due to more typical broadbanded wave fields:
Here and c are spectral functions of k, and is the directional wavenumber spectrum defined as
Wave attenuation in sea ice is orders of magnitude larger than wave attenuation in open water (e.g., Squire et al. 1995). Assuming a steady state, and that sea ice is the dominant cause of modification of the wave field in the MIZ (neglecting bottom pressure terms and wave setup/setdown), the divergence of the wave radiation stress is approximately equivalent to the stress applied to the sea ice by the wave field:
Here i and j are indices of summation, is the i, jth component of , and is a unit vector in the jth direction.
In the example of plane waves propagating in the x direction, where all wave energy attenuation is due to dissipation by sea ice, the stress applied to the ice by the wave field simplifies to
c. Ice internal stress
Following the work of Uzuner and Kennedy (1976), MC theory has been used to describe the compression of floating ice including in ice jams, ice ridging, and ice rafting. MC theory states that the horizontal failure stress of a granular material subject to a vertical stress can be written as
where ϕ is the internal friction angle of the material. Assuming that an ice raft is floating in hydrostatic equilibrium, gravity provides the vertical forcing holding the raft together, via both the weight of the ice above the water surface and the buoyancy of the ice below the surface. The normal stress in the vertical direction at any vertical level z above the water level is equal to the weight per unit horizontal area of ice above that level:
where n is the porosity of the ice raft, the fraction of free space filled with air or water between the individual ice pieces that make up the floe jumble, is the density of the ice, and is the freeboard of the ice jumble. Similarly, at locations under the waterline, the effective vertical normal stress holding the ice jumble together is equal to the integrated buoyancy force per unit area of ice below that location (e.g., Uzuner and Kennedy 1976; Dai et al. 2004),
where is the draft of the floe jumble. It is important to note that these stresses are not absolute pressures but rather the stresses holding individual pieces of the ice raft together into one mass. The average vertical normal stress is then
where is the total ice thickness. Substituting into Eq. (18) gives an expression for the mean horizontal stress, over horizontal scales much larger than ζ, supported by the ice jumble. Multiplying that mean horizontal stress by the ice thickness then gives the horizontal force per unit length (in the cross-wave direction) G required to compress the ice raft further:
3. An equilibrium model
In equilibrium, external forcing by wave radiation stress, winds, and currents is approximately equal to the ice internal stress. At the MIZ under compressive forcing, the wave momentum flux compresses the floating pancakes until they begin to raft upon one another. This rafting continues, increasing the thickness of the compressed ice “jumble” until the stress required to compress the ice further exceeds the force applied by the wind, current, and wave stresses.
It is illustrative to consider an example of unidirectional waves traveling toward an ice edge that is orthogonal to the direction of wave propagation (Fig. 1). When in equilibrium, at any point within the ice, the ice internal stress is balanced by the integral of the external stress applied to the ice between that point and the ice edge [cf. Eq. (3)]. At ice fetch , the compressive force per unit length of ice front is then
Rearranging gives an expression for ice thickness as a function of compressive stress:
Neglecting wind and current stresses (valid for the experiments discussed in this work, but not universally) at any arbitrary distance from the ice edge, the ice thickness can then be calculated based on the incoming wave field and the local wave field by
The maximum thickness to which the ice can be compressed by wave forcing is then found by setting the in-ice wave energy to zero in Eq. (28) (all wave momentum has been applied to the ice):
This means that the maximum ice jumble thickness has a linear dependence upon significant wave height .
The Bic Winter Experiment (BicWin) 2016 and 2017 field experiments were carried out during successive Februaries, in and around the Baie du Ha! Ha! (BdHH), Parc National du Bic, Quebec, Canada. BdHH is an approximately 1-km-wide and 2-km-long embayment, located on the south shore of the lower St. Lawrence estuary and open to the east-southeast (see Fig. 2). Although the majority of the estuary is covered in ice in February on average (Galbraith et al. 2017), the area immediately to the west of BdHH up to the head of the Laurentian Channel is a sensible heat polynya induced by tidal uplifting and mixing of warm Atlantic waters near a sill (Saucier et al. 2003). The presence of this generally ice-free area means that during westerly wind events, approximately 80 km of fetch are available for wave development, and those waves then propagate directly into the bay. However, because of its geometry and relatively shallow depth, BdHH itself typically fills with ice via either advection or thermal freezing. These unique environmental conditions make the bay an excellent natural laboratory for studying wave–ice interactions.
Data from three sampling periods on three different days are presented here: 15 February 2016 (day 15), 26 February 2017 (day 26), and 27 February 2017 (day 27). During each sampling period, concurrent measurements of waves, ice thickness, currents, and winds were available. In all cases, the MIZ formed during the course of individual wind/wave events with durations of 1–2 days. The MIZ was composed of a highly irregular mixture of broken floes, pancakes, slush, and frazil (see Fig. 3). When the measurements were taken, the ice pieces had not congealed into a solid mass and were stacked into constantly rearranging jumbles many layers thick. Thermal ice growth rates were estimated and in all cases found to produce ice from one to two orders of magnitude thinner than the observed MIZ ice thickness.
The lower St. Lawrence estuary has a tidal range of approximately 4 m and strong tidal currents. These have important implications for both wave propagation in the relatively shallow BdHH and current stress on the MIZ, meaning that accurate measurements of both were necessary for this work. Bathymetry for BdHH was taken from Canadian Hydrographic Service chart 1223, Chenal du Bic, and approaches. Tidal levels were obtained from the Fisheries and Oceans Canada tide station 2995, Bic. All measurements were georeferenced in space and time, meaning that the depth for any individual measurement could be calculated by interpolation of the bathymetry and tide level, in space and time, respectively.
Before the winter freeze-up, the sea bed of the experiment site was instrumented with a Nortek acoustic wave and current (AWAC) ADCP, to measure hourly averaged currents and wave spectra, and three RBR pressure sensors. Additionally, a time-lapse camera overlooking the bay was mounted on Pic Champlain (48°19′47″N, 68°50′04″W). Georeferenced imagery from that camera was used to track the ice edge (e.g., Fig. 4a), giving the positions of all instruments relative to that ice edge.
During the experiment, an ice canoe (Lavoie and Genest 2012) was used for daily excursions to the MIZ to deploy wave buoys and to make ice thickness measurements (Accot et al. 2014; Didier et al. 2014). The vessel is approximately 8.5 m long, 1.22 m wide, with a lightship mass of 120 kg. It is human powered and capable of sliding over or through any type of icescape, including brash, highly deformed floes, or open water. The large cargo capacity and its ability to navigate the MIZ without disrupting the ice field made it ideal for this experiment.
Ice thickness was measured manually; holes were drilled in the ice floe jumble, and a meter stick with a hook on the end was placed in the hole. When the hook caught on the bottom side of the bottommost floe in the jumble, the vertical distance from the hook to the top side of the top floe was taken as the ice thickness. All measurements were geolocated and time stamped using a handheld GPS. The ice thickness was highly irregular; in some cases, over horizontal scales of just a few meters, the variability of the jumble thickness exceeded the mean thickness. Because of this irregularity, multiple measurements were required to estimate the mean thickness. For the 2016 experiment (day 15), manual measurements were not taken, and instead the effect of the ice on the dispersion relation was used to estimate thickness (see appendix A).
Wave measurements were made by placing specially designed buoys on the ice for each daily deployment. The buoys were developed by IFREMER and consisted of a data logger, a GPS, an Iridium satellite modem, and a wave sensor. Three different wave sensors were used: 1) SBG Systems (Carrières-sur-Seine, France) Ellipse-N GPS inertial navigation system, 2) Vectornav (Dallas, Texas) VN-100 inertial motion unit, and 3) ST (Geneva, Switzerland) LIS3DH accelerometer. The sensitivity of these instruments varied significantly, with the noise floor for the SBG buoys being two orders of magnitude lower than that for the ST buoys (see Fig. 5). For this reason, the least sensitive buoys were placed closest to the ice edge, where the wave signal was strongest, and careful attention was paid to the noise floors of each instrument during the analysis of wave energy. Data were analyzed in 20-min segments. Standard processing techniques (e.g., Longuet-Higgins et al. 1963; Herbers et al. 2012) were used to extract the energy spectra and first four Fourier coefficients of the directional distributions, giving mean direction and spectral spreading. Wave energy E was defined as the integral of the energy spectra over the valid range of frequencies. In the 2017 experiment, incoming wave energy was measured using wave-mode processing of the acoustic data from the bottom-mounted AWAC. During the 2016 experiment, frazil ice blocked the AWAC acoustic beams near the surface, corrupting the wave measurements. Incoming wave energy was instead estimated using the spectral parameterization of Elfouhaily et al. (1997), based on the peak period and wind forcing. That parameterization was tested for the similarly forced 27 February 2017 data, and the estimated was within 3% of that observed by the AWAC.
Over the two experiments, three deployments were usable, having strongly wave-forced and accessible MIZs. Mean conditions for those deployments are given in Table 1, approximate locations of the experiments are shown in Fig. 2, and the spatial layouts are illustrated in Fig. 4. In the 2016 deployment (Fig. 4a) the wave field entered the ice obliquely, and was defined as the distance the wave field had traveled through the ice from the ice edge to the buoy. In the 2017 deployments (Figs. 4b,c), the instrumentation was placed in a region where the waves propagated in a direction approximately orthogonal to the ice edge. Ice fetch was then defined as the shortest distance from a buoy to any location along the ice edge.
Wind speeds and directions for the experiment were taken from the Environment Canada weather station on Ile Bicquette (WMO ID 71385). They were recorded in hourly averages, which were then linearly interpolated to sample times. Wind stress was calculated by applying those measurements to Eq. (4). The values reported in Table 1 are the components of that wind stress vector in the on-ice direction.
Because of the potential stress on the ice due to tidal currents in the experimental location, current forcing was measured. Although it was unplanned, all three datasets were taken near slack tide, meaning that current forcing was low. During the 2016 experiment, the ice edge was relatively close to the AWAC location (see Fig. 4a), where water velocity was recorded by the AWAC in hourly averages. The sampling period was near slack water, and the measured current was approximately 0.16 m s in a direction of 49° true. Stress on the ice due to this measured current was calculated using Eq. (5), assuming stationary ice. Since the current flow direction was nearly parallel to the ice edge, the on-ice component of the current was N m, where the negative value means that the current was flowing from the ice toward the open water. This stress was more than two orders of magnitude smaller than the wind forcing (see Table 1) and was consequently neglected for the 15 February 2016 data. The ice edge during the 2017 experiment was much farther inside BdHH, so the currents recorded at the AWAC (0.09 m s at 220° on 26 February and 0.11 m s at 70° on 27 February) would not be expected to be representative. Instead, currents within the ice during the 2017 experiment were measured using an ice-mounted, downward-looking Nortek Aquadopp high-resolution (HR) profiler, with three orthogonal beams. During the sample periods, the currents inside the bay were variable, and the mean current magnitudes were used to calculate the stress given in Table 1. In both 2017 datasets, the current stress on the ice was more than two orders of magnitude less than the wind stress. Appendix B provides a more detailed discussion of current forcing in general in BdHH, including over the tidal cycles leading up to the events discussed here. It is worth noting that in the cases studied here, current forcing is negligible. However, that is not true universally; indeed, in the river ice-jam literature from which this work is somewhat derived (e.g., Uzuner and Kennedy 1976), current forcing is dominant.
Example energy spectra of waves propagating into the ice on 15 February 2016 are given in Fig. 5. The buoys are sufficiently far into the ice that only the energy peak remains, with higher frequencies having been already attenuated below the instrument noise levels. Consequently, it was not possible to measure the expected frequency-dependent attenuation with this buoy deployment. Nonetheless, Fig. 5 shows a clear reduction of wave energy, and a narrowing of the peak of the spectrum, from near isotropy toward the ice edge to near unidirectionality deep within the ice.
Figure 6 shows the 20-min average wave energy at each buoy, plotted as functions of . Exponential functions of the form
were fit to the data and are plotted in each figure, and the associated decay rates α are given in Table 1. In Figs. 6a and 6b, when the exponential profile is extended toward the ice edge, it reaches the level of the incoming wave energy at some point , which is before the ice edge, meaning that the exponential profile cannot continue to the ice edge. Clearly then, the actual energy decay profiles between the outermost buoys and the ice edge are not known.
The data processing in this work includes a simplification that treats the wave fields as being composed of unidirectional plane waves with wavelengths equal to the peak wavelength and wave energies . This means that the radiation stress tensor was computed from Eq. (13) instead of from Eq. (14). Because real wave fields are spread directionally and broadbanded, this simplification can cause the forcing applied to the ice to be overestimated; momentum for waves traveling at some oblique angle to the ice edge was applied directly toward the ice. Directional spectra were available for the 2017 datasets, allowing to be computed for the incoming wave field using Eq. (14) with realistic directional spreading and bandwidth. For both days, the offshore spectral peaks were narrow, meaning that compressive forcing on the ice differed only slightly from the plane wave assumption. In both cases, the ice thicknesses computed from the directional spectra were approximately 5% lower than those from the plane wave assumption. Directional spectra were not available for the 2016 data.
The mechanisms for wave attenuation in sea ice have been studied for many years, but the lack of extensive field measurements has made it difficult to determine which mechanisms are dominant under which conditions. In this work, three attenuation mechanisms are likely candidates: internal friction in the ice jumble, scattering of ice floes, and turbulent kinetic energy dissipation in the water column. This work has been based on the assumption that the internal ice friction mechanism was dominant, and it is important to note that choice has significant implications. When ice floes scatter waves, the wave momentum is reflected, and conservation of momentum requires that the momentum applied to the ice is equal to the negative of the change in momentum in the waves. This means, for example, unidirectional waves that are reflected in exactly the opposite direction from which they arrived actually impart twice the momentum to the ice than they would have if they were simply absorbed by the ice. If the full directional spectrum is measured, then Eqs. (14) and (16) correctly treat the wave stress due to reflections. However, if waves are reflected rather than absorbed by the ice, Eq. (17) gives a stress of only 50% of the true value. For this work, attenuation due to scattering was not thought likely for two reasons: 1) waves are primarily scattered by ice floes larger than the wavelength, whereas in the measurements taken here, the MIZ was composed of a slushy mixture of small ice chunks, snow, and pancake floes with maximum diameters less than approximately 10%–15% of the peak wavelengths; 2) when the scattering mechanism is dominant, the wave spectrum spreads with distance into the ice, eventually reaching isotropy. However, in all cases studied here the spectra narrowed with distance into the ice, suggesting a dissipative mechanism (e.g., Fig. 5b). Dissipation of wave energy in the water column can be due to either direct viscous dissipation or turbulent kinetic energy dissipation. The attenuation rates observed were too high for viscous dissipation to be a possibility, but TKE dissipation could have been important. When wave energy is lost to TKE dissipation, the wave momentum is transferred to the water below the ice instead of to the ice itself, meaning that Eq. (16) would overestimate the stress applied to the ice. Until the contribution of TKE dissipation to attenuation can be quantified, the modeled ice thickness estimates given here must be treated as an upper bound.
b. Wave and wind forcing: Scaling MIZ extent
For the 2017 datasets, the wave propagation direction was very nearly orthogonal to the ice edge. External forcing was calculated using Eq. (25), with wave stress from Eq. (16), wind stress from Eq. (4), and current stress from Eq. (5). Ice thickness was then calculated using Eq. (27). In all cases, the current forcing was from two to three orders of magnitude lower than the wind forcing (see Table 1) and has been neglected. The compression due to external forcing on 26 February 2017 is shown in Fig. 7a, and the resulting ice thickness is shown in Fig. 7b. That day was chosen as an example because it was the case with the highest relative importance of wind stress to wave stress. Despite this, the wind forcing caused an increase in compressive stress of less than 12% and an increase in estimated ice thickness of less than 5%. The estimated ice thickness is not additive; for example, with exclusively wind forcing the ice thickness at m would be 25 cm, and with exclusively wave forcing it would be 79 cm, but with combined wind and wave forcing the ice thickness would only increase to 83 cm.
The increase in ice compression due to wave stress depends on the wave attenuation, and once all wave energy has been attenuated by the ice, the compressive forcing due to waves reaches a constant maximum value. The wave data in Fig. 7 display this asymptotic behavior in both stress (Fig. 7a) and ice thickness (Fig. 7b), where the ice thickness approaches . Compression due to wind stress increases linearly with , resulting in an ice thickness that increases as , neither of which having a theoretical maximum value. This suggests a characteristic length scale over which the wave effects are important for setting the ice thickness; the distance from the ice edge over which compressive forcing from wave stress is greater than or equal to that from wind stress. Equivalently, it is the distance at which integrated wind stress is equal to the total wave radiation stress forcing on the ice:
For the data collected here, that means that wave forcing is dominant over the first 1–4 km from the ice edge. In deep water, Eq. (31) can be written as
Using typical Southern Ocean values, m s, m, gives the region over which waves set the ice thickness of approximately km. However, obviously conditions vary widely, and wind/wave combinations are frequently observed that would cause this MIZ definition to vary by more than an order of magnitude larger or smaller.
c. Comparison with measured ice thickness
Since in all experimental cases, xice ≪ XMIZ, the effects of wind (and currents) have been neglected from the following analysis. Measured and modeled ice thickness are given in Fig. 8. The data have been plotted in nondimensional coordinates to allow intercomparison between experiments; the ice thickness plotted on the ordinate have been divided by , calculated using Eq. (29). The nondimensional ice fetch χ is calculated by
where is the distance from the ice edge at which an exponential profile fit to the in-ice data would intersect the off-ice wave energy (cf. Fig. 6).
Assuming an exponential attenuation profile of the form given in Eq. (30), and substituting into Eq. (28) using the linear deep water dispersion relation, gives an expression for ice thickness as a function of distance from the ice edge:
which is the theoretical ice thickness curve plotted in Fig. 8. That the modeled ice thickness values estimated from the buoy measurements follow the theoretical curve is unsurprising; all three sample days showed clear exponential decay (solid black lines in Fig. 6). The fit of the measurements to the theoretical curve in Fig. 8 depends upon , which is a function of ice and water density, gravity, and the porosity and MC internal friction angle of the ice jumble. The porosity used here is , which corresponds to the porosity of packed shaken spheres and has previously been used in the ice jumble literature (Hopkins and Tuhkuri 1999). The internal friction angle was estimated by varying ϕ to minimize the least squares difference between the nondimensionalized measured and modeled ice thickness. This minimization was performed using the thickness measurements from the 2017 datasets and resulted in a value of °, with an uncertainty of approximately ° due to an estimated 1-m uncertainty in water depth. This is lower than the internal friction angles of approximately 42°–58° typically given in the sea ice literature (e.g., Mellor 1980; Hopkins and Tuhkuri 1999). However, closer inspection reveals that the majority of those laboratory experiments used ice jumbles formed from floes with a single thickness and often uniform horizontal size, in stark contrast to the slushy, multiscale floe jumbles observed in this work. The measurements presented here are, to the authors’ knowledge, the first field measurements of the Mohr–Coulomb internal friction angle in a “real world” MIZ. They are more consistent with the estimated 30°–35° for snow (e.g., Podolskiy et al. 2015), although it should be noted that the range of reported internal friction angles for snow varies from approximately 6° to 58°.
Individual measurements of jumble thickness were highly variable owing to small-scale ice irregularities, and their considerable scatter can be seen in the small crosses in Fig. 8. However, after solving for ϕ, the difference between the bin-averaged ice thickness measurements and the theoretical values was less 22% and less than 6% for all bins with more than one measurement being averaged. The standard deviation of the measured thickness divided by the modeled thickness was 0.28, highlighting the scatter in the data.
During the 2016 field experiment, two differences were present: 1) the wave direction was not directly on ice, and so it was necessary to use only the component of the stress in the on-ice direction in the ice-thickness calculation [cf. Eqs. (13)–(16)]; 2) no direct thickness measurements were taken. However, it was possible to estimate a mean ice thickness using the measured wave dispersion relation between two synchronized buoys. This thickness measurement is described in appendix A. In that case (yellow diamond in Fig. 8), the measured thickness was approximately 6% higher than the modeled thickness. However, the uncertainty was owing to the uncertainty in water depth.
In all cases, the MIZ thickness was more than an order of magnitude larger than thermal growth would have been capable of producing in the time scale of the MIZ creation (1–2 days). This means that mechanical compaction of the MIZ was the only possible way for the observed thickness to be attained. Because of the differences in magnitude between the wave forcing and the wind and current forcing in the regions studied, wave forcing is the only mechanism capable of setting ice thickness to within the variability of the measurements.
Despite the large uncertainties, the consistency between measured ice thickness and modeled ice thickness across all three experiments, particularly regarding equilibrium ice thickness, is encouraging.
d. Exponential attenuation
Exponential decay is the most commonly observed form for wave energy attenuation in the MIZ (Wadhams et al. 1988; Kohout et al. 2014, etc.), and the results presented here are generally consistent with that. Figure 6 shows wave energy as a function of ice fetch for the three deployments. In each panel the wave energy decays approximately exponentially deep inside the ice where the buoys are located. However, for two of the three cases, the exponential profile cannot extend all the way to the ice edge. Two mechanisms are potentially responsible: 1) Although stationary over the time periods of individual wave records, new frazil and pancake ice was actually still accreting to the ice edge, meaning that the outer ice was not yet fully compacted in equilibrium. 2) The ice nearer to the edge is known to be thinner than the ice deep in the pack, and wave attenuation is known to depend strongly on ice thickness (e.g., Kohout et al. 2014). This means that wave energy likely does not attenuate strictly exponentially near the ice edge or, more generally, in any conditions where the ice thickness varies substantially. In the MC ice thickness model presented here, ice thickness, wave energy attenuation, and the wave field dispersion relation are coupled. Given a functional form for the dependence of the wave energy attenuation rate on ice thickness, it would be possible to solve for the ice thickness profile. Unfortunately, a robust theory for wave energy attenuation in the near-MIZ remains elusive. Mechanisms responsible for wave attenuation depend on the type of ice, which is highly heterogeneous in the MIZ and in particular on the relationship between ice thickness, wavelength, and floe size, which is not fully captured by any existing model formulation. Furthermore, it seems likely that no single mechanism adequately captures the dynamics of the broken, jumbled ice in the near-MIZ (Sutherland and Gascard 2016).
In this work, a model for ice thickness in the MIZ, based on wave radiation stress, is proposed. As waves are attenuated by sea ice, they transfer momentum to the ice. That momentum provides a compressive force that is resisted by ice internal stress. Mohr–Coulomb theory allows the estimation of ice thickness from that internal stress, creating the link between the wave field and the ice thickness.
The model was validated with wave buoy and ice thickness measurements taken in two field campaigns in the St. Lawrence estuary. Modeled ice thickness values were generally within approximately 6% of bin-averaged thickness measurements; however, because of the highly irregular nature of the MIZ ice, the scatter in the individual measurements was large (standard deviation 28%). Additional uncertainties are present in the modeled thickness owing to differing attenuation mechanisms. The parameter space available during these experiments was limited, and it will be interesting to compare these results with larger-scale measurements in the Arctic or Antarctic MIZ.
Using this model, it is possible to calculate the maximum ice thickness due to wave forcing in the MIZ, based solely on the incoming wave field. That maximum thickness scales linearly with significant wave height . The relative importance of wind and wave forcing can be used to estimate the distance over which the waves are dominant for setting ice thickness. That distance was found to scale with and [Eq. (32)]. Waves are typically the dominant forcing mechanism from over the first order 1 km to order 100 km from the open water, suggesting that a fully coupled wave–ice–ocean numerical model in which the wave radiation stress is included may produce significantly different results near the MIZ and improve forecasting skills.
Satellite remote sensing of ice thickness is very difficult and error prone, but satellite-derived measurements of , wave direction, and are widely used and much less uncertain. The model presented here provides a method to use these more feasible parameters to estimate the ice thickness and extent of the wave-forced MIZ over large-scale polar ice margins. Such estimates would be extremely valuable for global climate observations as well as for marine operators in polar waters.
This work was partially supported by EU-FP7 project SWARP under Grant Agreement 607476, CNES project WAVE-ICE (PS), the project WAVESCALE under the “Laboratoire d’Excellence” LabexMER (ANR-10-LABX-19) cofunded by a grant from the French government under the program “Investissements d’Avenir” (PS), the Network of Centers of Excellence MEOPAR (DD), the Québec-Océan Strategic Cluster (DD), and a NSERC Discovery Grant 402257-2013 to DD. We are grateful to the staff at the Parc National du Bic for their accommodation and support during the experiment. We wish to thank Paul Nicot and Marion Bandet at UQAR for their technical and field support. We also thank Michel Hamon, Olivier Peden, and Sebastien Prigent at IFREMER, who developed and built the wave buoys. The measurements would not have been possible without the extremely hard work of our fellow ice-canoe voyagers, led by Renaud McKinnon and including Jérémy Baudry, Gwenaëlle Gremion, Sébastien Dugas, Cédric Chavanne, Claude Tremblay, Guillaume St-Onge, and Hélène Frenette.
Dispersion Relation of Waves in Ice-Covered Water
The linear dispersion relation for waves propagating in water with depth h and covered by a thin elastic plate (sea ice) with thickness ζ can be written as (e.g., Liu and Mollo-Christensen 1988)
Here g is gravitational acceleration, is the water density, is the ice density, Y is the Young’s modulus of the floating ice, s is Poisson’s ratio for the ice, and P is the compressive stress in the ice pack.
In this work, the Young’s modulus was set to zero, as floe jumbles were not rigid over the scales of interest. The pressure term, derived from the MC internal stress [Eq. (23)], was also found to have a negligible effect ( at the scales of interest). With these simplifications, the dispersion relation retains the effects of ice cover in the mass-loading term only, and can be written
The associated phase speed is then
and the group speed is
In the 2016 dataset, no direct measurements were made of ice thickness, and instead it was derived from the dispersion relation. Using the correlated signals (e.g., vertical acceleration) of two nearby sensors is a well-known method for measuring wave dispersion in sea ice (e.g., Fox and Haskell 2001; Sutherland and Rabault 2016). The phase difference of a wave signal with wavenumber vector , measured at two instruments, a and b, separated by the vector , can be written
The phase can be calculated as
where is the cross-spectra of the signals recorded at each of the two buoys.
During the 2016 sample period, the two SBG buoys were placed approximately 21 m apart, aligned with the mean wave direction. The buoy sampling times were GPS synchronized, the buoy motions were transformed into earth coordinates, and the buoy tilts in the direction of wave propagation were extracted. The cross spectra of the tilt signals from the two buoys were taken, and the phase was calculated following Eq. (A6). Equation (A5) was then solved for k. The results are plotted in Fig. A1. To calculate the ice thickness, ζ was varied in Eq. (A2) to minimize the least squares difference between the theoretical (thick black dashed line in Fig. 6) and observed (red circles in Fig. A1) dispersion relations. The best fit was found with an ice thickness of m, with uncertainties of approximately m due to uncertainties in the estimate of water depth. This ice thickness estimate corresponds to the yellow diamond in Fig. 8. A major advantage of using the dispersion relation to estimate ice thickness is that it provides an average thickness over a relatively large area of ice.
Forcing due to Tidal Currents
These experiments were conducted in a region with tidal currents. As such, it was necessary to accurately measure the current forcing on the MIZ in order to ensure that the wave forcing was indeed dominant. The ice edge during day 15 was relatively far from the shoreline, meaning that it was within the large-scale tidal flow in the basin (Fig. 2). The ice edges during days 26 and 27 were inside BdHH, a region sheltered from the main flow, with much lower current speeds.
Figure B1 shows hourly surface currents measured by the AWAC current meter deployed during BicWin 2016. On day 15, the AWAC was within approximately 200 m of the ice edge, off ice from the buoy deployment location. The main axis of the tidal current was nearly parallel to the ice edge, meaning that the on-ice current component—the component contributing to compression of the MIZ—was generally small compared to the maximum current speed. In the 5 days preceding the day 15 deployment, the maximum on-ice flow of 0.18 m s was recorded at 2100 UTC 13 February 2016, approximately 45 h before the day 15 data were taken. The compressive stress due to that current, following Eq. (5), was N m. This is only approximately 17% of the stress applied by the wind during the observation period (cf. Table 1) and has been neglected. Furthermore, the current forcing during the sample period on day 15 was very near zero (Table 1), and during the 6 h directly preceding the day 15 sampling, the current flow was in the off-ice direction, acting to decrease compressive forcing in the MIZ.
The ice edges on days 26 and 27 were inside BdHH (Fig. 2). Since no continuous measurements of current velocities inside BdHH were available, the analysis performed above for day 15 was not possible. Instead, the dependence of surface current on tidal phase and amplitude was calculated using all available current measurements from inside, and near to, BdHH. This analysis provides insight into whether strong currents at different tidal phases than those of the experiments would affect the MIZ in the days leading up to the experiments.
Tidal phase and amplitude were calculated using data from the AWAC-mounted pressure gauge (2016 and 2017), and the Department of Fisheries and Oceans tide gauge, St. Lawrence Global Observatory (SLGO) Rimouski (2015). The Hilbert transform of the surface elevation due to the tide was computed; the tidal phase was then taken to be , and the tidal range was taken to be .
In situ currents were measured using an ice-mounted Aquadopp profiler. Data from six different days are given by the filled circles in Fig. B2; each circle indicates a 120-s average horizontal velocity at a depth of approximately 20 cm below the ice. In all cases shown, the location of the Aquadopp was within 200 m of the position of the day 26 and 27 ice edges. The measurements show that the maximum under-ice mean currents were less than 7 cm s, and the corresponding stress applied to the ice, following Eq. (5), was less than N m−2, meaning that the stress applied by the current was more than one order of magnitude smaller than the wind stress (see Table 1).
Since in situ measurements were not available for all phases of the tide, surface currents were also estimated using imagery from a time-lapse camera installed on Pic Champlain (351-m elevation). The field of view for this imagery was near the mouth of BdHH, approximately 1 km from the day 26 and 27 ice edges. Particle image velocimetry (PIV) was applied to georectified images taken every 120 s showing free-drifting ice floes. Floe motions were tracked, and their displacement between consecutive images was used to estimate their velocity. A correction for wind forcing on the floes was then applied in order to estimate the component of their velocity due to currents. The velocity magnitudes at three different locations are plotted in Fig. B2. Maximal currents on 8 March 2015, during which there was a nearly complete tidal cycle with a tidal range of 3.12 m (near spring tide maximum), reached approximately 0.3 m s. This corresponds to a stress of 0.14 N m−2, which is less than the wind stress during the events studied. Moreover, this value is likely to be an overestimate of current velocities at the ice edge, particularly because of tidal forcing, as wind and waves were present during this period and their effects on floe motion could not be entirely corrected. Furthermore, the experiments of days 26 and 27 occurred farther inside BdHH, where continuity suggests that the currents should be weaker.
The weak current magnitudes over the entire tidal phase suggest that current forcing was reliably less significant than the wind forcing during the experiments. This also means that, in the near-MIZ region studied, current forcing was far smaller than forcing by the wave radiation stress.