Elastic properties of floating sea ice from air-coupled flexural waves

Air-coupled flexural waves appear as wave trains of constant frequency that arrive in advance of the direct air-wave from an impulsive source travelling over a floating ice sheet. The frequency of these waves varies with the flexural stiffness 10 of the ice sheet, which is controlled by a combination of thickness and elastic properties. We develop a theoretical framework to understand these waves, utilizing modern numerical and Fourier methods to give a simpler and more accessible description than the pioneering, yet unwieldly analytical efforts of the 1950’s. Our favoured dynamical model can be understood in terms of linear filter theory and is closely related to models used to describe the flexural waves produced by moving vehicles on floating plates. We find that air-coupled flexural waves are a robust feature of floating ice-sheets excited by impulsive sources 15 over a large range of thicknesses, and we present a simple closed-form estimator for the ice thickness. Our study is focussed on first-year sea ice of ~20-80 cm thickness in Van Mijenfjorden, Svalbard, that was investigated through active source seismic experiments over four field campaigns in 2013, 2016, 2017 and 2018. The air-coupled flexural frequencies for sea-ice in this thickness range are ~60-240 Hz. While air-coupled flexural waves for thick sea-ice have received little attention, the higher frequencies associated with thin ice on fresh water lakes and rivers are well known to the ice-skating community and have 20 been reported in popular media. Estimation of ice physical properties, following the approach we present, may allow improved surface wave modelling and wavefield subtraction in reflection seismic studies where flexural wave noise is undesirable. On the other hand, air-coupled flexural waves may also permit non-destructive continuous monitoring of ice thickness and flexural stiffness using simple, relatively inexpensive microphones located in the vicinity of the desired measurement location, either above the ice-sheet or along the shoreline. In this case, naturally forming cracks in the ice may be an appropriate impulsive 25 source capable of exciting flexural waves in floating ice sheets in a passive monitoring context.

from an impulsive source travelling over a floating ice sheet. The frequency of these waves varies with the flexural stiffness of the ice sheet, which is controlled by a combination of thickness and elastic properties. We develop a theoretical framework to understand these waves, utilizing modern numerical and Fourier methods to give a simpler and more accessible description than the pioneering, yet unwieldly analytical efforts of the 1950's. Our favoured dynamical model can be understood in terms of linear filter theory and is closely related to models used to describe the flexural waves produced by moving vehicles on 15 floating plates. We find that air-coupled flexural waves are a real and measurable component of the total wavefield robust feature of floating ice-sheets excited by impulsive sources over a large range of thicknesses, and we we present a simple closedform estimator for the ice thickness based on observable properties of the air-coupled flexural waves. Our study is focussed on first-year sea ice of ~20-80 cm thickness in Van Mijenfjorden, Svalbard, that was investigated through active source seismic experiments over four field campaigns in 2013, 2016, 2017 and 2018. The air-coupled flexural wave for the sea-ice system 20 considered in this study occurs at a constant frequency thickness product of ~48 Hz.m. Our field data includes ice ranging from ~20-80 cm thickness with corresponding air-coupled flexural frequencies from 240 Hz for the thinnest ice to 60 Hz for the thickest iceThe air-coupled flexural frequencies for sea-ice in this thickness range are ~60-240 Hz. While air-coupled flexural waves for thick sea-ice have received little attention, the readily audible, higher frequencies associated with thin ice on fresh water lakes and rivers are well known to the ice-skating community and have been reported in popular media. The 25 results of this study and further examples from lake ice, suggest the possibility of Estimation of ice physical properties, following the approach we present, may allow improved surface wave modelling and wavefield subtraction in reflection seismic studies where flexural wave noise is undesirable. On the other hand, non-contact estimation of ice thickness air-coupled flexural waves may also permit non-destructive continuous monitoring of ice thickness and flexural stiffness using simple, relatively inexpensive microphones located in the vicinity of the desired measurement location, either above the ice-sheet or 30 along the shoreline. While we have demonstrated the use of air-coupled flexural waves for ice thickness monitoring using an active source acquisition scheme, In this case, naturallynaturally forming cracks in the ice are also shown as a potential 1 Introduction 35 The term "air-coupled flexural wave" was coined by Press et al. (1951) to describe wave trains of constant frequency, varying with the flexural stiffness of a floating ice sheet, that arrive in advance of the pressure waves produced by an explosive source.
The coupling in the case of "air-coupled flexural waves" is set up between pressure waves in air and flexural waves in the solid that have phase velocity equal to the speed of sound in air. Flexural waves propagating in a floating ice sheet are a class of guided elastic waves analogous to bending waves in rods or beams and the Lamb waves of a free plate 40 . Greenhill (1886) published one of the earliest mathematical descriptions of flexural waves in a floating ice sheet and also . Greenhill (1886) published one of the earliest mathematical descriptions of flexural waves in a floating ice sheet and also discussed the concept of coupling between waves in air and water, which was further developed in Greenhill (1916). Several studies have used the dispersion of ice flexural waves to estimate ice elastic parameters, of these we highlight in particular the early work of Ewing and Crary (1934) and the study of Yang and Yates (1995) that 45 advanced the use of transform methods in this field. Passive seismic studies of flexural wave dispersion are also emerging as an effective means to continuously monitor sea-ice properties using array based wavefield transform methods (Moreau et al., 2020a) and noise interferometry and Bayesian inversion of dispersion curves (Moreau et al., 2020b).

50
Despite their commendable efforts to develop a theoretical foundation for the air-coupled flexural wave (Press and Ewing,

A convergence of fields
Coupling between air and surface or guided waves is not limited to the case of a floating ice sheet and can be anticipated in all cases where the phase velocity of the medium below the air can attain the speed of sound in air (Haskell, 1951;Novoselov et al., 2020;Press and Oliver, 1955). For example, one of the strongest air waves ever recorded 65 (Haskell, 1951;Novoselov et al., 2020;Press and Oliver, 1955). For example, one of the strongest air waves ever recorded was produced by the volcanic eruption of Krakatoa in 1883, which travelled around the globe three times. It was not until the concept of air-coupled surface waves emerged that the extent of tidal waves associated with the eruption could be adequately explained, where previous explanations had to invoke multiple coincidental local earthquakes at different locations to explain the observed arrivals of gravity waves in the ocean . 70 The correlation of the major tide-gauge disturbances with the arrival of the first or second airwave from the Krakatoa eruption has since been explained by waves propagating in a realistically coupled atmosphere-ocean system (Garrett, 1970;Harkrider and Press, 1967).

75
The idea of coupling between air and surface waves also exists in other fields, under different terminological guises. In the field of structural acoustics, the frequency at which the free bending wave becomes equal to the speed of sound in air is called the critical frequency and is closely related to the coincidence frequency, where the speeds of the free and forced bending waves are equal (Renji et al., 1997). An understanding of the coupling between air waves and bending waves in structures 80 therefore permits the design of advanced components like the honeycomb sandwich composite panels that are used in aerospace applications, where one may wish to minimise the vibration response of the panels to acoustic excitation (Renji et al., 1997).
At the coincidence frequency, a plate excited to vibration by travelling acoustic waves is acoustically transparent and maximally transmitting (Bhattacharya et al., 1971). A fascinating manifestation of this property is well known in the world of wild ice-skating. Chasing perfectly smooth, newly frozen ice, typically in the ~4-5 cm thickness range, these skaters flex the 90 ice with their bodyweight and crack it with their skates. The acoustic waves produced by this cracking propagate over the ice surface and are maximally transmitted due to constructive interference with ice flexural waves at the frequency where the phase velocity of the flexural waves in the ice is equal to the speed of sound in air. When the ice is thin, this occurs in the audible frequency range and produces distinctive tones, that to the air-coupled flexural waves of Press et al. (1951), occurring within a different part of the frequency/flexural stiffness spectrum.

A general theory -moving loads on floating plates
We now consider the closely related topic of moving vehicles travelling over floating ice sheets. The topic of moving loads o n floating plates has received considerable attention due to the importance of roads on floating sea, river and lake ice that are 100 traversed by vehicles from snow scooters to semi-trailers (Takizawa, 1988;Van der Sanden and Short, 2017;Wilson, 1955).
Aircraft take-off and landing on floating ice (Matiushina et al., 2016;Yeung and Kim, 2000) or large man-made floating structures (Kashiwagi, 2004) has also received significant attention. A major focus of many of these studies is the so-called critical load speed, where vehicle speed coincides simultaneously with the flexural phase and group velocities and the ice deflection can grow very large over time (Wilson, 1955). In general, these studies are geared towards understanding the 105 critical load speed so that vehicular transit at that speed may be avoided (e.g. Schulkes and Sneyd, 1988). Further studies have shown that the build-up of large enough deflections to break the ice may be limited by 2D spreading (Nugroho et al., 1999), viscous effects (Wang et al., 2004, acceleration/deceleration through the critical speed (Miles and Sneyd, 2003) and/or nonlinearities (Dinvay et al., 2019). Interestingly, a complementary field of study also exists that focuses on leveraging the resonant flexural waves produced by moving hovercraft to enhance their effectiveness as icebreakers (Hinchey and Colbourne, 110 1995;Kozin et al., 2017). Submarines surfacing in ice covered waters may also seek to weaken or break the ice before surfacing by creating flexural gravity waves (Kozin and Pogorelova, 2008), which has been studied via an adaption of the theoretical framework to include moving loads in the water column under a floating plate.
We propose that the air-coupled flexural wave phenomenon can be understood as a special case of the more general 115 paradigm of moving loads on floating plates, where the load speed is equal to the speed of sound in air. Furthermore, we hope to highlight the close physical relationship between phenomena spanning from Arctic seismic experiments to wild ice-skating, ice-roads, floating runways, structural acoustics and the eruption of Krakatoa. Our study was also motivated from the pragmatic standpoint that the monochromatic air-coupled flexural wave frequency is readily estimated from single sensor timeseries and can be related directly to ice thickness and rigidity. While other studies 120 have focussed on measuring the dispersion of ice flexural-waves in order to estimate ice physical properties (e.g. DiMarco et al., 1993;Moreau et al., 2020a;Yang and Yates, 1995), we present an additional alternative that is straightforward to implement and effective for both point and line explosive seismic sources.
In this study we investigate a series of active-source seismic experiments conducted in the innermost part of Van Mijenfjorden, 125 on the island of Spitsbergen, in the high-Arctic Svalbard archipelago. The experiments were conducted during the spring season, when sea-ice is best developed, in 2013, 2016, 2017 and 2018, as illustrated in Figure 1. Detonating cord of the type "Nobelcord", containing 40g Pentrit (PETN) per meter, was laid on the ice surface in various lengths or coiled into point like charges, to provide the seismic source. An example of the detonation of a point charge is shown in Figure 2. These explosive charges produce a strong air-wave that propagates radially over the ice surface, in addition to energy that travels horizontally 130 through the ice, or downwards through the water column where it may be reflected back to the surface by layers of contrasting acoustic impedance. The sum of all of these modes of propagation gives the complex seismic wavefield that was recorded using line arrays of geophones installed on the ice surface. For the majority of the experiments, the geophones were arranged in-line with the source, but we also include oblique arrangements from the 2013 campaign that employed a cross type array (see Figure 1). 135 The main purpose of these experiments was to test different acquisition designs for reflection seismic surveying of sub-seabed sediments, as reported by Johansen et al. (2019), an objective which is made difficult by the flexural wave energy propagating through the ice. Various source and receiver arrangements were tested including airguns, hydrophones and ocean bottom seismometers. However, in this study we focus on the subset of experiments employing explosive sources and vertical-140 component gimballed geophones deployed on top of the ice. This acquisition setup is operationally the simplest under field conditions, but leads to the strongest recordings of ice flexural waves. While the flexural wave energy is regarded as bothersome coherent noise from the context of sub-seabed reflection surveying, it constitutes the primary signal in the present study and is used to estimate the ice thickness.

General wavefield solution
In this study we approximate the air-ice-water system by a thin elastic plate resting on an incompressible inviscid fluid of finite 155 depth, H, that corresponds to the widely used "simplest acceptable" mathematical model advocated by Squire et al. (1996). The plate extends infinitely along the horizontal x and y axes and the vertical axis, z, is positive downwards is the plate flexural stiffness, E is Young's modulus, h is the plate thickness, σ is Poisson's ratio, ρw is the water density, ρl is the plate density, g = 9.81 m.s -2 is the acceleration due to gravity and ( , , ) is the applied external spatiotemporal force. The solution to Eq. (1) gives useful insight into the dynamics of the air coupled flexural wave and an illustrative example is given in section 5.1. 170 In order to solve for the spatiotemporal deflection of the ice surface, we apply the 3D Fourier transform (FT) defined as for an arbitrary spatiotemporal function ( , , ). The wave number vector = [ ] is decomposed in the x and y directions, 2 = | | 2 = 2 + 2 , and i denotes the imaginary unit. By applying the FT defined in Eq.
(3) highlights that the floating ice plate simply acts as a linear filter on an arbitrary spatiotemporal input signal. This formulation is also attractive because the physical significance of the different terms in the denominator are clearly preserved: 4 expresses the 180 The general solution for the resulting spatiotemporal deflection wavefield, ( , , ) , is then given by the full threedimensional inverse FT of ( , , ) as 185 . (4)

Ring delta function source for radial symmetry
To exploit spatial symmetry and reduce the dimensionality of the problem, we approximate the explosive source as a ringshaped Dirac delta pulse expanding outwards from the origin with radial symmetry. This allows us to reduce the spatial dimensionality of the problem by using the radial coordinate = √ 2 + 2 and representing the propagating air wave source 190 using the ring delta function, where is the speed of sound in air. By performing a partial Fourier transform of Eq. (5) with respect to time t, we obtain the spatio-frequency representation of the source as 195 so the spatio-frequency deflection wavefield can be written as where 0 ( ) = ( ) 4 + − ℎ 2 − ( ) coth ( ).
The spatiotemporal deflection wavefield in radial coordinates then becomes 200 When evaluating Eq. (9) numerically, we substitute 0 ( ) → 0 ( + ), where > 0, which avoids dividing by zero and gives a tuneable heuristic damping parameter for wave attenuation and dissipation.

Dispersion relation
The dispersion relation for the propagating wavefield in the ice sheet is given by the wavenumber and frequency pairs that 205 satisfy the equation The solution of Eq. (10) gives the general dispersion relation (e.g., Squire et al., 1996) This dispersion relation is complete in the sense that it retains all physical mechanisms included in the dynamical model Eq. 210 (1), and it is well known from previous studies (e.g. Greenhill, 1886;Squire et al., 1996). At typical vehicle speeds and wavelengths that have been the main focus of moving load on floating plate studies, the plate acceleration may be safely neglected (e.g. Schulkes and Sneyd, 1988;Wang et al., 2004) giving the approximate dispersion relation Another common assumption is to assume infinite water depth, H, which causes the hyperbolic tangent term to approach unity. 220 However, when the load moves at the speed of sound in air, plate flexural waves with much shorter  compares the full dispersion relation (Eq. (11)) with the approximation that neglects water depth and plate acceleration (Eq. (12)). Air-coupled waves are excited where the phase velocity is equal to the speed of sound in air, shown as the horizontal red line in Fig. 3. With finite water depth H (black line in Fig. 3), we see that only flexural waves are coupled to the air wave.
Gravity waves have been observed in landfast sea ice (Sutherland and Rabault, 2016) and may theoretically couple to the airwave for the case of an infinite fluid (blue dashed line in Fig. 3), but the coupled wave is likely unobservable in practice. The 235 air-coupled flexural waves arrive in advance of the air wave because the group velocity, ⁄ , is larger than the phase velocity for a given frequency (magenta dashed line in Fig. 3). The air-coupled flexural frequency is insensitive to water depths >3.2 m, in this example, but is significantly affected by the plate acceleration term. Including the effect of plate acceleration increases the air-coupled flexural frequency by ~23% using physical properties relevant to our field data (see Table 1). While studies of moving vehicles on floating plates may safely ignore the effect of plate acceleration, we show that is important to 240 include it when the considered load is moving at the speed of sound in air.

Methodology
A key attraction of air-coupled flexural waves is the ability to use the timeseries recorded by a single sensor to estimate ice thickness. This is possible because, as illustrated by our dynamical model, the phase velocity of the air-coupled flexural wave is equal to the speed of sound in air. This is similar to studies of ice roads and runways, where the fact that the speed of the 255 vehicle is known raises the potential for routine monitoring of ice flexural rigidity using a single sensor (Squire et al., 1988).
The frequency of the air-coupled flexural wave is estimated from the high amplitude, monochromatic segment of the timeseries that directly precedes the air-wave arrival.

Air wave arrival and velocity estimation
Since the arrival time of the air wave increases linearly with horizontal offset between the seismic source and receiver, we can 260 estimate its arrival time and velocity using the linear Radon transform, also referred to as the slant-stack or − transform (Yilmaz, 2001). When the independent variable is source to receiver offset, energy with a constant velocity follows a linear Formatted: Font color: Auto (Yilmaz, 2001). When the independent variable is source to receiver offset, energy with a constant velocity follows a linear trajectory. Under the linear Radon transform, energy corresponding to a specific velocity and origin time is stacked along this trajectory to form a single point in the transform space. For a given explosive charge (see Figure 4a), we compute the linear 265 Radon transform for offset sorted geophone signals and estimate the velocity (cair) and intercept (tint) that correspond to the air wave by picking the local maximum of the transform magnitude (see Figure 4b). The estimated arrival time of the air wave at a given geophone, , is then given by, = + , where d is the source to receiver offset.

Thomson multitaper estimation of air-coupled flexural frequency
The frequency of the air-coupled flexural wave is estimated from segments of the recorded timeseries which precede the arrival 270 of the air-wave (see Figure 4c). We take the derivative of the timeseries to suppress low frequencies and linearly increase the gain of high frequencies and apply a Hamming window taper. For each timeseries, the Thomson multitaper power spectral density (Thomson, 1982) is then estimated using a time-half bandwidth product of two and a zero-padded Fourier transform length of 4096. Thomson's estimator utilizes a set of orthonormal data tapers called discrete prolate spheroidal sequences to control spectral leakage and stabilize the estimate through an inherent variance reduction (Hanssen, 1997). The air-coupled 275 flexural frequency (see Figure 4d) is then estimated from the multitaper power spectral density maximum. We only estimate the air-coupled flexural frequency for source to receiver offsets > 125 m because higher-velocity sub-seabed reflections arrive around the same time as the air wave and create too much interference at nearer offsets (see Figure 5).

Closed-form estimator of ice thickness
Once we have estimated the frequency of the air-coupled flexural wave ( ) we can estimate the ice thickness using the 285 dispersion relation (Eq. (11)) rearranged as, where = ⁄ is the wavenumber corresponding to the air-coupled flexural wave. The inclusion of the plate acceleration term means that we are left with a non-linear relation for the plate thickness, h, in terms of the plate elastic properties and the constant flexural frequency and wavenumber. We may either estimate the ice thickness by nonlinear numerical optimisation 290 of Eq. (13), using e.g. the bisection method, or we may rearrange Eq. (13) to give the cubic polynomial form, which has one real, positive root corresponding to the plate thickness h. Since Eq. (14) has a standard cubic form, a closedform analytical solution can be derived. We followed the modified Cardan's solution method outlined by Nickalls (1993) (15) all give identical results. 300

Results
Air-coupled flexural waves appear to be a robust part of the wavefield excited by explosive sources on floating sea-ice and were observed for all four field campaigns, in spite of thin, variable snow covers. Furthermore, Figure 5 illustrates that since the air-coupled flexural wave precedes the arrival of the air wave, it is recorded equally well for both point and line sources. This is significant because much of our field data employed line sources in order to attenuate the trailing low-frequency ice 305 flexural waves that are an unwanted noise component from the perspective of reflection seismic surveying (Johansen et al., 2019).

Solution of the full dynamical model
It is only strictly necessary to consider the solution to the dispersion relation in order to estimate ice properties from flexural wave observations. However, it is also beneficial to consider the solution to the dynamical equation (Eq. (1)) 315 in order to demonstrate that the key physics involved in the propagation of air-coupled flexural waves are represented by the model. To this end, we evaluated Eq. (9) numerically using an Inverse Fast Fourier Transform (IFFT) with N=2 19 discrete points, a temporal sampling interval, = 2.5 × 10 −4 , and a numerical attenuation/damping parameter, = 25. Figure 6 shows a portion of the modelled radially symmetric, expanding wavefield produced by a propagating ring delta loading that we use as an approximation of the air-wave expanding outwards from a point explosive charge and pushing 320 downwards on the ice sheet. We see that the ice sheet is deflected downwards underneath and behind the air wave. The negative deflection decays smoothly and extends to a position around 300 m behind the load, although this is outside the figure view.
Ahead of the air wave, the ice surface is elevated and we see a high amplitude, constant frequency wave-train that corresponds to the air-coupled flexural wave, which decays gradually with distance ahead of the load due to the inclusion of numerical damping in the model. 325 To compare the solutions of the full dynamical model result with field data, we calculate the velocity response of the plate at a given offset by numerical time differentiation of the output from Eq. (9), and compare that directly to th e recorded geophone 330 timeseries at the same offset. We find that local velocity estimates calculated from our solutions to the dynamical model are very similar to the real waveforms recorded by geophones, as illustrated by the example in Figure 7. This increases our confidence that we understand the underlying physics that describe the excitation of air-coupled flexural waves in floating ice sheets.

Ice thickness estimation from air-coupled flexural waves
As summarised in Table 1, ice thicknesses calculated using Eq. (15) from the estimated air-coupled flexural frequencies and 340 air-wave velocities show very good agreement with ice thicknesses measured in boreholes. Notably, this good agreement was achieved across all four field seasons with a single set of elastic properties, that are realistic for sea ice. This is partly because ice thickness is an extrinsic property that has a dominant effect on the air-coupled flexural frequency (stemming from the fact that ∝ ℎ 3 ), i.e., the 20-80 cm thickness range from this study corresponds to 400% variation in air-coupled flexural frequency. By contrast, variation of intrinsic material properties, i.e., Young's modulus (1.7-5.7 GPa) and Poisson's ratio 345 (0.33-0.39) as reported by Timco and Weeks (2010) for first year sea ice would give ~60% and ~2% variation in air-coupled flexural frequency, respectively. The combined effects of water and ice density variation over the ranges given in Table 1 would only produce a variation of ~0.3% in air-coupled flexural frequency, while water depth has no effect for depths >3.2 m.
Temperature variation over the range -20 to +2 °C would cause the speed of sound in air to vary from 318-332 m/s, resulting in ~8% variation in air-coupled flexural frequency. The apparent speed of sound in air at the measurement location is also 350 affected by horizontal wind velocity and gusts may lead to variation of +/-5 m/s over timescales of minutes (e.g. Franke and Swenson Jr, 1989). However, as described in Section 4.1, we estimate the apparent velocity of the air wave for each shot, which should ensure that our air coupled flexural thickness estimates are independent of changes in wind and temperature over time.

355
The Yong's modulus of 2.5 ± 0.2 GPa given in Table 1 was constrained by observing that the median air-coupled flexural wave thickness estimate for the 2013 field season would fall outside the range of borehole measured thicknesses for Young's modulus outside this range. The fact that this Young's modulus gives thickness estimates that fit with borehole measurements across all four field seasons indicates that the bulk, effective elastic properties of first-year sea ice in Van Mijenfjorden are relatively consistent from year to year. A somewhat higher Young's modulus of ~4 GPa was reported for Van Mijenfjorden 360 by Moreau et al. (2020a), but the deviation from the present study is likely attributable to the protected physical setting of their study site at the estuarine Vallunden Lake. It is also important to highlight that all elastic properties we assume represent bulk, apparent values corresponding to the thin isotropic plate assumption, while real floating ice sheets are a composite sandwich material consisting of multiple ice layers of varying strength (Timco and Weeks, 2010). season to include measurement of precipitation. We assume the air-coupled flexural frequency recorded at a given geophone is related to the ice thickness at the location of the geophone, as this gave the strongest correlation with borehole thickness measurements.

field season
The ice thickness estimates for the 2013 field season are summarised in Figure 8 and show good agreement with borehole measurements. There is somewhat greater spread in the results on the Y line which is likely related to the fact that the majo rity of the source positions were inline with the X line and therefore oblique to the Y line (see Figure 1). The oblique geometry reduces the precision of using the linear Radon transform to estimate the arrival time and velocity of the air wave, which 385 becomes smeared in the transform domain. Air temperatures remained consistently <0 °C and the thickness estimates are quite similar over the ~3 day field campaign, though there is some indication of ice growth with time.

field season
The results for the 2016 field season are illustrated in Figure 9. Some outliers are present in the flexural wave thickness estimates which are caused by interference of other spectral components for the extracted timeseries segments. This interference can be due to, e.g., other wavemodes such as sub-seabed reflections, snow scooter engine noises, or instrument noise due to wind. Despite these outliers, the median thickness estimates agree well with the range of drilled thicknesses (see 400   Table 1), although the along profile locations of the boreholes was unfortunately not recorded for this field campaign. with the wavelength of flexure. An alternative hypothesis that the air-coupled flexural wave is controlled by the path averaged ice properties between the source and geophone was tested by plotting the thickness estimates at the midpoint between source and receiver. While we do not include the figure for brevity, the result is a flat profile with estimates ranging from ~45-75 cm and a similar median of ~57 cm. This does not reflect the real thickness variation recorded by the borehole measurements, leading us to prefer the hypothesis that the frequency of the air-coupled flexural wave is controlled by the ice thickness and 415 flexural rigidity in the vicinity of the geophone. Indeed, we found the borehole measured spatial variation in ice thickness could only be approximated by assuming the air-coupled flexural wave estimates lie at least 98% towards the geophone (along the straight line from source to geophone).
We also observe that the air-coupled flexural wave thickness estimates significantly overestimate the ice thickness drilled in 420 the two boreholes closest land. While further field data would be needed to investigate this phenomenon in detail, it may be attributable to the effect of the land on the lateral boundary condition of the ice sheet (see Appendix 1). The infinite ice sheet modelled in the present study does not include a lateral boundary, though it may be possible to heuristically approximate this effect by defining a spatially variable, effective elastic modulus that increases near land.

field season
The air-coupled flexural wave thickness estimates for the 2018 field season are summarised in Figure 11. The thickness estimates agree well with the range of ice thickness measured in boreholes (see Table 1), though the along profile locations of the boreholes were unfortunately not recorded for this field campaign. The most striking feature we observe for the 2018 campaign, is an increase in ice-thickness of up to ~10-15 cm over the course of the campaign (a period of ~4 days). This rapid 435 increase in thickness is attributable to the weather event that occurred 26-27 th February, where >0 °C temperatures and rain led to a significant accumulation of fresh water on top of the sea ice. This event was significant enough that fieldwork was not possible on the 27 th February. The rapid decrease in air temperature on the 28 th then promptly caused the accumulated fresh surface water to freeze, a process that occurs much more quickly than basal ice accumulation due to heat loss through the overlying ice layer, salt rejection and eventual freezing of sea water. 440 (1) They are produced by the downward pressure of the air wave moving over the ice 460 surface, i.e., the moving loads on floating plates interpretation as elaborated in section 3. (2) A plate excited by a broadband pulse propagates flexural wave energy at a range of frequencies, but sound transmission from the plate is highest at the coincidence frequency (Bhattacharya et al., 1971) producing a resonance effect.
(3) Alternatively, the air-coupled flexural wave can be interpreted as a class of leaky Lamb waves where flexural wave energy radiates into the air, regulated by the phase matching condition 465 and occurring most efficiently at the grazing angle, when the phase velocity is equal to the speed of sound in air (Brower et al., 1979;Kiefer et al., 2019;Mozhaev and Weihnacht, 2002). The observation of air-coupled flexural waves and absence of dispersive ice flexural waves on the far side of open water leads by Hunkins (1960) is consistent with the interpretation that, for explosive sources on or above the ice, air coupled flexural waves are generated continuously by the downward pressure exerted by the air wave. The correspondence between experimental data and the solution of the full dynamical model (see 470 Figure 6) gives further support to this interpretation.
Our results indicate that air-coupled flexural waves are affected by spatial (see Figure 10) and temporal ice thickness variation (see Figure 11). A high degree of spatial variability in ice thickness has been reported for first year sea ice in Van Mijenfjorden (Sandven et al., 2010) and in the Artic in general (Wadhams et al., 2006). Spatial resolution is likely proportional to the air-475 coupled flexural wavelength, though further field data including high-resolution ground truth measurements would be required to elaborate this proportionality in detail. It is also important to highlight that the air-coupled flexural frequency is affected by the elastic properties of the ice, it's flexural rigidity, not just it's thickness. By assuming a constant set of elastic parameters our thickness estimates may be considered "effective thicknesses" corresponding to the assumed elastic parameters and subject to the thin isotropic plate assumption. However, we also expect that the "effective thickness" for assumed elastic properties is 480 still highly relevant when assessing the load bearing capacity of floating ice, even if it were to deviate from the true thickness.
A small percentage of outliers are present in the thickness estimates, that we attribute to spectral contamination caused by the simultaneous arrival of other wave modes with the air-coupled flexural waves. Such interference was minimised by discarding the nearest offsets where high-velocity sub-seabed reflections arrive at the same time as the air wave (see Figure 5). However, 485 other noise sources, such as snow scooter traffic, are difficult to avoid in real-world data. In the rare cases that traffic noise is recorded at the same instant as the air-coupled flexural waves, the resulting spectral interference can prevent accurate estimation of the air-coupled flexural frequency. Snow cover also likely decreases the efficiency of coupling between ice and air, though this did not appear to play a major role for the four field seasons investigated, where gimballed z-component geophones were deployed directly on top of the thin snow covers that were present. 490

Air-coupled flexural wave data acquisition
In this study, we present data acquired with linear arrays of geophones and active seismic sources that were primarily acquired to test the feasibility of seismic reflection surveying on floating sea ice. It is important to emphasize that it should be possible to employ much simpler equipment to record air-coupled flexural waves.
Indeed, from our results we would expect that a simple microphone, sensitive in the relevant frequency range and located in 495 the vicinity of the desired measurement, either above the ice-sheet or along the shoreline may be sufficient. In addition to low equipment cost, the potential to monitor ice thickness using a microphone positioned on land could be beneficial from the perspective of the potential to monitor ice thickness using a microphone positioned on land could be beneficial from the perspective of environmental monitoring, particularly early in the freezing season when the ice is too thin 500 to be safely traversed.
The use of microphones to record air-coupled flexural waves can be illustrated by a phenomenon that is familiar to the iceskating community. On thin, floating ice, skates striking or cracking the ice can excite sonorous tones that vary in pitch according to the ice thickness (Lundmark, 2001). These tones correspond to air-coupled flexural waves, simply shifted to 505 higher frequencies because of the thinner, stiffer fresh-water ice. We use the audio track from the well-known National Geographic short film "How Skating on Thin Ice Creates Laser-Like Sounds" (Rankin, 2018) to illustrate this. The frequency of the strong monochromatic component corresponding to the air-coupled flexural wave is ~725 Hz for the scene from 0:16 to 0:33, as shown in Figure 12. Using Eq. (15) we calculate that the ice was ~4.3 cm thick, assuming Young's modulus of 8.5 GPa, Poisson's ratio of 0.33, water density of 1000 kg/m 3 , ice density of 917 kg/m 3 , water depth of at least 0.3 m and sound 510 speed of 329 m/s corresponding to air temperature of -3.9 °C. This is in line with the measured ice thickness presented in the film and illustrates the potential to use air-coupled flexural waves recorded by a simple microphone to estimate ice thickness. Detonating cord seismic sources were used for the experiments presented in this study, but other impulsive sources 520 acting on a floating ice-sheet also have the potential to excite air-coupled flexural waves. Hammer blows are an alternative active source, while crack propagation is a potential passive source that produces a sudden energy release capable of exciting propagating elastic waves in the plate and acoustic pressure waves in air.
Ice quakes produced by the sudden cracking of floating ice are a well-known phenomenon (Kavanaugh et al., 2019;Moreau Ice quakes 525 produced by the sudden cracking of floating ice are a well-known phenomenon (Kavanaugh et al., 2019;Moreau et al., 2020b;Olinger et al., 2019;Ruzhich et al., 2009) and are typically related to the build-up of stress by thermal, tidal and wind forces . This phenomenon has also been studied from the structural non-destructive testing perspective, e.g., by Haider and Giurgiutiu . This phenomenon has also been studied from the structural non-destructive testing perspective, e.g., by Haider and Giurgiutiu 530 (2018) who model the acoustic emission of axis symmetric circular crested Lamb waves excited by crack propagation in a steel plate.
In order to illustrate that ice quakes are capable of producing air-coupled flexural waves, we went out and recorded a series of icequakes around midday, Jan-03 2021 at Storvatnet, Kvaløya (Norway) using the built-in mic of a Sony a6500 digital camera. 535 This is a relatively low-quality microphone with no shielding from wind noise, but the monochromatic air-coupled flexural waves can still be clearly discerned (in addition to the broadband impulses of the ice crack ruptures). The ice was not drill ed, but its thickness was estimated to be 10-20 cm by visual assessment of vertical cracks running through the clear, black ice.
The air-coupled flexural frequency of ~195 Hz indicates the thickness was 16 cm, using the same physical properties for lake ice as introduced earlier. The midday occurrence of these ice-quakes, when air temperature was relatively high, suggests they 540 were caused by thermal expansion stresses (Ruzhich et al., 2009). The results of this study and the examples shown in Figure   12 and Figure 13 suggest that the non-contact, passive recording of air-coupled flexural waves using a single microphone may provide an additional, alternative method of passive flexural wave ice-thickness estimation. This complements previously demonstrated passive seismic monitoring methods employing array based wavefield transform approaches (Moreau et al., 2020a) or noise interferometry and Bayesian inversion (Moreau et al., 2020b). The use of inexpensive single, non-contact 545 sensors is an attractive aspect of the air-coupled flexural wave method, particularly that a microphone could be placed on land when the ice is too thin to safely traverse. However, further study is needed to reveal the extent to which the temporal occurrence of ice-quakes producing air-coupled flexural waves and reduced air-ice coupling due to thick snow covers may limit the practical acquisition season.

Conclusion
Air-coupled flexural waves were found to be a robust feature of first-year sea ice excited by 560 explosive sources over a range of ice thicknesses up to ~80 cm. The physics of these waves can be understood from a theoretical perspective that has largely developed through the study of moving vehicles on floating plates. The dynamical model that we favour is straightforward to understand in terms of linear filter theory, and we derived a closed-form solution of the dispersion relation that relates ice thickness directly to the air-coupled flexural wave frequency and air-wave velocity.
We tested the proposed closed-form ice thickness estimator extensively on field data from four field seasons, and found a 565 remarkably good agreement with in-situ borehole measurements. The phenomenon of air-coupled flexural waves is relatively familiar to the wild ice-skating community, where thin ice produces air-coupled flexural waves that are readily audible . Here we build on the pioneering efforts of Press and Ewing (1951), developing a more accessible theory supported by modern numerical and Fourier transform methods, to show that the same phenomenon also occurs for much thicker sea-ice, simply shifted to a lower frequency regime. We have mainly presented data that was acquired for 570 the primary purpose of reflection seismic profiling, but a key benefit of air-coupled flexural waves is that they can also be recorded with very simple equipment, such as a single microphone located Formatted: Superscript either above the ice-sheet or along the shoreline. Ice quakes produced by natural cracking excite air-coupled flexural waves for freshwater lake ice, which raises the potential that passive recording of air-coupled flexural waves may also be possible for sea ice. 575

Appendix 1 -Possible effect of land on ice sheet lateral boundary condition
Two outlier points were observed for the 2017 field season, where ice with drilled thickness of 30-35 cm near land had aircoupled flexural frequency similar to ice of 50 cm drilled thickness further seaward (see Figure 10). To illustrate that this may be attributable to a lateral boundary condition effect, we approximate an ice sheet of finite extent as a collection of arbitrarily small, uniformly loaded, circular plates. The ice along the shoreline is represented by assuming a clamped boundary condition, 580 while the ice far from shore is considered simply supported by the neighbouring plate elements. We may then find the thicknesses of a simply supported and clamped plate that lead to equal maximum tangential stresses. From Hearn (1997), the maximum tangential stress, occurring at the center of the simply supported plate, , and the clamped plate, , is given by = 3 2 8ℎ 2 (3 + ) and = 3 2 8ℎ 2 (1 + ).
Here, is the uniformly distributed load, is the plate radius, is Poisson's ratio, ℎ and ℎ are the thicknesses of the simply 585 supported and clamped plates, respectively. Equating the two expressions gives ℎ 2 ℎ 2 = (1 + ) (3 + ) so that ℎ ≈ 0.63ℎ for = 0.33. Therefore, it is anticipated that ice 50 cm thick far from shore will experience the same maximum tangential stress under load as ice that is 0.63 × 50 = 32 cm thick at the shoreline (consistent with the outlier points from the 2017 field season). We may similarly equate the maximum deflections of the plates, which also occur at the plate 590 centers. Hearn (1997) gives the relevant expressions for the maximum deflection of the plates as where E is the Young's modulus. By equating the two deflections and simplifying we find , which again gives ℎ ≈ 0.63ℎ for = 0.33. Since both tangential stress and strain are equal for a clamped plate with 63% 595 of the thickness of a corresponding simply supported plate, the boundary condition effect can also be understood as a change in effective elastic modulus. This supports the heuristic approximation of the lateral boundary condition effect by a spatially variable, effective elastic modulus. While this theoretical description is simplistic since, e.g., the fluid-loading is ignored, it highlights that the boundary condition of a finite plate can play a significant role. It would be beneficial to explore this effect in greater detail in future studies, particularly with regard to the length scale involved in the transition from clamped to simply 600