An integrated multi-instrument methodology for studying marginal ice zone dynamics and wave-ice interactions

Wave-driven fragmentation is a key mechanism shaping the Marginal Ice Zone (MIZ). Capturing this process is therefore essential for improving sea ice models, which currently do not fully capture the complex interactions between the forcing imposed by waves and the nonlinear dynamics of the resulting sea ice breakup and deformation. To investigate these interactions, we introduce a comprehensive multi-instrument dataset from a field campaign in the MIZ of the St. Lawrence Estuary, Canada, designed to characterize wave propagation and mechanical properties of sea ice under natural forcing conditions. The dataset integrates synchronized measurements from geophone arrays, wave buoys, smartphones configured as motion sensors, and Unmanned Aerial Vehicles (UAVs), all collected during coordinated deployments across diverse ice types and sea states.

Seismic data, recorded with geophone arrays, enable estimation of the ice thickness and elastic properties via active and passive wavefield analyzes. Concurrently, wave buoys and smartphones capture ocean wave characteristics including amplitude, wavelength, and attenuation near ice edges. UAV imagery is processed with advanced methods to detect vertical ice displacements with sub-centimetre sensitivity, allowing extraction of wave dispersion relations in different ice conditions. Preliminary analyzes demonstrate strong agreement between independent measurement methods, validating the dataset's quality. This multi-sensor approach offers unique opportunities to improve our understanding of wave-ice interactions, wave attenuation, and fracture dynamics in situ, thus offering a valuable resource for the sea ice and oceanographic research community to gain insight in wave-induced ice break-up mechanisms under natural conditions.

1 Introduction

Understanding the dynamics of sea ice in the changing climate represents a major challenge requiring multidisciplinary and international collaborations, especially in light of the prospect of a seasonally ice-free Arctic. Twenty years ago, this was projected to occur around 2100 (Parry et al., 2007), with a summer minimum sea ice extent declining at a rate of 10 % per decade (Comiso et al., 2008). Recently, Kim et al. (2023) brought this timeline forward, projecting a seasonally ice-free Arctic between 2030 and 2050, regardless of greenhouse gas emission scenarios. This alarming acceleration alone is enough to justify the need to improve our understanding of sea ice dynamics. Compounding this prediction is the reality that the remaining ice in the Arctic is younger, thinner, more fragile and therefore more likely to be fragmented by ocean waves (Stroeve et al., 2012; Stroeve and Notz, 2018; Thomson, 2022; Horvat, 2022). This weakening feedback impacts the area where ocean swells or wind-generated waves can penetrate and therefore fragment the ice cover, known as the Marginal Ice Zone (MIZ). In the Arctic, the MIZ is projected to extend relative to the area covered by consolidated ice (Aksenov et al., 2017; Rolph et al., 2020; Song et al., 2025).

In the MIZ, waves break the ice into fragments, or floes, of smaller sizes. The associated evolution of the Floe Size Distribution (FSD) modifies the sea ice rheology (Dumont, 2022) – how the ice cover deforms under geophysical forcings – and triggers nonlinear feedback processes, for example by increasing melt rates and dynamical interactions with ocean currents and eddies (Horvat et al., 2016; Gupta et al., 2024). Besides that, spatial variations in ice and snow thickness modify the way light is filtered through sea ice, which influences the albedo, as well as phytoplankton bloom, both in the Arctic (Ardyna et al., 2014) and the Antarctic (Taylor et al., 2013). In this context, accurate monitoring of sea ice thickness and rheology is crucial for understanding the fragmentation of the ice cover by waves, in the view of forecasting its future state in polar to mid-latitude regions.

Large-scale monitoring of sea ice extent and concentration relies on low resolution ($\mathcal{O}\sim {\mathrm{10}}^{\mathrm{3}}-{\mathrm{10}}^{\mathrm{4}}$ m) satellite-borne microwave radiometers, while ice types, form, and thickness proxies can be derived or inferred at higher resolution from Synthetic Aperture Radar (SAR) imagery. Typical satellite-derived estimations of the ice thickness result essentially in an apparent, averaged value at the kilometre scale, because its estimation is sensitive to many local-scale factors, such as ice density, melt ponds, surface roughness, snow cover or concentration from CryoSat2 or IceSat2. In the absence of ground truth data, these lead to large uncertainties that are not compatible with mid-term predictions of the rapid evolution of sea ice. Therefore, obtaining accurate, multi-scale estimates of ice parameters will establish a new benchmark for sea ice monitoring, with impacts on many aspects of oceanographic and climate research (Kern et al., 2018).

Recently, observations of the interaction between waves and the marginal ice zone (MIZ) have been made possible – for example, in Svalbard using arrays of OpenMetBuoys (Rabault et al., 2024), or in Baffin Bay and the St. Lawrence Estuary (Canada) using unmanned aerial vehicles (UAVs) to video-record sea ice break-up caused by waves generated by the Canadian Coast Guard Ship (CCGS) Amundsen (Dumas-Lefebvre and Dumont, 2023a). The latter has proven to be a new and promising methodology for understanding wave-induced sea ice fracture. However, the quantitative assessment necessary to build, test, and develop predictive models and well-informed parametrization cannot be achieved without combining these aerial data with a characterization of the mechanical properties of the ice.

In this paper, we introduce a multi-instrument dataset that combines aerial video recordings with data measured on the ice with a variety of instruments: geophones, wave buoys and smartphones controlled remotely. With specific data processing for each instrument type, we show that it is possible to quantitatively assess ice breakup while simultaneously measuring its thickness and mechanical properties.

2 Instruments and methods

The BicWin 2024 field campaign is a series of multi-instrument acquisitions aimed at measuring simultaneously the mechanical properties of sea ice and the motion of ice induced by short-period (T<5 s) ocean waves. To this end, we have combined aerial observations (hereinafter referred to as D, for drone), seismic recordings (G for geophones), wave buoys (B) and smartphones (S). Each type of instrument is summarized in Table 1, and is described in detail in a dedicated section thereafter.

Several deployments were performed in February and March 2024 in the Saguenay Fjord as part of the Horizon glacé scientific expedition onboard the CCGS Amundsen, on Lake St-Mathieu (Fig. 1b) and in the surroundings of Baie du Ha! Ha! in the Bic National Park (Fig. 1b), located on the south shore of the St. Lawrence Estuary (Canada). Data acquisitions were performed at different dates and on different types of sea ice covers. Various configurations were used for each instrument type, to adapt the measurements to the local conditions. Table 1 summarizes the instruments used on each date, as well as the corresponding ice conditions and the qualitative sea state. We use the word wave for characterizing the sea state to indicate a significant wave height of 50 cm or larger, and calm otherwise.

Figure 1 (a) Map showing the experiment sites in the Saguenay Fjord (stations 101 and 106), Lac Saint-Mathieu and the coastal area around the Bic National Park, Rimouski, Canada. (b) A close-up view of Baie du Ha! Ha! and Anse-à-Mercier, which naturally collect sea ice pushed by dominant westerly winds. Bathymetric data are taken from the Canadian Hydrographic Service non-navigational (NONNA-10) product. The typical tidal range is about 3 m, with maximum values around 4 m.

Figure 2a shows a typical situation map with multi-instrument deployed on 26 February 2024 in Baie du Ha! Ha! between 18:00:00 and 19:30:00 UTC. The background is a geo-referenced drone picture, and the symbols mark the positions of the buoys, geophones and smartphones. In this configuration, smartphones and wave buoys were placed as close as possible from the ice edges, to record wave attenuation and fracture events. Geophones were also deployed on a line of 45 m to characterize the ice properties, while the drones recorded the ice state and the ice motion at several locations. Figure 2b shows the associated timeline, illustrating time concomitance of several types of recordings.

Table 1 Field experiments conducted during the 2024 BicWin campaign, with the corresponding ice type and sea state. See Fig. 1 for site location. Instruments: D for drone, B for buoys, G for geophone and S for smartphones.

Figure 2 (a) Situation map of the 26 February around 19:00:00 UTC. Symbols present the different instruments deployed on the field. The scalebar corresponds to a horizontal distance of 10 m. (b) Deployment timeline of the various instruments.

2.1 Geophones

The use of seismic waves to characterize the mechanical properties and the thickness of sea ice dates back to the late 1950s with pioneering works by Anderson (1958) and Hunkins (1960). However, such studies have remained rare because of the heavy logistics required in hostile, polar environments; there are therefore very few reports on in situ measurements of these parameters. Since the 2010s, however, new generations of instruments which significantly reduce the burden of seismic measurements have emerged: autonomous geophones can be installed on sea ice and record the seismic wavefield for several weeks (Moreau et al., 2020; Voermans et al., 2023), microphones can measure the air-coupled flexural wave to infer the ice thickness (Romeyn et al., 2021), and optical fibers can be deployed on floating ice for dense sampling of the wavefield (Nziengui-Bâ et al., 2022). Combined with appropriate signal processing and inversion strategies, the aforementioned authors have demonstrated the potential of these instruments to monitor the ice properties from passive seismic recordings. In itself, this achievement is of great value for future deployments of instruments able to transmit data remotely, thus providing direct measurements of the ice pack properties without the need of human intervention. During the 2024 BicWin campaign, our aim was twofold: (i) to implement a protocol to rapidly characterize sea ice from active and/or passive measurements, for all ice conditions encountered including thin, porous or soft ice; (ii) provide reference seismic waveforms to validate those extracted with other instruments such as video recordings from unmanned aerial vehicles (see Sect. 2.2), or wave buoys 2.3. In the next section, we describe the seismic array configurations used in our experiments. For our acquisitions, we used SmartSolo IGU-16HR 3C 5 Hz smart sensors, with a sampling frequency of 1000 Hz.

2.1.1 Data acquisitions

Linear arrays

Linear arrays are very convenient for extracting the dispersion relation of waves propagating in a medium. Dispersion curves relate the wavelength or wavenumber to the wave frequency, from which the phase and group velocities of the various propagating modes can be calculated at each frequency and inverted for estimating the material properties. Figure 3a shows three successive linear arrays arranged in a “U” shape. Each array contains 16 geophones separated by 3 m thus spanning a total length of 45 m. On each array side, we excite mechanical waves with local active sources, positioned at 5, 8 and 11 m away from the first and last geophones of the array, respectively. Each line is excited using 6 source positions, denoted 101 to 106 (Fig. 3b). Having multiple source positions is required for applying the methodology described in Serripierri et al. (2022). Three types of waves can propagate in a thin solid plate: shear waves, longitudinal waves, and out-of-plane waves denoted E (for east), N (for north) and Z (for vertical) in the context of a thin, horizontal ice sheet, with the north defined here as pointing in the direction of the linear array (see Fig. 3b). In order to facilitate signal processing for extracting information out of these three different modes, sources are systematically produced at each location (101 to 106) following the procedure below, always in the same “ZEN” order.

1. out-of-plane (Z) waves are generated by one person who jumps three times on the ice, each jump being a few seconds apart, in order to let waves propagate beyond the array.

2. shear (E) waves are generated by hitting the ice with a hammer in the direction perpendicular to the linear array.

3. longitudinal (N) waves are generated by hitting the ice with a hammer in the direction of the array.

At each location 101 to 106, nine sources (three in each direction) are performed. Given the large variability of ice conditions in time and space, we designed this procedure to acquire sufficient seismic signals to characterize the ice in the shortest possible time. Each acquisition with a linear array typically takes 20 min for two people, including deployment, recording and removal of the geophones. We show in Sect. 3.1 that we can extract from these recordings the averaged effective ice thickness, the Young's modulus, the Poisson's ratio and the density of the ice at each array position.

Figure 3 (a) Map showing the geographical location of the geophones for the four acquisitions on 11 February 2024: three acquisitions using linear arrays with 16 geophones (green circles), and one acquisition using a squared array with 12 geophones (blue circles). (b) Relative positions of geophones within linear arrays with respect to the position of the active sources used on each side of the array. The six positions are numbered 101 to 106. Geophones are shown as circles while sources are shown as black squares.

Two-dimensional arrays

Two-dimensional (2D) arrays enable the use of more advanced array processing techniques, such as beamforming or tomographic inversions of the ice parameters (Rost and Thomas, 2002). In these experiments, we used quadrangular arrays to test strategies for tomographic inversions of the ice properties. To this end, we used active seismic sources, following an approach similar to the linear arrays source procedure. The sources were located on a circle having the same centre as the array. Figure 3a shows an example of a 2D array deployment on 11 February 2024. Contrary to the linear array, the 2D array has not been deployed every day, but only when the ice was strong enough to support active sources along the circle. In practice, we first performed the three linear array acquisitions and then proceeded with a 2D array when time allowed it and when the ice cover allowed it. We were able to deploy 2D arrays on 10, 11 and 26 February, as well as on 6 March. On 23 and 26 February, 2D diamond-shaped arrays were deployed, for passive seismic recordings. This geometry was chosen to optimize the sensitivity and the robustness of ice properties estimates in the presence of directional noise, such as the motion induced by a swell or residual wind waves.

2.2 Unmanned aerial vehicles

Aerial observations of sea ice have been made in the Arctic and Antarctic regions since the 1970s. First used to perform reconnaissance flights to evaluate sea ice conditions and enable safer ship navigation, manned aircraft progressively became a powerful tool to study sea ice formation, evolution, interactions with surface waves (Martin and Kauffman, 1981; Wadhams, 1975) and ice floe size distribution (Toyota et al., 2006). Aerial observations provide a wide range of fields of view, with a good temporal resolution, with field validation by the aircraft crew. However, operating a manned aircraft often requires large infrastructure and reasonable weather conditions and often implies prohibitive exploitation costs.

Satellites offer an operational platform for sea ice remote sensing and space-borne imagery is used to monitor sea ice conditions on large spatial scales. Ardhuin et al. (2015) for instance managed to estimate swell propagation in the marginal ice zone using Synthetic Aperture Radar (SAR) data from the Sentinel-1A satellite. However, they underlined that significant uncertainties remain in the directional spreading of the swell. Moreover, satellite observations of sea ice often lack temporal resolution, and concomitant field observations are often missing.

The recent development of lightweight unmanned aerial vehicles (UAVs) has enabled a high spatio-temporal resolution of aerial observations at low operational costs. Moreover, UAVs can be equipped with a large range of light-weight sensors such as RGB cameras, LiDAR, multispectral, thermal, or atmospheric sensors. Over the last decade, their use has spread in many fields of cryoscience to monitor snow cover, permafrost, glaciers, or even wild life evolution (Gaffey and Bhardwaj, 2020). Recently, UAVs were used to estimate sea ice floe size distribution (Zhang and Skjetne, 2015) in the Arctic. Dumas-Lefebvre and Dumont (2023b) also used an UAV to observe sea ice break-up by ship-generated waves and analyze the resulting FSD, while Sutherland and Gascard (2016) used an UAV equipped with a scanning LiDAR to measure surface waves propagation and attenuation in a marginal ice zone.

Data acquisition

During this field campaign, we used two DJI Mavic 3 Pro to perform aerial observations. This UAV model can perform both photography and video acquisitions. The resolution varies from 12 to 20 MP. All videos were captured using the Hasselblad camera with a resolution of 3840×2160 pixels at a frame rate of 30 Hz. The 4K-video format has been calibrated using the Camera Calibrator package developed on Matlab (Bouguet, 2022) for each drone. The calibration is based on the inference of the camera optical parameter from multiple images of a checker-board taken at various angles and distances. Looking at the deformations of the checker-board pattern at different locations on the sensor, we deduce both the camera focal length and the distortion induced by the camera lenses. To perform this calibration at a distance comparable to the field conditions, we replaced the checker-board by the tilling of a large building hall at Sorbonne University. We detected several rectangular tiles arrays on different pictures to get the camera characteristics. The measured focal length is $f=\mathrm{2696}\pm \mathrm{15}$ px, and no significant optical radial distortion was found.

In the field, the flight information of the drone, including its location and orientation, and the camera state are monitored during each flight. Flight parameters are recorded either in the subtitle file associated with each video or in flight record files generated by the UAV controller. In practice, we use the drone GPS position (ϕ_D,λ_D), its vertical height H_D with respect to the water or sea ice elevation, the camera angle α to the horizontal and the yaw angle ψ defined as the angle with respect to the geographic north. The vertical height H_D is measured thanks to a high precision barometric sensor, leading to a precision of about 0.1 m, which may go up to 0.3 m in windy conditions.

During this campaign, we used UAV for (i) documenting fieldwork activities (doc), (ii) monitoring sea ice conditions using orthophotography (ortho), (iii) characterizing wave propagation in sea ice (waves), and (iv) observing wave-induced sea ice fracture (FRAC). During documentation and orthophotography operations, pictures or videos are collected without any particular restrictions on the camera pitching angle or drone height. On the other hand, movies of wave-ice interaction were captured during stationary flights above the ice edge at various heights ranging from 60 to 200 m.

The main purpose of each recordings defines the category of operation which are used to classify all images and videos recorded. The files associated to a given operation are gathered in a single folder named after the category of the performed operation (doc, ortho, waves, FRAC). A prefix and a suffix, corresponding respectively to the index over all types of operations and to the index over the operations of the same category, are added to the folder name.

2.3 Wave buoys

We used autonomous wave buoys to measure low frequency (<25 Hz) motions of ice, combining accelerometer, gyroscope, and magnetometer sensors. We deployed 6 high precision wave buoys (HPWB) developed by Peter Sutherland and his team (Veras Guimarães et al., 2018). We also tested the capacity of smartphones to be used as autonomous wave buoys (SWB for smartphone wave buoy): 6 to 30 smartphones were deployed for this purpose.

Ice motion recording with wave buoys pursued three goals. The first is to measure the wave propagation under ice as a function of time or space-dependent factors such as wind conditions, ice mechanical properties, water depth or distance to the ice edge. The second goal is to cross-validate ice motion recorded by the UAV. The third one is to test the reliability and accuracy of smartphone buoys in various wave conditions against HPWB.

We deployed buoys in arrays of various geometries to estimate wave dispersion and attenuation in the ice. Most of the data acquisition was performed using lines of wave buoys, oriented approximately in the main direction of propagation of the wave motion to focus on wave attenuation. We also deployed the instruments in other shapes, typically triangular shapes, to analyze the wave directionality. HPWB, SWB were also often deployed at the same location in several occasions, together with geophones and UAV recordings for multi-instrument cross-validation and integration.

Data acquisition

Details about the sensors and the control of the HPWB wave buoys are described in Veras Guimarães et al. (2018), while the details on the smartphone fleet are described in Zhang et al. (2025). We recall here only the main sensor characteristics. Each buoy measures simultaneously its acceleration in the three directions of space, angular velocity and orientation. The GPS signal eventually records the buoy position over time, with an horizontal accuracy of a few meters. The typical acceleration accuracy is 0.01 m s⁻². HPWB are synchronized in time through the GPS chip, while the SWB has no automatic synchronization. SWB clocks are synchronized mechanically with a typical difference of 10 ms between different SWB. We synchronize the SWB buoys with the SWB from the cross correlation function in time of a simultaneous recording of the ice motion at the same location. The smartphone temperature was maintained above 0 °C thanks to a 3D-printed insulated box that keeps the phone temperature at about 15 °C above the external temperature. We did not observe any alteration of the smartphone functionalities due to the cold.

Figure 4 Sample of the vertical acceleration a_z measured simultaneously by a wave-buoy (B2) and a smartphone (T1) on 26 February 2024. The comparison shows quantitative agreement.

Figure 4 shows a vertical acceleration signal a_z as a function of time recorded on 26 February 2024, comparing an HPWB (orange line, B2) and an SWB acquisition (blue line, T1). We observe typical wave packets, propagating under ice with a period T∼5 s. We found a quantitative agreement between the two sensors, validating the measurement of waves with smartphone sensors. The typical wave vertical acceleration amplitude, $a_z^{\prime }\sim \mathrm{0.1}$ m s⁻² here is 10 times larger than the noise level.

3 Results

This section provides some preliminary analyzes, showing promising results and demonstrating the potential of the synergy between all instruments.

3.1 Seismic data analysis

Figure 5 shows an example of seismic waveforms acquired on 6 March at the time of one vertical source at position 102, during a linear array acquisition. In the vertical channel, the flexural wave stands out, with high frequencies propagating faster than low frequencies. In the horizontal channels, however, the motions associated to shear and longitudinal waves are mixed. As the waveforms likely contain contributions from the three modes of propagation, it is therefore necessary to compute the dispersion curves to understand which modes are present. To compute the dispersion curves, we use the processing described in Serripierri et al. (2022). This processing takes advantage of the wavefield coherency through the array for the different sources. In short, it consists of applying a singular value decomposition (SVD) to the frequency spectra of all sources, before calculating the wavenumber spectra with a Fourier transform on the space variable. This significantly improves the signal-to-noise ratio (SNR) in the dispersion curves (Minonzio et al., 2010).

Figure 5 Example of 3C data recorded after a source was made in the vertical direction on 26 February 2024 at position 102. Each trace is normalized by its maximum in the recording. Black, blue and red colours correspond to the vertical (Z), shear (E) and longitudinal (N) components, respectively.

The wavenumber vs. frequency spectra shown in Fig. 6 were obtained by applying this processing to all combinations of source directions and geophones channels. Since waves propagate only in one direction through the linear array, there is no negative wavenumbers. Knowing the direction of propagation, we unwrap the spectrum beyond the spatial aliasing limit. This technique was applied to produce Fig. 6, which explains why the main branch is repeated vertically several times. Typically, the Fourier transform in the range 0 to 2 rad m⁻¹ is repeated between 2 and 4 rad m⁻¹. We eventually identify the flexural wave branch on a larger bandwidth, which decreases the uncertainty on the parameter extraction.

The polarization of the flexural wave is in the sagittal plane (defined by the normal to the ice surface and the wavenumber), and its energy is mainly on the vertical direction (Moreau et al., 2020). Hence, we expect the energy of the flexural waves to appear mainly on the spectra of the vertical channel (Z), with residual signal on the longitudinal (N) channel. The polarization of the longitudinal wave is also in the sagittal plane, but it is dominated by the contribution of the N direction. Hence, one expects the longitudinal wave to appear mainly in spectra obtained from channel N, with residual signal on the spectra calculated from channel Z. Conversely, these two types of waves should only be generated by vertical or longitudinal sources. The shear horizontal (SH) wave, on the other hand, is polarized only along the direction E. Therefore, one expects to retrieve the spectral signature of shear waves from the channel E analysis. Similarly, the shear wave should only be generated by sources in the shear (E) direction.

Spectra in Fig. 6 reveal that flexural waves are generated even with a E source, and it is visible on both N and Z channels. The longitudinal waves also appear in the spectrum of source E and channel N. We also note that the shear wave is recorded on channel E for all types of sources. To conclude, the method we use to generate different wave types from different source polarization, by either jumping onto the ice or hitting the ice with a hammer in a given direction was not fully selective. Nevertheless, the flexural, longitudinal and shear waves characteristics are contrasted enough to be separated during post-processing operations. The clear identification of the three wave types proves the efficiency of our protocol.

To recover the ice thickness, Young's modulus, Poisson's ratio and density from the dispersion curves, we use the forward model by Squire et al. (1996), based on the thin plate approximation, combined with Markov Chain Monte Carlo sampling of the parameter space, following Serripierri et al. (2022). However, instead of applying this methodology to dispersion curves obtained separately from the seismic waves generated in two opposite directions, we apply it to the average of the dispersion curves from the two propagation directions. The reason to do this is twofold. First, since there are no significant gradients of ice properties within the arrays, the averaged dispersion curves provide better SNR in the frequency-wavenumber spectra. Second, this averaging cancels undesired phase shifts that could potentially result from fluid layers with non-constant thicknesses Moreau et al. (2014), for example water layers trapped in the ice, or the snow layer. Inversion results are given in Table 2.

Figure 6 Wavenumber-frequency spectra of the acquisitions with the second linear array on 26 February 2024, for all combinations of source directions and geophone channels. Spectra were unwrapped by repeating the part of wavenumbers between 0 and 2 rad m⁻¹ from 2 to 4 rad m⁻¹, in order to identify the flexural wave on a larger bandwidth.

Table 2 Sea ice thickness h, Young's modulus E, Poisson's ratio ν and density ρ_ice inverted from seismic acquisitions of all linear arrays, compared with ice thickness h_dh measured from drill holes.

3.2 AUV data analysis

The aerial pictures are geo-rectified using a backward projection method. To do so, we use the pinhole camera model (Dawson‐Howe and Vernon, 1994) depicted in Fig. 7a. In this model, the camera optic system is simplified to a camera pupil entrance and a sensor plane, separated by a distance f corresponding to the camera focal length. In this model, all light rays entering the camera are assumed to cross at the pupil entrance, located on the optical axis of the camera. We also assume that the surface to be reconstructed is horizontal. This approximation is reasonable for ice recording motions, as the distance from the camera to the surface is much larger than the out-of-plane motions of the ice. Using geometrical transformations, we associate to each pixel P on the camera sensor a point M on the ice surface. We attribute a pixel coordinate system (x_p,y_p) to the sensor plane with the origin taken as the image top-left corner. We note (x₀,y₀) the coordinate of the camera sensor centre. We also attribute a metric coordinate system (X,Y) to the ice plane, with the origin taken at the point corresponding to the image centre. The coordinates (X,Y) of a point M on the ice surface are related to the coordinates (x_p,y_p) of a pixel P on the camera sensor by

$$\begin{array}{}\text{(1)}& \begin{array}{rl}X& ={\displaystyle \frac{(x_p-x_{\mathrm{0}})H}{f\sin \mathit{\alpha }+(y_p-y_{\mathrm{0}})\cos \mathit{\alpha }}}\\ Y& =-{\displaystyle \frac{\left(y_p-y_{\mathrm{0}}\right)H}{\sin \mathit{\alpha }\left[f\sin \mathit{\alpha }+\left(y_p-y_{\mathrm{0}}\right)\cos \mathit{\alpha }\right]}},\end{array}\end{array}$$

where α is the viewing angle with respect to the horizontal, and H the camera altitude with respect to the ice plane (see Fig. 7). To calibrate the reconstruction technique, we deployed a star-shaped array of 13 red boxes easily detectable on the drone image in Baie du Ha! Ha! on 2 March 2024. We use decameters to set the edge length of each equilateral triangular mesh to 10 m. We performed a set of seven stationary flights with different camera angles α and drone heights H. We then recorded a short video and extracted a picture of the star shape for each set of parameters (H,α). A picture captured at a height H=37.7 m and an angle α=44.8° is depicted in Fig. 7b. We manually detected each smartphone box position on the images as indicated by the red circles of Fig. 7b. The image distortion of the star-shaped array resulting from the inclined view is clearly visible.

Figure 7 (a) Schematic of the pinhole camera model. The camera is filming a plane from a height H, with an angle α to the horizontal. Point M in the field plane corresponds to a point P on the camera sensor. (b) Image of smartphones deployment captured by the Hasselbald 4K camera objective. The UAV is at a height H=37.7 m above sea ice, the camera angle to the horizontal is α=44.8°. Red circles correspond to the star-shaped array of smartphones, manually detected on the image. The plane (x_p,y_p) corresponds to pixel coordinates on the camera sensor. (c) Backward projection of image (b) on the sea ice surface where the camera sensor centre is chosen as the centre of the (x,y) coordinate system.

Using the geometrical model, we associate the pixel position of the detected smartphones to their positions (X,Y) in the metric coordinate system, neglecting the curvature of the earth. Figure 7 shows the geo-rectified image together with the smartphone positions. The star formed by the smartphone array is now regular. We applied the projection algorithm to pictures taken from different positions, and we extracted the average edge length d of the star pattern for each parameter set (H,α). We found $d=\mathrm{10.6}\pm \mathrm{0.2}$ m, compared with a target value of 10 m. The systematic error of about 5 % may be attributed to the uncertainty on the UAV height H. Note that the uncertainty will also be a function of space, in particular for small α values as we expect a larger uncertainty near the horizon line. From the geo-rectified image, the field of view can be geo-referenced by using the UAV GPS position at the record time and the drone orientation angles from the flight record. In practice, the uncertainty on the two angles introduce additional errors on the GPS location, of about 5 m for H∼100 m. Whenever possible, the detection of several geo-referenced objects such as geophones, wave buoys or smartphones in the field of view provides a better localization.

Figure 8 (a) Frame extracted from a drone flight performed over the ice edge in Baie du Ha! Ha! on 26 February 2024 with a vertical viewing angle at an altitude of 165.7 m. (b) Demodulated wave field for a frequency f=0.212 Hz. (c) Space-time spectrum $\widehat{V_x}(k,f)$ of the corrected velocity field where $k=\left|\right|\mathbf{k}\left|\right|$. The white dash line is the dispersion relation of open water waves with a water depth h_w=5.9 m. Both axes and colormaps are in logarithmic scale. (d) Linear dispersion relation f(k) extracted from the space-time spectrum. The open water dispersion relation is represented by the red dashed line. Both the x and y axes are logarithmic. The spectrum magnitude is represented by the colorbar. This figure was obtained using UAV Mesange (23waves012) on 26 February 2024 (48.34951° N, 68.81569° W) at $t_{\mathrm{0}}=\mathrm{20}:\mathrm{35}:\mathrm{59.647}$ UTC.

The drone recording of wave motions can be used to extract the dispersion relation of surface waves in the presence of ice. Here we illustrate the procedure on one drone recording, performed on 26 February 2024 in Baie du Ha! Ha!. The video was taken with a vertical view at an altitude of 165.7 m for 3 min and 12 s, at a frame rate of 30 Hz. Figure 8a shows one image in physical scale. Open water appears on the left side, while continuous fast ice is visible on the right side. Small compacted ice fragments stand in between. After geo-rectification, we define a metric coordinate system with the origin at the bottom left corner. Thanks to the texture of sea ice fragments, we use a digital image correlation (DIC) algorithm (Pan et al., 2009) to efficiently compute the velocity field ${\mathbf{V}}^{*}(x,y,t)=V_x^{*}(x,y,t){\widehat{\mathbf{e}}}_x+V_y^{*}(x,y,t){\widehat{\mathbf{e}}}_y$ seen from the drone camera point of view. DIC computes the local displacement between two images. To do so, it detects the motion of local patterns, using square boxes of dimensions W×W pixels. This correlation method is performed using the software PIVlab, developed in Matlab by Thielicke and Stamhuis (2014). To optimize the performance of the DIC, we compute the displacements between non-successive frames, separated by a number of frames Δt. In this test case, we used W=32 pixels and Δt=6 frames.

The resulting velocity field is then post-processed to eliminate spurious velocity values. As the UAV is constantly adjusting its position to compensate the surrounding air flow, an additional field ${\mathbf{V}}^{\mathrm{drone}}(x,y,t)$ appears in the computed velocity field V. Yet, this field differs from the surface wave field as the drone translational and orientational motions can be fitted by a 2D quadratic field$V_i^{\mathrm{drone}}(x,y)=a_{\mathrm{0}}+a_{\mathrm{1}}x+a_{\mathrm{2}}y+a_{\mathrm{3}}xy+a_{\mathrm{4}}{x}^{\mathrm{2}}+a_{\mathrm{5}}{y}^{\mathrm{2}}$. We correct the total displacement field by fitting and subtracting for each time step a quadratic vector field V^drone(x,y). We attribute the remaining velocity field to the field associated to the wave propagation, noted $\mathbf{V}(x,y,t)=V_x(x,y,t){\widehat{\mathbf{e}}}_x+V_y(x,y,t){\widehat{\mathbf{e}}}_y$. This field contains a broad range of waves frequencies from 0.1 to 0.6 Hz. For each frequency, we extract a complex field $\widehat{\mathbf{V}}(x,y,f)$ filtered at one frequency. Figure 8b shows the real part of the filtered wave field, for f=0.212 Hz where one can see wave fronts parallel to the ice edge with a constant wavelength of about λ=30 m over the area where sea ice is made of slush and small fragments.

Using multidimensional Fourier transform, we compute the space-time spectrum of the x-component of the velocity field $\widehat{V_x}(k,f)$, where $k=\left|\right|\mathbf{k}\left|\right|$. We then perform an average over the directions of the wave vector k, to build the velocity spectrum $\widehat{V_x}(k,f)$ presented in Fig. 8c. The spectrum exhibits a main branch corresponding to the dispersion relation of linear gravity waves in finite depth, ω²=gktanh (h_wk). The water depth was obtained using publicly available information about the local bathymetry (see Fig. 1) and the tidal elevation with respect to the chart datum. The presence of sea ice does not alter significantly the dispersion relation. However, we observe wave attenuation, as shown on the demodulated field (Fig. 8b). From the two dimensional spectrum, we extract the dispersion relation by computing the maximum of $\widehat{V_x}(k,f)$ for each frequency. The position of these maxima is refined by performing a second order polynomial fit with respect to the wavenumber k. We eventually obtain the dispersion relation as shown in Fig. 8d.

Smartphone sensors data analysis

During the field campaign, we observed large spatial and temporal variations of the wave amplitude. As an illustration, Fig. 9a shows the vertical acceleration a_z measured in 6 locations indicated in the situation map of 26 February (Fig. 2a). The colorbar encodes the distance to the ice edge, whose position was inferred from a manual detection from the UAV video recordings. The wave amplitude decreased by a factor of about 2 over a 60 m distance. This spatial attenuation can be attributed to the wave damping by the ice sheet, which is significant for both continuous and fragmented ice.

Figure 9 Wave propagation under ice. (a) Time series of vertical acceleration recorded with smartphone wave buoys (SWB) installed in a marginal ice zone in a line perpendicular to the ice edge on 26 February 2024. The colorbar indicates the distance from the first buoy. (b) Phase diagram of the wave conditions measured with the smartphones for the three days (21, 23 and 26 February 2024). Data are represented as a function of the local water depth h and the wave amplitude A. Colorbar codes the peak frequency of the associated wave spectrum.

Wave conditions experienced during the Bicwin 2024 campaign are eventually summarized on the phase diagram of Fig. 9b for the three days during which significant wave heights were observed, namely 21, 23 and 26 February. The water depth, obtained from the bathymetry and sea level data, ranges from 0.5 to 3 m while the wave amplitude varies from 2 mm to 20 cm. The peak wave frequency was mainly set by the wind conditions, with $f_{\mathrm{p}}=\mathrm{0.16}\pm \mathrm{0.02}$ Hz for 21 February, $f_{\mathrm{p}}=\mathrm{0.32}\pm \mathrm{0.05}$ Hz for 23 February and $f_{\mathrm{p}}=\mathrm{0.21}\pm \mathrm{0.03}$ Hz for 26 February.

4 Conclusions

We introduce a multi-instrument dataset designed to study mechanical wave propagation in fast ice, including its interaction with ocean swell and the propagation of seismic waves. This dataset was recorded at several locations of the St. Lawrence Estuary, Canada. Combined with original processing methodologies, it provides access to the complete process of sea ice breakup driven by gravitational waves, from the characterization of the ice's mechanical properties to the moment of breakup.

From seismic recordings, we extract the ice thickness, Young's modulus, Poisson’s ratio, and density. The accuracy of the ice thickness estimations is particularly high, while the mechanical parameters are within the typical range for sea ice, although we observe substantial variability in Young's modulus depending on the meteorological conditions during ice formation – ranging from very soft, porous ice (E=1.2 GPa) to rigid ice (E=5.3 GPa). However, in cases of very soft ice where attenuation dominates, or in partly consolidated ice where seismic waves cannot propagate, some acquisitions did not allow for parameter extraction. This is expected, since flexural motion cannot be transmitted through cracked ice.

From video recordings of aerial units, or arrays of inertial motion units (wave buoys and smartphones), we extract hydro-elastic wave propagation excited by ocean waves. We developed a processing method to extract the vertical velocity field in continuous sea ice, assuming no horizontal displacement and a flat surface. We validate this approach by comparing the local displacement with IMUs measurements. It should be noted that, when the ice breaks, horizontal motion becomes possible between ice floes, which is current a limitation of our method. However, we intend to tackle this limitation by combining drone recordings with stereo particle image velocimetry. We hope that this will also provide a new way to measure longitudinal wave propagation, not only the flexural wave. We extract the main wave characteristics from the wave spectrum (amplitude, wavelength, attenuation). The method presents a sub-centimeter sensitivity on the vertical displacement, and allows access to the wave characteristics remotely. We demonstrate that the dispersion of ocean waves is recovered. The method can be applied to various types of sea ice – from very soft fast ice to rigid fast ice – and lake ice as wells as various sea conditions, from calm to agitated. This method could also extend seismic noise interferometry approaches to video recordings on sea ice, which would be a groundbreaking advance for sea ice monitoring.

The measurements of ice properties and wave characteristics are consistent with our direct observations in the field. Using the synergy of the multi-instrument data, and by renewing yearly the acquisitions, we aim to resolve the surface deformation prior and after break-up events, in both space and time. We will thereby track the propagation of fractures in the ice cover, in order to determine the rupture threshold associated with ocean waves forcing.