Laboratory study of wave-induced flexural motion of ice floes

Wave-ice interactions involve complex physical processes. Well-designed laboratory investigations are indis- pensable for studying these processes. In the present study, laboratory experimental data on saline ice obtained during the HYDRALAB + project: Loads on Structure and Waves in Ice (LS-WICE) are analyzed. Here, we devise a cross-validation method that reduces the uncertainty in estimating the wavelength of surface gravity waves from wave elevation measurements made by closely located and equidistant sensors. Both experimental and numerical case studies show that the new method produces accurate results (normalized error smaller than 5%). Furthermore, experimental case studies show that the elasticity of ice lengthens the waves within the ice cover compared with the open-water wavelength, and predictions of the wavelength from the elastic-plate model are concordant with the experimental values for waves beneath intact ice sheets. In addition, we apply multivariate analysis methods to identify flexural modes of ice floes under wave actions. The analysis results suggest that multiple flexural modes exist in the motion. The results produced by using the Morlet wavelet and Prony ’ s method confirm the presence of second-order harmonics in the motion. Part of the nonlinearity likely originates from ice-ice interactions in addition to some contributions from nonlinearity in the waves.


Introduction
Wave-ice interactions have gained much attention in recent decades. This increased interest is attributed to the increase in commercial activities in the Arctic (Melia et al., 2016) and ongoing climate change in the Arctic (Stopa et al., 2016;Li et al., 2019b) and Antarctic (Alberello et al., 2019b;Vichi et al., 2019). Multiple physical processes are involved in wave-ice interactions. On the one hand, waves weaken ice by inducing ice-ice collisions (Li et al., 2020a), transporting warm water to ice-covered regions and bending ice. In particular, waves fracture ice (Collins et al., 2015;Dolatshah et al., 2018;Kohout et al., 2014) so that the exposed areas of ice increase, thereby accelerating melting Squire, 2020). On the other hand, waves contribute to the growth of pancake ice by sweeping frazil ice together (Shen et al., 2001;Smith and Thomson, 2019). From the perspective of waves, ice alters the wave dispersion relation Meylan et al., 2018), scatters waves Bennetts and Williams, 2015), and dissipates wave energy (Nelli et al., 2017;Voermans et al., 2019;Toffoli et al., 2015;Smith and Thomson, 2019).
To study these complex wave-ice interactions, controlled laboratory experiments are needed to isolate physical mechanisms of interest from other physical processes that exist in nature. Moreover, we need to apply proper data analysis techniques in order to gain an in-depth understanding of the measured processes before reaching solid conclusions. In wave-ice interaction studies, the most important wave parameters are wave amplitude and wavelength. The former is vital for studying wave attenuation (Nelli et al., 2017;Toffoli et al., 2015), while the latter is essential in examining the wave dispersion relation (Mosig et al., 2015;Sree et al., 2018;Shen, 2019). Nevertheless, estimating wavelength accurately is quite challenging and often associated with considerable uncertainty arising from small separations between the wave measurement sensors  and the evanescent modes induced by waves at the edges of ice floes Squire, 1990, 1991). To reduce the uncertainty of wavelength estimates by using wave elevation measurements from closely located and equidistant sensors, we propose a cross-validation method that comprises two pairs of methods: (1) the genetic algorithm (Kutz, 2013) and Prony's method (Hu et al., 2013), and (2) cross-spectrum analysis  and intersite phase clustering (Cohen, 2015). To the best of our knowledge, this is the first time that both Prony's method and intersite phase clustering have been applied in a wave-ice interaction study. In cross-validation, we take the mean of the least different wavelength estimates produced by the methods in each method pair. The final result obtained by the cross-validation method is the average of the estimates from the two method pairs. This proposed method inherently mitigates noise by least-squares fitting and singular value truncation, and removes outliers by discarding maximum and minimum values in the set of estimated wavelength. Both experimental and numerical examples are presented to demonstrate the high-quality estimates of wavelengths obtained by using our proposed method.
The experimental case studies presented in this paper also illustrate the effect of the elasticity of ice on modifying the wavelength of waves when ice cover is present and the performance of the elastic-plate model (Wadhams, 1981(Wadhams, , 1986 to describe the dispersion relation of waves that pass under intact ice sheets.
Another focus of this paper is the bending motion of ice floes under wave actions. Multivariate data analysis methods such as proper orthogonal decomposition (POD, Feeny and Kappagantu (1998); ) and smooth orthogonal decomposition (SOD, Chelidze and Zhou (2006); Gedikli et al. (2020)) are applied to identify the flexural modes of ice floes. We apply these methods to analyze an extended experimental dataset compared to the one investigated by Li et al. (2019a). These two methods are also employed to reconstruct the linear response of the ice floe to satisfy the linear assumptions of some theoretical models for studying the hydroelastic response of floating mats (e.g. Montiel et al., 2013aMontiel et al., , 2013b. In comparison with the data-driven approach used here for studying the bending motions of saline ice floes, Meylan et al. (2015) used theoretical (potential flow theory in combination with thin plate theory) and experimental model to study the flexural motions of two types of plastic model ice floes. Meylan et al. (2015) found that the response of ice floes to linear wave forcing is virtually linear, despite the presence of the observed nonlinear phenomenon involved in wave-floe interaction, such as overwash and slamming. In the present study, the nonlinearity in the flexural motion of ice floes is also investigated, and the source of the nonlinearity is systematically examined. This paper is organized as follows. We start with the description of the experiments in Section 2. Thereafter, we describe the proposed method to estimate wavelengths by using measurements from closely located and regularly spaced sensors and the methods to quantify wavelengths using data collected from sensors with varying separations. In Section 4, POD and SOD are presented. Section 5 provides the results produced by using the new method. Results obtained by applying POD, SOD and their discussions are presented in Section 6. The last section (Section 7) summarizes our findings.

Experiments
The experimental data presented in this study were collected during the Loads on Structure and Waves in Ice (LS-WICE) project (Tsarau, 2017). The experimental part of this project was performed in the Large Ice Model Basin (LIMB) at Hamburg Ship Model Basin (HSVA) in Germany. The whole experimental campaign comprises in total six test series. In this study, we present the analysis of the data from test series #1000 and one subset of test series #3000, i.e., test group #3100, in which the displacements of multiple points on ice floes are recorded and are suitable for studying the bending motion of ice floes.
The setups of test series #1000 and test group #3100 are shown in Fig. 1. The ice thank is 72 m long, 10 m wide and 2.5 m deep. On the lefthand end of the ice tank, a flap-type wavemaker occupies the width of Fig. 1. Experimental test setup for test series #1000 (intact ice sheet) and test group #3100 (fragmented ice cover). (a) and (c) represent test series #1000, (b) and (d) represent test group #3100 (adapted from Li et al. (2019a)). The leading ice edge, locations of ultrasound and pressure sensors in test series #1000 are shown in (c). The floe of interest in test group #3100 is highlighted in magenta and is labeled with #A in (b) and (d). The markers on floe #A are numbered consecutively as #3, #4, ⋯, #8 along the wave propagation direction in (d).
the tank (Figs. 1a,b). On the right-hand end of the ice tank, a wooden parabolic beach is installed to dissipate the incoming waves. Adjacent to the wavemaker is an 18-m-long open-water region used for generated waves to mature.
The ice cover in test series #1000 is an intact ice sheet that is approximately 49 m long and 9.8 m wide (Fig. 1a). Figs. 1a,c illustrate the locations of pressure sensors and ultrasound sensors in test series #1000 along the ice tank. Ultrasound sensors (S1 and S2) and pressure sensors (P3 ~ P9) that measure wave motions were placed at varying intervals along the length of the tank for the purpose of quantifying the wavelength and wave attenuation in ice cover accurately. In the tests, the combination of Qualisys cameras and Qualisys markers tracked the three-dimensional positions of the points on the ice. Five Qualisys markers with longitudinal distances 1.5 m apart were attached to the ice in test series #1000 (Fig. 1a).
In test group #3100, a fragmented ice cover, which is composed of multiple rectangular ice floes (3 m long and 1.63 m wide), was produced by cutting an intact ice sheet that is different from the ice sheet in test series #1000 (Figs. 1b,d). The positions of pressure sensors (P2 ~ P9) are displayed in Fig. 1b. There were 12 markers that were attached to the ice surface 0.5 m apart to monitor the positions of the points on several ice floes. The floe and markers of interest in test group #3100 are highlighted in magenta in Figs. 1b,d and labeled as #A (ice floe) and numbered between #3 ~ #8 (markers) (Fig. 1d). An overview of the experimental setup for test group #3100 can be found in Appendix G, Experimental parameters in test series #1000 and test group #3100 are listed in Table 1. Prescribed wave periods in test series #1000 are 1.6 s and 2 s, whereas the target wave periods in test group #3100 range from 1.1 s to 2 s. In these two sets of tests, the target-generated waves are linear (wave steepness is less than 0.05). The ratio of target incoming wavelength over ice floe length (L ow /L ice ) varies between 0.0814 and 0.126 in test series #1000, the same parameter ranges between 0.630 and 2.05 in test group #3100.
Before running the tests, the properties of the ice (including elastic modulus) and water are measured. Table 2 summarizes the measured properties of ice and water in test series #1000 (intact ice sheet) and test group #3100 (fragmented ice cover). The properties of ice and water between test series #1000 and test group #3100 are approximately the same with the exception of the elastic modulus, the value of which in the test series #1000 is less than half of that in test group #3100. For brevity, reader is referred to Evers (2017); Quality Systems Group of the 28th ITTC (2017); Li and Lubbad (2018) for the methods and procedures to measure the various properties of ice and water.
For completeness, we provide an overview of the LS-WICE project and the pertinent publications here: open water test (Tsarau, 2017), test series #1000 about wave-induced ice fracture (Herman et al., , 2018, test series #2000 and #3000 regarding the dynamics of ice under wave actions and the evolution of waves in fragmented ice cover (Cheng et al., 2017;Li and Lubbad, 2018;Herman et al., 2019;Cheng et al., 2019;Li et al., 2020a), test series #4000 concerning wave-ice-structure interactions , and test series #5000 about wave loads on structure in open water (Tsarau, 2017).

Methods to estimate wavelength
In this section, we propose a method to estimate wavelength from wave elevation measurements by cross-validating the wavelength estimates from pairs of methods. One method pair is composed of the genetic algorithm (Kutz, 2013) and Prony's method (Hu et al. (2013); this method pair is abbreviated as GPR), and the other method pair consists of intersite phase clustering (ISPC, Cohen (2015)) and cross-spectrum analysis ; this method pair is denoted by ICS). For brevity, details on these methods to estimate the wavelength are provided in Appendix A. The cross-validation approach is presented in Subsection 3.1.

Estimating the wavelength by using measurements from closely located and equidistant sensors
The procedures to quantify the wavelength based on cross-validating the wavelength estimates by applying various method pairs are illustrated in Fig. 2. Firstly, we illustrate the cross-validation approach applied on GPR (shown in the left panel of Fig. 2). The first step is forming a matrix of wavelength estimates (LM 1 ) from the elevation measurements at closely spaced and equidistant points produced by using Prony's method and the genetic algorithm. Then, LM 1 reduces to L M 1 when the columns of LM 1 with minimum and maximum values (e. g., the first and third columns) are deleted, i.e., when the outliers are rejected. Thereafter, two columns of L M 1 with the smallest normalized difference are chosen, i.e., the uncertainty is alleviated, giving matrix  Thereafter, the estimates of the wavelength using GPR (L GPR ) based on the cross-validation approach are obtained by taking the mean of the three largest elements in L M 1 . Following the same procedures of applying cross-validation approach on GPR, we employ ICS to obtain wavelength estimate (L ICS ), as demonstrated in the right panel of Fig. 2. The final step in cross-validation approach is finding the mean of L GPR and L ICS , hence L m .

Estimating the wavelength from sensors with varying separation
The wavelengths under the ice cover in test group #3100 have been quantified in Cheng et al. (2019) by analyzing the pressure measurements taken by pressure sensors P2 ~ P9. For test series #1000, Herman et al. (2018) only presented the wave amplitude obtained by converting the pressure measurements to wave elevations without showing any analysis results for the wavelengths by using pressure measurements. Analogous to Cheng et al. (2019), we applied the MATLAB function xcorr to estimate the wavelengths by using all available sensor pairs in test series #1000 (S1, S2 and P3 ~ P9). Wavelength estimation is prone to uncertainty due to small separations (less than half of a wavelength) between sensors . In Cheng et al. (2019), the median value of the set of estimated wavelengths for a specific test run is taken as the final result. As opposed to Cheng et al. (2019), we employ an unsupervised learning algorithm, k-means clustering (Brunton and Kutz, 2019, and references therein), to remove outliers first and then take the mean value of the set of estimated wavelengths for a specific test run as the final estimated wavelength. The number of clusters in k-means clustering is determined by the kink-finding method (Murphy, 2012), similar to the L-curve method for Tikhonov regularization (Mueller and Siltanen, 2012). The optimal means of the clusters are found by randomly initializing the values of means of the clusters multiple times, and correspond with the values that minimize the objective function involved in k-means clustering (James et al., 2013). In addition to using k-means clustering, another approach is to take the average of the wavelengths estimated that lie within the range of 0.95 ~ 1.05 of the corresponding wavelength in the open-water condition according to Herman et al. (2018) and the references therein. Hereafter, this approach is referred to as the wavelength criterion approach.

Methods to analyze the flexural response of ice floes under wave action
While the previous section introduces a new method to estimate wavelengths from wave elevation measurements, here we will apply two advanced multivariate analysis techniques, POD and SOD, to identify the flexural modes of ice floes induced by wave forcing. We find that these multivariate analysis techniques successfully reveal the hidden nonlinearities in the bending response of the ice floes.

Proper orthogonal decomposition (POD)
Both POD and SOD can be formulated as an optimization problem. POD searches for an orthogonal coordinate system with a reduced number of coordinates that can maximize the expression of the variance in the data. According to Feeny and Kappagantu (1998); Ilbeigi and Chelidze (2017); ; , in POD, proper orthogonal modes (POMs) represent spatial modes, proper orthogonal coordinates (POCs) denote the temporal modes, and proper orthogonal values (POVs) represent the proxy of kinematic energy of each mode (or subspace dimension). Hereafter, the modes and subspace dimensions are used interchangeably throughout the paper.
The cumulative energy proportion of modes up to the i th mode relative to the total response can be evaluated as where λ j is the value of the j th POV and p indicates the number of total modes identified by using POD.

Smooth orthogonal decomposition (SOD)
In contrast, SOD seeks a nonorthogonal coordinate system that maximizes the ratio between the variance in the data and the variance in the first time derivative of the data in the nonorthogonal coordinate system (Chelidze and Zhou, 2006). In SOD, smooth orthogonal modes (SOMs) represent spatial modes, smooth orthogonal coordinates (SOCs) denote the temporal modes, and smooth orthogonal values (SOVs) indicate the smoothness of the corresponding SOCs (Ilbeigi and Chelidze, 2017;Gedikli et al., 2018;Li et al., 2019a).
In SOD, SOMs are linearly independent of each other (Li et al., 2019a), while POMs are orthonormal. After orthonormalization, orthonormalized SOMs are comparable to POMs (Chelidze and Zhou, 2006). SOCs are orthogonal to each other (Li et al., 2019a;. After normalization, the normalized SOCs are similar to the normalized POCs. Further details of POD and SOD are given in Appendix C.

Modal assurance criterion (MAC)
To compare the temporal modes and spatial modes given by POD and SOD, we introduce the modal assurance criterion (MAC) (Allemang, 2003;Seyed-Aghazadeh and Modarres-Sadeghi, 2016;Farooq and Feeny, 2008). The MAC provides a quantitative comparison between the experimental measurements and theoretical predictions of the flexural mode shapes and hence is complementary to the visual comparisons performed for experimental studies in other wave tanks (e.g. Montiel et al., 2013b;Meylan et al., 2015).
The MAC is defined as: where Ψ A and Ψ B are vectors (modes produced by POD and SOD) with the same dimension. The MAC number can be thought of as the square of the correlation between the modal vectors Ψ A and Ψ B . MAC = 1 when Ψ A and Ψ B are identical, and MAC = 0 when these two vectors are independent of each other. In all other cases, the MAC number can vary between zero and one.

Examine the nonlinear response of ice floes
The MAC number identifies the linear relationship between two modes/vectors. However, it does not provide any information regarding the nonlinearities in the response. Previously, Montiel et al. (2013b,a) applied a short-time fast Fourier transform (STFFT) to isolate the linear response of floes from the nonlinear response. In the linear response, ice floes oscillate at the same frequency as the incoming wave frequency. In contrast, the frequency of the nonlinear response of ice floes differs from that of the incident wave frequency. In this study, we apply Morlet wavelet (Cohen, 2014(Cohen, , 2019 as opposed to using the STFFT. The Morlet wavelet is preferable because there are fewer parameters to set and it has higher time-frequency resolution (Gramatikov and Georgiev, 1995;Zhang et al., 2003). For detailed procedures to separate the linear response from the nonlinear response, readers are referred to Montiel et al. (2013b).

Phase difference estimated by intersite phase clustering (ISPC)
Figs. 3a,b,c,d represent the vertical oscillation signals, the instantaneous phase angles, the instantaneous phase difference and the distribution of instantaneous phase difference in the polar coordinates between the motions of markers #7 and #8, respectively. Figs. 3a,b show that there seems to be a constant phase difference between the elevation measurements at markers #7 and #8, z 7 and z 8 . Fig. 3c shows that the variation in the phase difference is quasi-steady and fluctuates around 50 ∘ , which is demonstrated more evidently by the high concentration of the phase difference in Fig. 3d. Additionally, the amount of instant phase difference clustering |IPC | = 0.998 suggests the tight clustering of the phase difference (see the mathematical definition of |IPC| in Eq. A.9, Appendix A.2). In Fig. 3c, abrupt jumps are observed near the boundaries of the phase difference time history. This is attributed to that fast Fourier transform (FFT) and inverse fast Fourier transform (IFFT) involved in obtaining the analytic signals in ISPC assume periodic boundary conditions (Kutz, 2013;Bruns, 2004). Nonperiodic boundary of the measured signals induces ripple effect when transforming between time and frequency domains (Cohen, 2014). To alleviate this edge artefact, we remove a small portion (1.56%) of the phase difference time series near its boundary. (2019) where they used data measured by pressure sensors (k ps ) and the estimates of our new method where we use measurements collected from the Qualisys system (k m ) for test group #3100 (fragmented ice cover). Furthermore, Fig. 4 displays the theoretical prediction of the wavenumber made by the mass-loading model (k mass and k mass , Peters (1950); Weitz and Keller (1950)) and elastic-plate model (k plate and k plate , Wadhams (1981); Fox and Squire (1990), along with the open-water wavenumber k ow for test group #3100. k mass and k plate represent the predictions made by using all the available varying measurements of the ice properties, i.e., the ice thickness, the ice density, and the elastic modulus of ice. k mass and k plate denote the predictions obtained via the mean values of the measured ice properties.

Wavenumbers estimated for test group #3100 (fragmented ice cover)
According to the results presented in Fig. 4, the wavenumbers estimated using the GPR and ICS approaches based on the vertical oscillation measurements of the markers match well with those from Cheng To validate the k-means clustering approach, we apply k-means clustering to filter the wavelength obtained by analyzing the pressure measurements from test runs #3110 ~ #3150 (L ow /L ice ≥ 1.02, k ow a ow ≤ 0.0308) and take the mean of the cluster close to the open-water wavelength. We find that normalized relative differences are within 3.2% with respect to k ps from Cheng et al. (2019).

Wavenumbers estimated for test series #1000 (intact ice sheet)
Fig. 5 shows the wavenumber under ice cover versus the open-water steepness (k ow a ow ) for test series #1000 (intact ice sheet) and two test runs in test group #3100 (fragmented ice cover). As illustrated in Figs. 5b,c, the wavenumber estimates by using Qualisys position measurements (k m ) are close to those based on pressure measurements (k ps ).
Normalized difference of k m relative to k ps lies within 1.05%. Moreover, the wavenumber estimates (k ps and k m ) are close to the predictions from the elastic-plate model (k plate and k plate , Wadhams (1981); Fox and Squire (1990)) in test series #1100 (intact ice sheet).
To test our proposed cross-validation approach to estimate the wavelength, we create numerical signals based on the features extracted from 75 measured vertical displacements of ice floes in waves in three different ice/wave tanks. The results of using the proposed method to analyze these numerical signals show that the normalized error of estimated wavelength is within 5% relative to ground truth, see details in Appendix B.

Effect of ice cover properties on the change of wavenumber
A comparison between Fig. 5a and Fig. 5b shows that the wavelength increases (or the wavenumber decreases) in test run #3110 (fragmented ice cover) relative to test series #1000 (intact ice sheet) according to the theoretical predictions made by the elastic-plate model (k plate and k plate ) and the wavelength estimated by using wave motion measurements (wavenumbers k ps and k m ). The reason is that the elastic modulus of the ice in test group #3100 is at least 2.3 times higher than that in test series #1100, while the other properties of the ice (density, thickness) and the water depth in the ice tank are similar between test series #1100 and test group #3100, see details in Table 2. The increase in the wavelength with larger elasticity is also seen when comparing Fig. 5c and Fig. 5d. An increase in the wavelength under ice cover with the elasticity of ice is also demonstrated in Collins et al. (2017).
To validate the applicability of the physically measured elastic modulus in the elastic-plate model for test series #1000 (intact ice sheet)  and test group #3100 (fragmented ice cover), we apply the empirical formula proposed by Sakai and Hanai (2002)(see Eqs. E.1 ~ E.3) to estimate the equivalent elastic modulus, as also performed in a recent experimental study (Parra et al., 2020). The empirical formula includes a nondimensional length scale denoted by IF, which includes parameters such as ice thickness, the length and the elasticity of ice (see Eq. E.2 in Appendix E for the mathematical definition of IF). When IF is greater than 0.6, the equivalent elastic modulus equals the physical elastic modulus. Using Eq. E.2, we find that IF is 1.82 and 0.61 for test series #1000 (intact ice sheet) and test group #3100 (fragmented ice cover), respectively. Hence, the usage of the physical elastic modulus in the elastic-plate model is justified. This finding implies that the dispersion relation of waves below ice can be fully characterized by the elasticplate model in test series #1000 having a 49-m-long intact ice sheet and in test group #3100 having an ice cover composed of multiple 3-mlong rectangular ice floes. This implication is further confirmed by the good agreement between the wavenumber predicted by the elastic-plate model (when the physical elastic modulus is applied) and the experimental wavenumbers (see Figs. 4 and 5). The only exception is the experimental result of the test run #3160 (wave frequency of 0.91 Hz). This discrepancy can be attributed to the more pronounced effects of viscosity of ice and the nonlinear material properties of ice on modifying wavelength in high frequency wave regimes, which are not incorporated in the linear elastic-plate model (Sree et al., 2018. Moreover, the steady part of the signals is very short in test run #3160 (L ow /L ice = 0.630, k ow a ow = 0.0499). Therefore, test run #3160 is excluded from the analysis from now on. The experimental studies performed by Sakai and Hanai (2002) and Parra et al. (2020) also demonstrate that the elasticplate model describes the dispersion relation of waves below intact ice sheets well.
Furthermore, the waves shorten and lengthen (i.e., the wavenumbers increase and decrease) in test series #1000 with an intact ice sheet (Figs. 5b,c) and test group #3100 having a fragmented ice cover (Fig. 4 and Figs. 5a,d) relative to the open-water wavelength, respectively.
According to earlier studies (Sakai and Hanai, 2002;Sree et al., 2018Sree et al., , 2020Herman et al., 2019;Collins et al., 2017, and references therein), regarding continuous ice cover and fragmented ice cover, whether waves shorten or lengthen is dependent on wave frequency, floe size, the viscosity, flexural rigidity and the thickness of ice; i.e., the inertia of the ice providing that ice density does not vary significantly.

Crossover frequency
The crossover frequency when the wavelength predicted by the elastic-plate model starts exceeding that in the open-water case for deep water conditions is where g is gravitational acceleration; ρ i is the density of ice; h i is the thickness of ice cover; and D i is the flexural rigidity of ice, which is defined as where ν is Poisson's ratio of ice, which is 0.3 here as in Cheng et al.
(2019); Kovalev and Squire (2020). The theoretical crossover frequency using the mean values of the ice properties is predicted to be 0.656 Hz and 0.523 Hz for test series #1000 (intact ice sheet) and test group #3100 (fragmented ice cover), respectively. These theoretical crossover frequencies are close to our observed results. By applying the graphical method to our results, we numerically determine the crossover frequency to be 0.656 Hz and 0.520 Hz for test series #1000 and test group #3100, respectively. This is why wavenumbers are larger than and less than the corresponding open-water values when f w = 0.625 Hz for test series #1000 and test group #3100, respectively (see Figs. 5c,d). Fig. 4, Figs. 5a,d and Figs. 5b,c illustrate that the wavenumber predicted by the mass-loading model (k mass and k mass , Peters (1950); Weitz and Keller (1950)) always exceeds the open-water wavenumber (k ow ) and deviates increasingly more from the open-waver value as the wave frequency increases. As exhibited in Fig. 4, the deviation in the predicted value by the elastic-plate model (k plate and k plate ) from the open-wave wavelength increases with the wave frequency when the wave frequency is higher than the crossover frequency. These two trends have also been revealed in previous studies (Sakai and Hanai, 2002;Collins et al., 2017;Herman et al., 2019). For completeness, the equations to calculate the wavenumber for open-water conditions, mass-loading model and elastic-plate model are listed in Appendix D.

Identified flexural modes of ice floes
Figs. 6, 7 and 8 illustrate the relative energy contribution of decomposed mode shapes (POVs) and SOVs, the temporal mode (TM) shapes, and the corresponding frequency components, respectively.  (b) and (c) are estimated as the mean of values found by using kmeans clustering and the wavelength criterion approach. The k mass value from the mass-loading model and the k plate value from the elastic-plate model are the predicted wavenumbers using all the available measured ice properties, such as varying ρ i , E i and h i . k mass and k plate are the predictions using the mean values of the measured ice properties. k m represents the mean value of the wavenumber estimated by using the GPR and ICS approaches.

POVs and SOVs
Using POD, Fig. 6a shows that the first three subspace dimensions account for 99.98% of the total energy for ice floe #A. The cumulative energy proportion of the subspace dimensions is calculated by summing the first i POVs and dividing by the sum of all 6 POVs, which is equal to the total number of markers on ice floe #A; see Eq. 1 for the mathematical definition.
We further examine the SOVs for the first six subspace dimensions (Fig. 6b). The first two subspace dimensions have the largest and the same SOVs. In the higher subspace dimensions, the corresponding SOVs are smaller than those in the first two subspace dimensions (Fig. 6b).

POCs and SOCs (temporal modes)
When the corresponding temporal mode shapes are investigated (see Fig. 7), we observe smooth temporal modes for the first two subspace dimensions (Figs. 7a,b). A frequency domain analysis of these temporal modes (Figs. 8a,b) indicates only one frequency component involved in these first two modes. Consecutively, we notice that the third mode is smooth and has a single dominant frequency (twice the frequency of the first two mode shapes), as illustrated in Fig. 7c and Fig. 8c.
Moreover, the temporal mode shapes in Figs. 7d,e,f vary greatly when compared with those in the lower subspace dimensions. A frequency analysis of these mode shapes in correspondence with higher subspace dimensions (fourth, fifth and sixth subspace dimensions) shows more high-frequency components (Figs. 8d,e,f).

Connections of POVs and SOVs with POCs and SOCs
The same SOVs in the first two subspace dimensions indicate that the first two temporal modes oscillate with the same frequency (see Appendix C, Eq. C.8). This finding is also reflected in Figs. 7a,b and Figs. 8a, b, in which only one frequency component is observed.
The third subspace dimension has smaller SOVs than the first two subspace dimensions. This finding suggests that the corresponding temporal mode is less smooth and contains several frequency components. However, the third mode has single dominant frequency as shown  in Fig. 8c.
Following the same procedures, we find much smaller SOVs (Fig. 6b), hence much less smooth temporal modes (Figs. 7d,e,f) and more frequency components for the fourth, fifth and sixth subspace dimensions (Figs. 8d,e,f) than for the first three subspace dimensions. This simple practice illustrates that the first three subspace dimensions give smooth time histories, but higher subspace dimensions (higher than the third subspace dimension) seem to be contaminated with noise.
The energy distribution in Fig. 6a suggests that the first two modes are the most dominant modes. Although the energy contribution of the third mode is low (Fig. 6a), the temporal variation in this mode suggests that this mode should be a linear addition to the first two modes.

POMs and SOMs (spatial modes)
The spatial modes identified by POD and SOD are illustrated in Fig. 9. The first two spatial modes resemble the first flexural mode (half sinusoidal), and the third mode is analogous with the second flexural mode (sinusoidal). For traveling wave motion, two POD modes (subspace dimensions) are needed to represent one physical flexural mode (Taira et al., 2017). The same is true for SOD (Li et al., 2019a). This is supported by the fact that both the first and the second temporal modes oscillate at incoming wave frequencies (Figs. 8a,b) and the first two SOVs are identical (Fig. 6b), i.e., the two modes have the same frequency (see Eq. C.8). Consequently, the fourth mode should correspond with the second flexural mode of the response of the ice floe induced by traveling waves. However, noise contaminates this mode, as illustrated in Fig. 7d and Fig. 8d. Hence, the fourth mode cannot be physically interpreted. In contrast, higher modes (i.e., the fifth and sixth modes) do not represent a physical mode because (1) there is no dominant frequency component in these modes (see Figs. 7e,f and Figs. 8e,f); (2) their energy contribution is negligible when compared with the first four modes (within 1%, see Fig. 6a). This finding is corroborated by the Nyquist criterion, that is, when  the number of measurement points in the spatial domain p ≥ 2n md + 1, n md spatial modes can be correctly identified (Mukundan et al., 2009;Seyed-Aghazadeh and Modarres-Sadeghi, 2016). In test group #3100, six markers were placed on ice floe #A (see Fig. 1b). Hence, in principle, two spatial physical flexural modes can be correctly extracted. By using the complex resonance approach to address the hydroelastic problem, Meylan and Tomic (2012) theoretically predicted the motion modes of an elastic plate floating atop fluid. The first two flexural modes produced by Meylan and Tomic (2012) resemble to the second and third modes identified by POD and SOD. Interestingly, although the fifth and sixth modes extracted by POD and SOD are subjected to some uncertainty due to the limited number of markers placed on ice floe #A and noise effect, their mode shapes are similar to the theoretical mode shapes presented in Meylan and Tomic (2012). In addition, the two free modes, which are not identified by POD and SOD, are related to heave and pitch motions, and they are only predicted by the theoretical model in Meylan and Tomic (2012).

Multivariate analysis results for all test runs
In test series #1000, the five Qualisys markers cover only a very small proportion of the intact ice sheet (L marker /L ice = 0.122, see Fig. 1a and Table 1). Thus the modes cannot be identified solely by data-driven approach in this test set. Below, we summarize the results of all test runs only in test group #3100 (fragmented ice cover). Table 3 shows the proportion of cumulative energy in the total vertical response of ice floe #A for test group #3100 (fragmented ice cover). For example, ENE 2 denotes the cumulative energy contribution of the first two modes, ENE 3 represents the cumulative energy contribution of the first three modes, and so on. The results for test run #3130 (L ow /L ice = 1.33, k ow a ow = 0.0197) are repeated here for comparison. The proportion of the linear part of the response generally decreases with increasing wave steepness. The first four modes account for more than 99% of the response.
A comparison between the SMs and TMs identified by POD and SOD using the MAC is listed in Table 4. The consistency in the SMs and TMs from POD and SOD drops with the increase in the subspace dimension, although these SMs and TMs are similar (MAC greater than 0.8).
SOD can also be used to estimate the incoming wave period when multiple sensors are placed on one ice floe (see Eq. C.8). The normalized error of the wave period identified by SOD by using position measurements of markers relative to the incident wave period ranges from 0.05 % ~ 0.8% for test runs #3110 ~ #3150 (L ow /L ice ≥ 1.02, k ow a ow ≤ 0.0308).

Quantifying the nonlinearity
To examine the source of higher-order effects, we investigated the motions of markers placed on floe #A using the Morlet wavelet timefrequency analysis (Cohen, 2019). As an example, the time-frequency analysis results of the vertical displacements of marker #7 are shown in Fig. 10, presented similar to Montiel et al. (2013b). The amplitudes of both the fundamental harmonic and the second-order harmonic are visible (black and light gray stripes in Fig. 10a). The temporal variations in the amplitude values of these two harmonics are seen more clearly in Fig. 10b than in Fig. 10a. |z 7 (2f w , t) | / | z 7 (f w , t)| is approximately 6.4% in the steady time window, which is not negligible. The results of |z 7 (2f w , t) | / | z 7 (f w , t)| for test runs #3110 ~ #3150 (L ow /L ice ≥ 1.02, k ow a ow ≤ 0.0308) for all the markers on ice floe #A are summarized in Fig. 12.
As discussed above, the first two subspace dimensions identified by using POD and SOD represent the first physical mode. When Prony's method is used, two complex exponentials are needed to represent one physical frequency component (Hu et al., 2013). Hence, we reconstruct the linear response of ice floes by adding only the first two subspace dimensions/complex exponentials. Thereafter, we apply the Hilbert transform (Bruns, 2004) to obtain the analytic signal. Lastly, we take the mean of the oscillation magnitudes extracted from the analytic signal of each marker in the steady time window to calculate the oscillation amplitudes of each marker.
To compare the reconstructed linear responses obtained by POD, SOD and Prony's method with those from the Morlet wavelet timefrequency analysis, we use the relative difference (RD) formula, which is defined as follows: where the superscript * denotes POD, SOD or Prony's method; the superscript M represents Morlet wavelet and Z Amp, 1 represents the amplitude vector for markers #3 ~ #8 corresponding to the incoming wave frequency. Fig. 11 illustrates that the relative difference in the fundamental amplitudes obtained by POD, SOD and Prony's method with respect to those extracted by the Morlet wavelet is negligible (less than 0.5%). Fig. 12 demonstrates that the nonlinearity (the ratio between amplitudes related to 2f w and the amplitudes corresponding to f w ) becomes more evident with increased wave steepness. Generally, the Morlet wavelet and Prony's method produce similar results except for zero values. The zero values from Prony's method indicate that Prony's method fails to identify the components that have a frequency equal to twice the wave frequency. At the same time, the corresponding ratios given by the Morlet wavelet are less than 1%. This finding suggests that the estimate of the amplitude related to 2f w is prone to uncertainty in these cases.

Investigating the source of nonlinearity
According to Fig. 4 and Figs. 5a,d, the wave lengthens in test group #3100 (fragmented ice cover). Considering this fact, along with the open-water wave steepness in test group #3100 being less than 0.05 and the attenuation of waves beneath the ice cover, linear wave theory is still applicable because the wave steepness is less than 0.05 Herman et al., 2018;Alberello et al., 2019a). Therefore, we do not expect a significant nonlinear contribution from the waves to the out-ofplane response of the ice floe, which is further substantiated by the analysis results of the measurements taken for the two open-water test runs with wave steepness values greater than 0.05 (the values of the  other open-water test runs are greater than 0.065) but closest to that of test run #3160 (k ow a ow = 4.99 × 10 − 2 ) with the largest wave steepness in test group #3100. Specifically, the analysis results suggest that z Amp, 2 m /z Amp, 1 m < 1.7% and the other higher-order components are negligible in the two open-water test runs. Additionally, z Amp, 2 m /z Amp, 1 m even exceeds 2.2% when the wave steepnesses are 2.68 × 10 − 2 and 3.08 × 10 − 2 for test runs #3140 and #3150, respectively (Fig. 12). Except for the higher-order effects arising from waves generated, scattering is another possible source of nonlinearity (Montiel et al., 2013b). However, because the open-water gap between the ice floes in test group #3100 (fragmented ice cover) is indiscernible based on video recordings, the contribution of the nonlinearity from scattering induced by the ice floe edge can be reasonably assumed to be negligible. Due to the small varying separations between ice floes, interfloe interactions induced by waves are frequent (Li et al., 2020a). As a result, interfloe interactions most likely contribute to the nonlinearity except for those from incident waves. See Appendix F for details about the open-water tests.

Conclusions
In this study, we proposed a novel method to extract wavenumber information from closely spaced and equidistant sensors. Using two experimental groups as an example, we showed that the relative difference between the results obtained from measurements taken by pressure sensors with variable separations and those from position measurements of Qualisys markers placed at a fixed interval is small (within 5%). The results of the numerical examples give a normalized error that is less than 5% when compared with the ground truth. The proposed method intrinsically reduces noise and rejects outliers that may arise from evanescent modes of flexural motion of ice floes. This method also mitigates the uncertainty due to small separations between sensor pairs, as mentioned in Cheng et al. (2019) and Sakai and Hanai   Fig. 10. Morlet wavelet time-frequency analysis of the vertical oscillations of marker #7 in test run #3130 (L ow /L ice = 1.33, k ow a ow = 0.0197). (a) The time-frequency plot of the oscillation amplitude of marker #7. (b) The temporal variations in the amplitudes of the harmonics corresponding to f w (blue line) and 2f w (red dash-dot line) are extracted from the black and gray stripes in (a), respectively. The time span between the two purple dashed lines in (a) and (b) denotes the steady time window. Fig. 11. Relative difference (RD) of the extracted fundamental amplitudes from POD, SOD and Prony's method relative to those from the Morlet wavelet for test runs #3110 ~ #3150 (L ow /L ice ≥ 1.02, k ow a ow ≤ 0.0308).

Fig. 12.
Ratios between the amplitudes corresponding to 2f w and the amplitudes corresponding to f w . The subscripts M and P denote the Morlet wavelet and Prony's method, respectively. The number in round bracket represents incoming wave steepness in each test run. (2002).
The investigation of the experimental cases also shows that the elasticity of ice increases the length of waves below the ice cover, and the elastic-plate model predicts congruent wavelength results with the experimental values for waves propagating beneath intact ice sheets.
The second part of this paper addresses the flexural motion of ice floes. Two multivariate analysis techniques, proper orthogonal decomposition (POD) and smooth orthogonal decomposition (SOD), are used to identify the flexural modes of ice floes induced by waves. The results show that both POD and SOD successfully identify the dominant first and second flexural modes. These results are consistent with the Nyquist criterion, in which case only two flexural modes can be correctly recovered when six measurement points are available on one ice floe. It is noteworthy that the identified second, third, fifth and sixth spatial mode shapes by POD and SOD are similar to the theoretically predicted first, second, third and fourth flexural modes as in Meylan and Tomic (2012). In addition, POD and SOD produced comparable results for temporal and spatial mode identification. The results yielded by POD suggest that the first flexural mode accounts for more than 94% of the total vertical response of ice floes. This finding implies that the higherorder harmonics do not contribute substantially to the vertical motion of ice floes. The ratio between the amplitude of second-order harmonics and that of fundamental harmonics (less than 8%) for each marker obtained by Morlet wavelet time-frequency analysis and Prony's method confirms the insignificant contribution of the nonlinear component to the total response. In addition, SOD correctly retrieves the incoming wave period (normalized error less than 0.8%). POD, SOD and Prony's method can also be used to reconstruct the linear response of ice floes. Hence, POD and SOD are useful tools to analyze the wave-induced bending motion of ice floes.
The analysis of the open-water tests shows that the ratio between the amplitude of the second-order harmonics and that of the fundamental harmonics is less than 1.7%, whereas the corresponding analysis results for the markers give values greater than 2.2%. Moreover, the varying gaps between the ice floes are very small; hence, the higher-order effects resulting from scattering are very limited. This finding suggests that an ice-ice interaction may be the source of the remaining nonlinear response.

Data availability
The dataset used in this study is available in Haase and Tsarau (2019).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that should have appeared to influence the work reported in this paper. mental work described in this publication was supported by the European Community's Horizon 2020 Programme through a grant to the budget of the Integrated Infrastructure Initiative HYDRALAB+, contract no. 654110. The lead author would like to thank the Hamburg Ship Model Basin (HSVA), especially the ice tank crew, for their hospitality, technical and scientific support and the professional execution of the test programme in the Research Infrastructure ARCTECLAB. The corresponding author would like to express gratitude towards colleagues from the LS-WICE project, including Hayley Shen and Meleta Truax from Clarkson University, Agnieszka Herman from the University of Gdansk, Karl-Ulrich Evers from HSVA, Andrei Tsarau and Sveinung Løset from the Norwegian University of Science and Technology and Sergiy Sukhorukov from Kvaerner AS. The lead author also thanks Dr. Lucas Yiew and Dr. Fabien Montiel for kindly sharing their experimental data.

A.1. Wavelength determined by using genetic algorithm and Prony's method (GPR)
To estimate the wavelength using the genetic algorithm, one needs to fit the measured/numerical wave data to extract the phase lag between neighboring sensors. In essence, fitting is an optimization problem. Fitting the measured or numerical signals with a sinusoidal function is to find a group of parameters that minimize the following cost function (Li et al., 2019a): where a (j) , f (j) , θ G, (j) are the wave amplitude, wave frequency [Hz] and phase [radians] in correspondence with the j th group; t is the time instant [s]; superscript G represents the genetic algorithm henceforth; and subscript k denotes the k th time instant. This optimization problem is solved by applying the MATLAB built-in function ga. Prony's method is one of the methods to decompose signals, which is analogous to the Fourier method. Prony's method linearly superimposes a series of damped complex exponential components to express signals (Hu et al., 2013;Rodríguez et al., 2018). For a discrete signal, the formula takes the form: where a j , θ j , α j and f j are amplitude, initial phase [radians], grow/decay rate [s − 1 ] and frequency [Hz] of the j th component; t k is the k th time instant; the superscript PR denotes Prony's method hereafter; r is the total number of components to be used, and this parameter is judiciously determined by applying truncated singular value decomposition (TSVD, see e.g., Mueller and Siltanen (2012)). When applying Prony's method for numerical traveling wave signals, we set r = 2, which represents two complex exponentials that correspond to the most dominant frequency component in the wave signals.
In this study, the noise-rejecting algorithm as presented in Hu et al. (2013) when implementing Prony's method is applied. To the best of our knowledge, this is the first time that Prony's method has been applied in a wave-ice interaction study. The wavelength is then estimated as: where L m represents wavelength estimated by using the measurements of markers #(m + 1) and #m. x m+1 and x m are the positions of the neighboring markers #(m + 1) and #m along the wave propagation direction. Essentially, fitting the signals with a sinusoidal function by using the genetic algorithm and using truncated singular value decomposition (Mueller and Siltanen, 2012) to determine the number of complex exponentials to reconstruct the signal in Prony's method (Hu et al., 2013) inherently reduces the noise embedded in the wave motion measurements.

A.2. Wavelength determined by applying intersite phase clustering and cross-spectrum analysis (ICS)
Cross-spectrum analysis has been used in earlier wave-ice studies to determine phase lags between two signals (Wang and Shen, 2010;Rabault et al., 2019). When compared with cross-spectrum analysis, intersite phase clustering (ISPC) is not weighted by the amplitude (Mormann et al., 2000;Cohen, 2015). In earlier studies, ISPC is also referred to as phase synchronization, phase coherence and phase-locking factor/ value (Bruns, 2004;Cohen, 2014). Here, we use ISPC in wave-ice interactions. The wavelengths estimated by cross-spectrum analysis and ISPC are: where the superscripts I and CS denote intersite phase clustering and cross-spectrum analysis, respectively; Δθ m I and Δθ m CS are the estimated phase lags between the signals recorded for markers #(m + 1) and #m by ISPC and cross-spectrum analysis, respectively. Δθ m I and Δθ m CS are estimated as: where S m+1, m is the cross-spectrum of z m+1 and z m . IPC is defined as where θ m+1 k and θ m k are the instantaneous phases of the analytic signals z H m+1 and z H m for vertical oscillation measurements z m+1 and z m , respectively. These analytic signals are obtained by the Hilbert transform. For brevity, readers are referred to Bruns (2004) for more details about the Hilbert transform.
Note that the magnitude |IPC| indicates the concentration of the instantaneous phase difference (θ m+1 k − θ m k ) (see, e.g., Fig. 3d). When the instantaneous phase difference is distributed uniformly, |IPC | = 0, while |IPC | = 1 when the phase difference is a constant.

Appendix B. Testing the cross-validation method to estimate wavelength by using numerical signals
To augment the database to test our proposed cross-validation method to estimate wavelength, we supplement the experimental data from the LS-WICE project with synthetic numerical data.
Using Prony's method to analyze 75 measured signals for vertical displacements of ice floe from test series #1000, test group #3100 and laboratory experiments on wave-ice interaction from the other two wave tanks (Montiel et al., 2013b;Yiew et al., 2019), we observe several frequency components embedded in the motion signals where the fundamental harmonic is dominant and the motion signals selected in the steady time window is not completely stationary (not shown). Hence, we propose following equations to generate numerical wave signals (see also Li et al., 2020b). , where j = 1, 2, ⋯, 10; ε represents white noise (0% for a smooth signal and 5% for a noisy signal); a wi is the oscillation amplitude of the first marker; α denotes the spatial attenuation coefficient; β = 50 for a smooth signal and β = 20 for a noisy signal. The sign of α j m is determined by the sample of the uniform distribution U(0, 1), where α j m > 0 when the sample is larger than 0.5, otherwise α j m < 0. L wi is the wavelength of the waves beneath ice cover.
The values of f w , a wi , α and L wi are taken from the literature, except L wi and a wi for cases #1 ~ 5 and #7, which are assigned with assumed physical values; these values are listed in Table B.1. d is the constant longitudinal separation (along the wave propagation direction) between adjacent markers. Note that we tailor the proposed method for Qualisys markers, and in principle, the method is applicable for other closely spaced sensors.   The normalized error (NE) relative to the ground-truth value L wi is within 4% by means of the GPR method, while the ICS method yields lower NEs (less than 2.5%). The mean of NE GPR and NE ICS , i.e., NE m , is less than 2.5%.
The results for noisy signals are presented in Fig. B.3. Similar to Fig. B.2, NE ICS is smaller than NE GPR . NE m lies within 4.5%. Although the GPR approach produces less accurate estimates of wavelengthsthan the ICS approach based on Figs. B.2 and B.3, a combination of these two methods is needed because only the mean of the results from these two methods accurately quantifies the wavelengths from the experiment (Figs. 4 and 5).

C.1. Proper orthogonal decomposition (POD)
Let the measurement matrix of wave elevations be: where n is the number of points sampled in the time domain and p represents the number of measurement points in the spatial domain. The optimization problem of POD can be translated into an eigenvalue problem as (Feeny and Liang, 2003): where the columns of the eigenvector matrix Ψ ∈ ℝ p×p give proper orthogonal modes (POMs), i.e., spatial modes or basis vectors of the orthogonal coordinate system to be found; diagonal elements of the eigenvalue matrix Λ are proper orthogonal values (POVs). In other words, POD fits the best ellipsoid to the data matrix in the least squares sense, where the length of the semiprincipal axes corresponds to the standard deviation in the projected points on the principal axes and represents the square root of POVs . Each POV (λ i ) is related to the energy in the corresponding mode (or subspace dimension); hence, POD provides a way to perform an energy-optimal order reduction. For practical calculations, we employ economical singular value composition (SVD) to solve the POD problem, which is expressed as (Chatterjee, 2000;Epps and Techet, 2010): where the columns of UΣ ∈ ℝ n×p represent proper orthogonal coordinates (POCs), the columns of the orthonormal matrix U ∈ ℝ n×p denote the temporal modes; the diagonal singular value matrix Σ ∈ ℝ p×p is related to Λ by Σ 2 = Λ; the columns of orthonormal matrix Ψ ∈ ℝ p×p represent spatial modes (Feeny, 2002;Amabili et al., 2003;Ilbeigi and Chelidze, 2017;Epps and Techet, 2010).

C.2. Smooth orthogonal decomposition (SOD)
SOD as an optimization problem can be translated into a generalized eigenvalue problem (Chelidze and Zhou, 2006): where Λ is the diagonal generalized eigenvalue matrix, the columns of Ψ are the smooth projection modes (SPMs) that represent the nonorthogonal coordinate system to be found and V z is the first time derivative of Z.
For stable computations, we apply the economical generalized singular value decomposition (GSVD), which is formulated as (MATLAB, 2019;Chelidze, 2014;Gedikli et al., 2020): where Ũ and Ṽ are orthonormal matrices; C and S are nonnegative diagonal matrices. A vector of generalized eigenvalues, i.e., the vector formed by diagonal elements of Λ , is calculated by elementwise division between diag(C T C) and diag(S T S) (Chelidze, 2014).
According to Chelidze (2014), smooth orthogonal modes (SOMs) that represent spatial modes are the columns of Φ ∈ ℝ p×p . Smooth orthogonal coordinates (SOCs) that represent temporal mode shapes are the columns of Q = ZΦ − T =ŨC ∈ ℝ n×p . SPMs and SOMs are related to each other by Ψ = Φ − T . Smooth orthogonal values (SOVs), which imply the smoothness of SOCs, are given by the generalized eignevalues, viz. the diagonal entries of Λ .
Smoother SOCs correspond to larger SOVs.
The oscillation period of each mode can be estimated as (Chelidze and Zhou, 2006): where λ i is the i th diagonal element of Λ .

Appendix D. Theoretical wave dispersion relations for wave-ice interactions
Wave dispersion relations corresponding to open-water conditions, mass-loading model and elastic-plate model are as follows. The open-water wave dispersion relation is (Holthuijsen, 2010): where ω w is the angular wave frequency in radians; k ow is the open-water wavenumber; and H is the water depth.
The wave dispersion relation that corresponds to the mass-loading model (Peters, 1950;Weitz and Keller, 1950) is: where ρ is the mass density in kg/m 3 . The subscripts i and w represent ice and water, respectively; k w denotes the wavenumber of the waves passing ice cover. The wave dispersion relation for the elastic-plate model (Wadhams, 1981(Wadhams, , 1986Fox and Squire, 1990) is expressed as: where D i is the flexural rigidity of ice.

Appendix E. Empirical equations for the equivalent elastic modulus
Sakai and Hanai (2002)  where E i and E eq are the physically measured elastic modulus and the equivalent elastic modulus when viewing whole ice cover as a homogeneous continuum, respectively. I i and I c are the length and characteristic length of the discrete ice floes, respectively.

Appendix F. Open-water tests
The water depth and water density for open-water tests are 2.51 m and 1005.5 kg/m 3 , respectively. The parameters of waves for open-water tests are summarized in Table F.1.

Appendix G. Overview of the experimental setup
After all the test runs completed in test group #3100, the top-view images of the ice cover were taken by using a camera fixed on a crane. These images were stitched by means of the FIJI image processing tool (Schindelin et al., 2012) and in-house MATLAB scripts, similar to Li et al. (2020a).