Sea ice thickness 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 of the ice sheet, which is controlled by a combination of thickness and elastic properties. We develop a theoretical framework 10 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 real and measurable component of the total wavefield of floating ice-sheets excited by impulsive sources and we present a simple closed-form estimator for the ice thickness based on 15 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 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 ice. While air-coupled flexural waves 20 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 results of this study and further examples from lake ice, suggest the possibility of non-contact estimation of ice thickness using simple, inexpensive microphones located above the ice-sheet or along the shoreline. While we have demonstrated the use of aircoupled flexural waves for ice thickness monitoring using an active source acquisition scheme, naturally forming cracks in the 25 ice are also shown as a potential impulsive source that could allow passive recording of air-coupled flexural waves.


Introduction
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 30 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. 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 35 and Crary (1934) and the study of Yang and Yates (1995) that 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).

40
Despite their commendable efforts to develop a theoretical foundation for the air-coupled flexural wave (Press and Ewing, 1951), the term has not been widely used since the initial investigation period in the 1950's and the study of Hunkins (1960).
This may be because the theory developed by Press and Ewing (1951) is cumbersome, being built-up using demanding analytical integration methods at a time when computing power to handle more convenient numerical methods did not exist.
We think this is unfortunate since the term "air-coupled flexural wave" gives a concise description of the physical mechanism 45 that produces these waves. Furthermore, their monochromatic frequency is attractive from an experimental standpoint, since it can be estimated from a single timeseries more readily than the precise time-frequency evolution of the highly dispersive ice flexural wave. Air-coupled flexural waves may therefore provide a simple and effective means to estimate the flexural rigidity of a floating ice sheet and to study variation in ice thickness, for the case where the elastic properties of the ice sheet are assumed or independently estimated. 50

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 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 55 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 . 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). 60 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 65 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).
Coupling between air and guided waves in engineered structures, like concrete slabs, is also relevant to non-contact applications of non-destructive testing (e.g. Harb and Yuan, 2018;Zhu, 2008).

70
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 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 75 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 have captured significant media attention (Griffin, 2018;Rankin, 2018). This phenomenon is analogous 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 80
We now consider the closely related topic of moving vehicles travelling over floating ice sheets. The topic of moving loads on floating plates has received considerable attention due to the importance of roads on floating sea, river and lake ice that are 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 85 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 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 90 (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, 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. 95 We propose that the air-coupled flexural wave phenomenon can be understood as a special case of the more general 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 100 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 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. 105

Study area and data acquisition
In this study we investigate a series of active-source seismic experiments conducted in the innermost part of Van Mijenfjorden, 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 110 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 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 115 in-line with the source, but we also include oblique arrangements from the 2013 campaign that employed a cross type array (see Figure 1).
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 120 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 verticalcomponent 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 125 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 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 with its origin at the upper undisturbed water surface. We impose the constraint that the fluid must satisfy the Laplace equation, ∇ 2 = 0, where 140 is the velocity potential in the fluid. In addition, we impose a non-cavitation condition at the ice-water interface, | = = , where b is the draught of the ice and a normal flow condition at the sea bottom, | = = 0, where ( , , ) is the vertical deflection of the plate's neutral surface. The linear spatiotemporal dynamics of the system is described by the partial differential equation (e.g., Squire et al., 1996) D∇ 4 + ℎ (1) 145 Here, D = ℎ 3 12(1− 2 ) 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.

150
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.
(2) to the spatiotemporal fields 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 bending forces in the plate, is due to the plate buoyancy, ℎ 2 is the plate acceleration and ( 2 ) coth( ) arises due 160 to the constraint of finite water depth.
The general solution for the resulting spatiotemporal deflection wavefield, ( , , ), is then given by the full three-dimensional inverse FT of ( , , ) as .

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 using the ring delta function, 170 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 so the spatio-frequency deflection wavefield can be written as 175 where 0 ( ) = ( ) 4 + − ℎ 2 − ( ) coth ( ).
The spatiotemporal deflection wavefield in radial coordinates then becomes 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 satisfy the equation 185 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.
(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.
However, when the load moves at the speed of sound in air, plate flexural waves with much shorter wavelengths are resonant 195 and it becomes important to retain the plate acceleration term. The common approximations that are valid at typical vehicle speeds may lead to significant inaccuracies when estimating physical properties of the floating ice sheet from estimates of aircoupled flexural wave frequency.

200
Figure 3 -Comparison of simplified dispersion relation neglecting water depth and plate acceleration (blue, Eq. 12) with the full dispersion relation including these effects (black, Eq. 11) for h = 0.74 m and sea-ice physical properties from Table 1. 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 205 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 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 210 >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 include it when the considered load is moving at the speed of sound in air.

Methodology 215
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 vehicle is known raises the potential for routine monitoring of ice flexural rigidity using a single sensor .
The frequency of the air-coupled flexural wave is estimated from the high amplitude, monochromatic segment of the timeseries 220 that directly precedes the air-wave arrival. The air-coupled flexural frequency can then be directly related to the physical properties of the floating ice sheet using the dispersion relation introduced in the previous section.

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 estimate its arrival time and velocity using the linear Radon transform, also referred to as the slant-stack or − transform 225 (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 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 230 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 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 235 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 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 240 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 dispersion relation (Eq. (11)) rearranged as, 250 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 of Eq. (13), using e.g. the bisection method, or we may rearrange Eq. (13) to give the cubic polynomial form, 255 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) to derive the following estimator,

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.

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 275 wave observations. However, it is also beneficial to consider the solution to the dynamical equation (Eq. (1)) 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 280 an approximation of the air-wave expanding outwards from a point explosive charge and pushing 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. 285 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 the recorded geophone 290 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 300 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 305 (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 310 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.

315
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 320 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). Here we give a summary of our air-coupled flexural wave ice thickness estimates for each field season and examine them in the context of borehole measurements of ice thickness and prevailing meteorological conditions, as represented by records at 330 the nearby Sveagruva weather station (The Norwegian Meteorological institute, 2020), which was upgraded prior to the 2017 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 335
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 majority 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 becomes smeared in the transform domain. Air temperatures remained consistently <0 °C and the thickness estimates are quite 340 similar over the ~3 day field campaign, though there is some indication of ice growth with time.

field season 350
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 Table 1), although the along profile locations of the boreholes was unfortunately not recorded for this field campaign. 355 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 flexural rigidity in the vicinity of the geophone. Indeed, we found the borehole measured spatial variation in ice thickness 370 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 the two boreholes closest land. While further field data would be needed to investigate this phenomenon in detail, it may be 375 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 385
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 increase in thickness is attributable to the weather event that occurred 26-27 th February, where >0 °C temperatures and rain 390 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.

6 Discussion
It is possible to interpret the origin of the air-coupled flexural waves in several equivalent ways: (1) They are produced by the downward pressure of the air wave moving over the ice 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 405 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 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 410 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 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 415 (Sandven et al., 2010) and in the Artic in general (Wadhams et al., 2006). Spatial resolution is likely proportional to the aircoupled 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 420 to the thin isotropic plate assumption. However, we also expect that the "effective thickness" for assumed elastic properties is 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 425 the nearest offsets where high-velocity sub-seabed reflections arrive at the same time as the air wave (see Figure 5). However, 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 430 geophones were deployed directly on top of the thin snow covers that were present.

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 435 simple microphone, sensitive in the relevant frequency range and located in 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 environmental monitoring, particularly early in the freezing season when the ice is too thin to be safely traversed.

440
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 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 445 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 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. 450 Detonating cord seismic sources were used for the experiments presented in this study, but other impulsive sources 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 460 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 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 (2018) who model the acoustic emission of axis symmetric circular crested Lamb waves excited by crack propagation in a steel plate. 465 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. 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 drilled, 470 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 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 475 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 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 480 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 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 495 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 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 500 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 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 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 505 sea ice.

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 510 small, uniformly loaded, circular plates. The ice along the shoreline is represented by assuming a clamped boundary condition, 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 + ). 515 Here, is the uniformly distributed load, is the plate radius, is Poisson's ratio, ℎ and ℎ are the thicknesses of the simply 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 centers. Hearn (1997) gives the relevant expressions for the maximum deflection of the plates as , 525 which again gives ℎ ≈ 0.63ℎ for = 0.33. Since both tangential stress and strain are equal for a clamped plate with 63% 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 530 in greater detail in future studies, particularly with regard to the length scale involved in the transition from clamped to simply supported behavior. Accurately describing this behavior would require the collection of more detailed field data in the zone near the shoreline in order to constrain the development of a more complete model, such as a full finite element simulation.

Code and data availability
Data and code used to produce this research can be shared upon request to the authors. 535

Author contributions
Development of theory and associated modelling was carried out by RR and AH. The field campaigns were administered and led by TAJ and BOR. Initial data preparation was conducted by BOR, while the air-coupled flexural wave data processing methodology was developed and implemented by RR. RR was responsible for analysing and visualising the data and writing the manuscript with contributions from all authors 540

Competing interests
The authors declare that they have no conflict of interest.

Financial support
This research has been supported by the University of Tromsø -The Arctic University of Norway, ARCEx partners and the 545 Research Council of Norway (grant no. 228107). The publication charges for this article have been funded by a grant from the publication fund of UiT The Arctic University of Norway.