A collection of wet beam models for wave–ice interaction

Theoretical models for the prediction of decay rate and dispersion process of gravity waves traveling into an integrated ice cover expanded over a long way are introduced. The term “wet beam” is chosen to refer to these models as they are developed by incorporating water-based damping and added mass forces. Presented wet beam models differ from each other according to the rheological behavior considered for the ice cover. Two-parameter viscoelastic solid models accommodating Kelvin–Voigt (KV) and Maxwell mechanisms along with a one-parameter elastic solid model are used to describe the rheological behavior of the ice layer. Quantitative comparison between the landfast ice field data and model predictions suggests that wet beam models, adopted with both KV and Maxwell mechanisms, predict the decay rate more accurately compared to a dry beam model. Furthermore, the wet beam models, adopted with both KV and Maxwell mechanisms, are found to construct decay rates of disintegrated ice fields, though they are built for a continuous ice field. Finally, it is found that wet beam models can accurately construct decay rate curves of freshwater ice, though they cannot predict the dispersion process of waves accurately. To overcome this limitation, three-parameter solid models, termed standard linear solid (SLS) mechanisms, are suggested to be used to re-formulate the dispersion relationship of wet beam models, which were seen to construct decay rates and dispersion curves of freshwater ice with an acceptable level of accuracy. Overall, the two-parameter wet beam dispersion relationships presented in this research are observed to predict decay rates and dispersion process of waves traveling into actual ice covers, though three-parameter wet beam models were seen to reconstruct the those of freshwater ice formed in a wave flume. The wet beam models presented in this research can be implemented in spectral models on a large geophysical scale.

1 Introduction

Mutual interaction between water waves and ice is a multi-physical problem, frequently occurring in polar seas where waves can penetrate ice covers, traveling over kilometers until they die out. The phase and group speeds of the resulting gravity-flexural waves advancing through the ice can be different from that of an open-water sea owing to the effects of forces caused by solid motions. There is a pressing need to understand the mutual effects of ice and gravity waves on each other due to the recent retreat of the sea ice in the Arctic (Stroeve et al., 2008; Comiso et al., 2008; Meier et al., 2013) and the emergence of large and powerful wind-generated water waves in the Antarctic (Young et al., 2011), which can affect the ice extent (Kohout et al., 2014) by breaking the ice. Specifically, the break-up of ice may potentially lead to evolution of larger waves. Hence, it is very important to understand the pattern of wave propagation in polar seas on larger geophysical scales. Wave modelers aiming to numerically simulate wave propagation in polar regions need to use proper formulations to calculate the ice-induced energy damping and the group speed of waves traveling into ice covers (Rogers and Orzech, 2013; Liu et al., 2020). These two can be helpful in the prediction of the amount of wave energy traveling into an ice cover, which can be calculated through spectral modeling of waves on a large geophysical scale, and if they are combined with sea ice break-up models, they can also provide us with evolution of marginal ice zones (e.g., in Kousal et al., 2022).

Mathematical modeling of the wave–ice interaction firstly received the attentions of researchers in the 19th century. The first model was developed by Greenhill (1886), who formulated harmonic linear motions in a fluid domain covered with an elastic beam. To build the model, he assumed that the ice extent was spanned over an infinite way, and the solid body had relatively small motions. This model lacks energy damping, but, instead, it can be employed in the prediction of the dispersion processes, the result of which is observed to be consistent with the physics of long integrated bodies covering water.

The Greenhill study is the kernel of the next generation of models established for the prediction of mutual effects of water waves and ice. Researchers looking into the wave–ice interaction developed models by modifying the original model of Greenhill. Early developments of models dates back to years between the 1950s and 1970s when scientists had become able to reach polar seas, recording the wave climate. The energy of waves traveling through the ice was observed to be reduced by sea ice (Robin, 1963; Wadhams, 1972) and the phase speed was observed to be affected by the ice. The Greenhill model lacked the former since it was developed for an elastic solid body.

Following the Greenhill's model, various linear mathematical models of wave–ice interaction have been developed. In some of them, interactions between water waves and finite length ice floes have been modeled. Using such an assumption, scattering and radiation problems can be addressed, and motions of a flexible ice floe can be found (e.g., Meylan and Squire, 1996; Smith and Meylan, 2011; Meylan et al., 2020). Such a problem has been solved in frequency or time domains (e.g., Meylan and Sturova, 2009; Hartmann et al., 2022), and also shares similarities with engineering problems within the field of marine hydroelasticity (e.g., ship hydroelasticity and seakeeping of very large floating structures; e.g., in Hirdaris et al., 2003).

Some other models highlight mutual interaction of water waves with an infinite length ice cover. Simply stated, ice cover is viewed as a continuum medium. When this approach is embarked, the ice-induced energy decay is found through finding the root of the dispersion equation. To address the energy decay rates of continuum models, scholars have mostly prescribed viscoelastic (e.g., Squire and Allan, 1977; Wang and Shen, 2010) or proelastic (Chen et al., 2019; Xu and Guyenne, 2022) behavior for the ice, or presumed that a thin layer of viscous fluid (e.g., De Carolis and Desiderio, 2002) can be representative of the ice behavior. In addition, in some other models, a damping term, which represents a fluid-based dissipation mechanism, is considered, though the ice is assumed to be elastic. Such an approach has been introduced by Robinson and Palmer (1990), and is known as the RP model. Researchers employing this approach have prescribed elastic behavior for the ice cover (Squire et al., 2009; Williams et al., 2013). This fluid damping term has also been combined with scattering models (Williams et al., 2017). To consider the scattering of water waves by an infinite length cover, or wave reflection by cracks or variation in the thickness, matching methods have been used (e.g., Fox and Squire, 1991; Kohout and Meylan, 2006; Kohout et al., 2007; Zhao and Shen, 2013).

Models developed for the prediction of decay rates and the dispersion process of continuously integrated ice have been mostly developed by assuming small displacement for the cover. Thus, they solve the solid dynamic problem by using an Euler–Bernoulli beam theory. Their applications have been seen to be limited. This can be due to the reason that they have mostly developed by simplifying the problem, neglecting some aspects of the fluid and solid motions. Studies performed in the recent decade provide a clear picture of this fact. In some studies concerned with the dispersion process of waves advancing in ice (or elastic) covers, formulated dispersion relationships were reported to reconstruct the dispersion plot with an effective value of rigidity (or Young's modulus), which is much smaller than what was measured in dry tests (Langleben, 1962; Sakai and Hanai, 2002; Cheng et al., 2019). In some other studies, different values of ice viscosity were seen to give the best fitting for energy damping of various frequencies, which may not agree with reality (Marchenko et al., 2021).

To overcome the limitations of available continuum models, the role of different mechanisms in energy dissipation and the parameters influencing the dispersion process of waves propagating in ice should be understood well. Reviewing the structure of all developed models (example of review papers: Squire et al., 1995; Squire, 2007; Zhao et al., 2015; Collins et al., 2016; Squire, 2020), one can conclude that energy dissipation is either assumed to be triggered by solid motions (e.g., Wang and Shen, 2010) or by fluid motion, such as fluid damping, overwash or shear stresses (e.g., Liu and Mollo-Christensen, 1988; Mosig et al., 2015; Herman et al., 2019; Huang et al., 2019). Furthermore, the dispersion process is dependent on ice flexural rigidity, thickness and its density in most of the models, though added mass can also affect the dispersion process. New continuum models can still be developed by assuming that the fluid-based and solid-based energy damping mechanisms emerge at the same time. The coexistence of solid-based and fluid-based dissipation mechanisms in formulations may help us predict the decay rates with fewer limitations and solve the fluid–solid problem for a more realistic condition. This means that, unlike the previous studies which presented dispersion relationship based on the RP model by assuming that ice is an elastic material, viscoelastic behavior of the ice can also be considered. This provides us with more options in building the dispersion relationship as viscoelastic behavior of the ice cover can also be formulated using different models (e.g., Kelvin–Voigt and Maxwell), which may allow for incorporation of viscosity in flexural rigidity.

The present paper aims to develop continuum wave–ice interaction models by hypothesizing that fluid forces and solid forces emerge simultaneously. Distinguishing the developed models from the previous studies, fluid forces, including damping and added mass, are hypothesized to emerge under a viscoelastic ice beam covering water, vibrating due to the wave forces. To understand the role of rheological behavior of ice, two different two-parameter solid mechanisms, linking stress and strains arising in the solid beam, are employed to establish these models. Accordingly, two new viscoelastic models with incorporation of fluid damping force are introduced. It is also aimed to provide a better understanding of how different rheological patterns can be implemented in the prediction of decay rate of the ice-covered water and seas. This paper is structured as follows. Models are developed in Sect. 2. In Sect. 3, results including wave height decay rates and dispersion curves are presented. At the first step, sensitivities of models to different parameters, shear modulus and viscosity, are analyzed. At the next step, predictions of models are compared against those of field and flume measurements. Afterward, two three-parameter solid models are also used to formulate dispersion relationships to investigate the way consideration of more complex rheological behavior can affect the results. Finally, a discussion on the ability of models in the prediction of decay rates and dispersion curves is presented. In Sect. 4, concluding remarks and suggestions are presented. Appendix A presents an example of dimensional data.

2 Models

2.1 Development of models

Consider a two-dimensional fluid domain containing water. The domain extent is stretched over an infinite length and has a depth of D. An ice sheet covers the upper surface of the domain. It implies that no air–water interaction occurs at all. A schematic of the domain is shown in Fig. 1c. Water waves propagate in this domain and their energy is dissipated over time. Assuming that wave height is decayed exponentially, wave height at a longitudinal position of x is formulated as

$$\begin{array}{}\text{(1)}& H\left(x\right)=H_{\mathrm{0}}\left(x\right)e^{-\mathit{\alpha }\left(x-x_{\mathrm{0}}\right)}.\end{array}$$

In Eq. (1), α is the decay rate of wave height, and H₀(x) is the wave height at x=x₀. In addition, as was explained in Sect. 1, the wavelength of ice-covered sea (the distance between two consecutive wave crests) can be different from that of open-water condition (see Fig. 1b). Scattering is not considered in the problem, as the cover is assumed to have an infinite length.

Let the fluid be ideal and irrotational. Hence, the fluid motion is represented by the potential field, which is indicated with Φ(x,t). Assuming that fluid has a linear cyclic motion with a frequency of ω, the potential field can be re-written in the frequency domain, as per

$$\begin{array}{}\text{(2)}& \mathrm{\Phi }(\mathit{x};\mathrm{t})=\text{Re}\left[\mathit{\varphi }\left(\mathit{x}\right)e^{-i\mathit{\omega }t}\right].\end{array}$$

The Laplace equation holds the fluid domain:

$$\begin{array}{}\text{(3)}& \begin{array}{lll}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\varphi }\left(\mathit{x}\right)=\mathrm{0}& -\mathrm{\infty }<x<\mathrm{\infty }& -D<z<\mathit{\xi }.\end{array}\end{array}$$

Here, ξ is the elevation of the upper layer of the fluid domain with respect to the calm water line.

The vertical component of velocity is zero at the seabed, which signifies that

$$\begin{array}{}\text{(4)}& \begin{array}{lll}{\partial }_z\mathit{\varphi }\left(\mathit{x}\right)=\mathrm{0}& -\mathrm{\infty }<x<\mathrm{\infty }& z=-D.\end{array}\end{array}$$

The solid body covering the upper layer of the water is assumed to be very thin, and its thickness is denoted with h. Assuming linearity, vertical motion of the upper layer is expressed as Re[ξ(x)e^−iωt]. In the absence of fluid forces, the solid layer follows the Euler–Bernoulli beam theory, as per

$$\begin{array}{}\text{(5)}& \begin{array}{ll}-{\mathit{\omega }}^{\mathrm{2}}\left({\mathit{\rho }}_{\mathrm{i}}h\right)\mathit{\xi }+\frac{G_E{h}^{\mathrm{3}}}{\mathrm{6}\left(\mathrm{1}-\mathit{\nu }\right)}{\mathit{\xi }}_{xxxx}=p& -\mathrm{\infty }<x<\mathrm{\infty },\end{array}\end{array}$$

where G_E is the dynamic shear modulus of the material, which shall be introduced later. p is the pressure acting on the beam, causing the vibration. Also, ρ_i is the ice density. The above equation is formulated for a dry beam (see Fig. 1a).

It is assumed that fluid forces emerge per unit area when the beam interacts with the water, which are given as

$$\begin{array}{}\text{(6)}& f={\mathit{\omega }}^{\mathrm{2}}A\mathit{\xi }+i\mathit{\omega }B\mathit{\xi }.\end{array}$$

Here, A is the added mass coefficient and B is the damping coefficient. The damping term was previously used for formulation of dispersion relationship of a pure elastic beam by Robinson and Palmer (1990). These two coefficients are phenomenological. Equation (5) can be extended to

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \left(-{\mathit{\omega }}^{\mathrm{2}}\left({\mathit{\rho }}_{\mathrm{i}}h+A\right)\right)\mathit{\xi }-i\left(\mathit{\omega }B\right)\mathit{\xi }+{\displaystyle \frac{G_E{h}^{\mathrm{3}}}{\mathrm{6}\left(\mathrm{1}-\mathit{\nu }\right)}}{\mathit{\xi }}_{xxxx}=p\\ & -\mathrm{\infty }<x<\mathrm{\infty }.\end{array}\end{array}$$

When the ice cover is broken, the body–body interaction may emerge, causing an extra dissipation mechanism, which can be implemented through consideration of an artificial damping term (Lu et al., 2010). But, such a damping mechanism is not considered in the present research as the ice cover is assumed to be integrated.

Equation (7) is an extended version of the Euler–Bernoulli beam model which is adopted for a beam interacting with water. Hence, the term “wet beam” is used with an aim to distinguish the model developed based on Eq. (7) from what was presented in Eq. (5). The time derivative of surface elevation is approximately ∂_zϕ, and the time derivative of pressure under the beam is approximately $-{\mathit{\rho }}_{\mathrm{w}}g\left({\partial }_z\mathit{\varphi }\right)+{\mathit{\rho }}_{\mathrm{w}}{\mathit{\omega }}^{\mathrm{2}}\mathit{\varphi }$. Thus, the boundary condition on the fluid–solid interface is formulated as

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \left({\mathit{\rho }}_{\mathrm{w}}g-{\mathit{\omega }}^{\mathrm{2}}\left({\mathit{\rho }}_{\mathrm{i}}h+A\right)\right)\left({\partial }_{\mathrm{z}}\mathit{\varphi }\right)\\ & -i\left(\mathit{\omega }B\right)\left({\partial }_z\mathit{\varphi }\right)+{\displaystyle \frac{G_E{h}^{\mathrm{3}}}{\mathrm{6}\left(\mathrm{1}-\mathit{\nu }\right)}}{\partial }_{xxxx}\left({\partial }_{\mathrm{z}}\mathit{\varphi }\right)={\mathit{\rho }}_{\mathrm{w}}{\mathit{\omega }}^{\mathrm{2}}\mathit{\varphi }\\ & -\mathrm{\infty }<x<\mathrm{\infty }.\end{array}\end{array}$$

Using the separation of variable, the solution of the Laplace equation can be established as the sum of e^ikxcosh kz. Therefore, the general form of the dispersion equation can be established as

$$\begin{array}{}\text{(9)}& k\tan \left(kD\right)={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}{\mathit{\omega }}^{\mathrm{2}}}{\frac{G_E{h}^{\mathrm{3}}}{\mathrm{6}\left(\mathrm{1}-\mathit{\nu }\right)}{k}^{\mathrm{4}}+{\mathit{\rho }}_{\mathrm{w}}g-{\mathit{\omega }}^{\mathrm{2}}\left({\mathit{\rho }}_{\mathrm{i}}h+A\right)-i\left(\mathit{\omega }B\right)}}.\end{array}$$

As was mentioned, G_E is the dynamic shear stress modulus of the material having cyclic motions, which depends on the mechanical behavior of the substance. It can be a complex number. Its real component is storage modulus and its imaginary component is the loss modulus. As explained before, solid ice cover can either be assumed to be elastic or viscoelastic. As such, models have been developed for both elastic and viscoelastic ice covers. As displayed in Fig. 1d, mechanical behavior of these materials can be represented by using a spring element (which demonstrates the elastic nature of the material) and a dashpot element (which demonstrates the viscous nature of the material).

For an elastic solid body, G_E (dynamic modulus) is given by

$$\begin{array}{}\text{(10)}& G_E=G,\end{array}$$

where G is the shear modulus of the material. Loss modulus of an elastic material is zero (there is no imaginary component).

For a viscoelastic material, two different two-parameter models can be used. The first one is Kelvin–Voigt (KV) and the second one is Maxwell. For the former, displacements of both elements (i.e., spring and dashpot) are similar, but the resulting forces are different. For a Maxwell material, however, forces emerging in elements are similar, and displacements are different.

Figure 1 A pictograph of the problem and the different patterns of rheological behavior considered for the ice layer. Panels (a) and (b) repetitively show a dry beam and waves propagation in an open-water sea. Panel (c) shows waves propagating in an ice-covered sea. Panel (d) shows different rheological behavior that can be considered for the ice.

For a KV material, the dynamic shear modulus can be written as

$$\begin{array}{}\text{(11)}& G_E=G-i\mathit{\omega }\mathit{\eta }.\end{array}$$

Here, η is the dynamic viscosity of the material (Serra-Aguila et al., 2019). The real and imaginary components of G_E are the storage and loss moduli, as explained before. The former is responsible for harmonic response, and the latter is responsible for energy dissipation.

For a Maxwell material, the dynamic shear modulus is given by

$$\begin{array}{}\text{(12)}& G_E={\displaystyle \frac{G\left({\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}\right)}{\left(\mathrm{1}+{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}\right)}}-i{\displaystyle \frac{G\left(\mathit{\tau }\mathit{\omega }\right)}{\left(\mathrm{1}+{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}\right)}},\end{array}$$

where $\mathit{\tau }=\mathit{\eta }/G$ is the relaxation time (Serra-Aguila et al., 2019). It is worth noting that dynamic viscosity can affect the storage modulus of Maxwell materials, and shear modulus can affect the loss modulus of Maxwell materials.

Dispersion relationships can be established for different materials. The first dispersion relationships is developed for water waves traveling into a pure elastic material as

$$\begin{array}{}\text{(13)}& {\mathit{\omega }}^{\mathrm{2}}=\left({\displaystyle \frac{G{h}^{\mathrm{3}}}{\mathrm{6}\left(\mathrm{1}-\mathit{\nu }\right){\mathit{\rho }}_{\mathrm{w}}}}{k}^{\mathrm{4}}+g-i\mathit{\omega }{\displaystyle \frac{B}{{\mathit{\rho }}_{\mathrm{w}}}}-{\displaystyle \frac{({\mathit{\rho }}_{\mathrm{i}}h+A)}{{\mathit{\rho }}_{\mathrm{w}}}}{\mathit{\omega }}^{\mathrm{2}}\right)k\tanh \left(kD\right).\end{array}$$

Another version of this relationship has been previously documented by Mosig et al. (2015). The only difference between Eq. (13) and the one presented in Mosig et al. (2015) is the added mass term, which may slightly increase the wavenumber, which shall be discussed later. This model is called RP. Note that the RP model (dispersion relationship of pure elastic material with zero added mass) was presented in some other studies (e.g., Squire et al., 2009, and Williams et al., 2013). But note that this model only considers one dissipative mechanism (fluid damping), though in the present research it is aimed to show that coexistence of two different dissipative mechanism (fluid-based and solid-based) help us predict the decay rate and dispersion relationship. The rest of the dispersion relationship formulated in the present research incorporates fluid-based and solid-based mechanisms. One of the reasons that the pure elastic model is presented in this research is to show how the results of viscoelastic models may differ from those of this model.

The second dispersion relationship describes the link between frequency and wavenumber for a KV material, which is found to be

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\mathit{\omega }}^{\mathrm{2}}& =\left({\displaystyle \frac{G{h}^{\mathrm{3}}}{\mathrm{6}(\mathrm{1}-\mathit{\nu }){\mathit{\rho }}_{\mathrm{w}}}}{k}^{\mathrm{4}}-i{\displaystyle \frac{\mathit{\omega }\mathit{\eta }{h}^{\mathrm{3}}}{\mathrm{6}(\mathrm{1}-\mathit{\nu }){\mathit{\rho }}_{\mathrm{w}}}}{k}^{\mathrm{4}}\right.\\ & \left.+g-i\mathit{\omega }{\displaystyle \frac{B}{{\mathit{\rho }}_{\mathrm{w}}}}-{\displaystyle \frac{({\mathit{\rho }}_{\mathrm{i}}h+A)}{{\mathit{\rho }}_{\mathrm{w}}}}{\mathit{\omega }}^{\mathrm{2}}\right)k\tanh \left(kD\right).\end{array}\end{array}$$

Equation (14) with B=0 is known as FS model (Mosig et al., 2015). The dispersion relationship of a Maxwell model is built as

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\mathit{\omega }}^{\mathrm{2}}& =\left({\displaystyle \frac{G\left({\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}\right)h^{\mathrm{3}}}{\mathrm{6}(\mathrm{1}-\mathit{\nu })\left(\mathrm{1}+{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}\right){\mathit{\rho }}_{\mathrm{w}}}}{k}^{\mathrm{4}}\right.\\ & -i{\displaystyle \frac{G\left(\mathit{\tau }\mathit{\omega }\right)h^{\mathrm{3}}}{\mathrm{6}(\mathrm{1}-\mathit{\nu })\left(\mathrm{1}+{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}\right){\mathit{\rho }}_{\mathrm{w}}}}{k}^{\mathrm{4}}+g-i\mathit{\omega }{\displaystyle \frac{B}{{\mathit{\rho }}_{\mathrm{w}}}}\\ & \left.-{\displaystyle \frac{({\mathit{\rho }}_{\mathrm{i}}h+A)}{{\mathit{\rho }}_{\mathrm{w}}}}{\mathit{\omega }}^{\mathrm{2}}\right)k\tanh \left(kD\right).\end{array}\end{array}$$

Each of the above equations provides us with the roots of the dispersion relationships. We choose the dominant root with the positive real part, which refers to the wave advancing in the solid body. The dominant root is

$$\begin{array}{}\text{(16)}& k=k_{\mathrm{i}}-i\mathit{\alpha },\end{array}$$

and it is found through using a numerical method with an initial guess. All the dispersion relationships have multiple roots. We seek the one found with an initial guess which is set to be very close to the wavenumber of the open-water condition. In Eq. (16), k_i is the wavenumber in the ice-covered sea, and α is the wave height decay rate, which were introduced before. If the initial guess is set to be much greater than the open-water wavenumber, abrupt changes in the α vs. λ_o may emerge. Readers who are interested in finding the root of the dispersion relationship are referred to Zhao et al. (2017) and Das (2022).

2.2 Non-dimensional representation

To analyze results and to use any of the models more easily, all parameters are normalized. Parameters are normalized using the Buckingham Pi theorem, enabling us to perform scaling law. For the first dispersion relationship, eight parameters are involved, but for the second and third relationships, nine parameters are involved. Therefore, five non-dimensional numbers are identified for the pure elastic model, and six non-dimensional numbers are identified for models developed for viscoelastic materials.

Non-dimensional numbers have been previously presented by different authors, concerted with the field wave–mud or wave–ice interaction (Jain and Mehta, 2009; Yu et al., 2019, 2022). The first non-dimensional number represents the non-dimensional wavelength in an open-water condition, which is given by

$$\begin{array}{}\text{(17)}& {\widehat{\mathit{\lambda }}}_{\mathrm{o}}={\mathit{\lambda }}_{\mathrm{o}}/h.\end{array}$$

The second non-dimensional number is the non-dimensional wavenumber of waves propagating in a covered sea condition, given by

$$\begin{array}{}\text{(18)}& {\widehat{k}}_{\mathrm{i}}=k_{\mathrm{i}}/k_{\mathrm{o}}.\end{array}$$

Here, k_o is the open-water wave number. The attenuation rate is normalized as

$$\begin{array}{}\text{(19)}& {\widehat{\mathit{\alpha }}}_{\mathrm{i}}=\mathit{\alpha }/k_{\mathrm{o}}.\end{array}$$

The third non-dimensional number is the elasticity per unit of mass, which can be formulated as

$$\begin{array}{}\text{(20)}& \widehat{G}=G/{\mathit{\rho }}_{\mathrm{i}}gh.\end{array}$$

In the present research, $\widehat{G}$ is called elasticity number (which is inspired by Yu et al. (2019)). The other non-dimensional number is

$$\begin{array}{}\text{(21)}& \widehat{\mathit{\eta }}=\mathit{\eta }/{\mathit{\rho }}_{\mathrm{i}}\sqrt{g{h}^{\mathrm{3}}}.\end{array}$$

The fourth non-dimensional number is the relative density of the ice and is calculated by

$$\begin{array}{}\text{(22)}& \widehat{\mathit{\rho }}={\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}}.\end{array}$$

The hydrodynamic damping coefficient is normalized by

$$\begin{array}{}\text{(23)}& \widehat{B}=B/{\mathit{\rho }}_{\mathrm{w}}\sqrt{gh}.\end{array}$$

The added mass coefficient can be normalized by

$$\begin{array}{}\text{(24)}& \widehat{A}=A/{\mathit{\rho }}_{\mathrm{w}}h.\end{array}$$

3 Results and discussions

Results of wet beam models are presented in five separate sub-sections. The first sub-section presents a discussion on the effects of different parameters on the dispersion and decay rate plots. The primary aim is to provide a better understanding of the sensitivity of models to mechanical behavior of material and fluid forces. The second and third sub-sections discuss the ability of models in reconstruction of the decay rate and dispersion process through comparing their results against field and flume measurements. A brief summary of these measurements is documented in Table 1. In the fourth sub-section it is attempted to understand whether any other solid model can be used to formulate the dispersion relationship or not. The final sub-section of results presents a deep discussion on models and their abilities by covering limitations of models.

Voermans et al. (2021)Wadhams et al. (1988)Meylan et al. (2014)Yiew et al. (2019)Sree et al. (2018)

Table 1 Cases studied in the present paper.

3.1 Effects of different parameters on decay rate and dispersion process

Figure 2 shows the normalized decay rates calculated by setting different values for the mechanical properties of the cover. Left and right columns respectively show the data found for KV and Maxwell materials. The dashed and dotted curves denote the results of the dispersion relation of pure elastic (PE) material (i.e., RP model with added mass contribution). The first and second rows show the effects of dynamic viscosity on the decay rate. The results presented in the first row correspond to a material with $\sqrt{\widehat{G}}=\mathrm{340}$ (large flexural rigidity) and the ones plotted in the panels of second row correspond to a material with $\sqrt{\widehat{G}}=\mathrm{3.4}$ (low flexural rigidity).

As seen, for a viscoelastic material with a larger flexural rigidity, $\widehat{\mathit{\alpha }}$ increases under the increase in ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ reaching a maximum value at a critical wavelength. With the increase in ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$, $\widehat{\mathit{\alpha }}$ decreases. For a KV material, the critical ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ is sensitive to the dynamic viscosity, growing with the increase in the dynamic viscosity. For a Maxwell material, the increase in the dynamic viscosity reduces the decay rate. Interestingly, the pure elastic model gives negligible decay rates at small dimensionless wavelengths (${\widehat{\mathit{\lambda }}}_{\mathrm{o}}<\mathrm{50}$). The dimensionless fluid damping used to calculate the decay rates is set to be 0.0032.

Figure 2 Effects of the dynamic viscosity (a–d) and elasticity number (e, f) on $\widehat{\mathit{\alpha }}$ vs. the ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves different solid models give. Dashed and dotted models denote the decay rates that a PE material (RP model with consideration of added mass force) predicates.

To provide a clear picture of the effects of the fluid damping on the decay rates, $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curve that the pure elastic model gives is also plotted. The decay rate of a pure elastic material peaks at a specific ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$. Interestingly, the decay rates of the viscoelastic materials converge to that of the pure elastic model, signifying that the contribution of solid-based energy damping vanishes at relatively long waves. Instead, the fluid-based energy damping is dominant over the range of long waves. Note that the peak observed in $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves constructed using KV and Maxwell models is due to the method used for normalizing the decay rate, and will not emerge if dimensional data are plotted. The dimensional data are presented later, and it will be shown that the decay rate of the pure elastic material peaks in short-wave range, though that of viscoelastic material does not.

Now we discuss the decay rates of covers with a lower flexural rigidity (second row). No critical ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ is observed in $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves the KV model gives for a low elasticity number. Decay rate of the KV model with low rigidity reduces with an increase in dimensionless wavelength. The increase in the dynamic viscosity can affect the decay rate, though its effects are noticeable at shorter waves. For longer waves; however, different values of $\widehat{\mathit{\eta }}$ give similar decay rates, which match with the those of a pure elastic material. This matches with what was observed for the larger rigidity (upper row). For a Maxwell material with low rigidity, a critical ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ emerges over the short wavelength range. Similar to KV materials, the decay rate of Maxwell materials with low rigidity are sensitive to the dynamic viscosity over a very short range of wavelength. The $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves found by setting different values for the dynamic viscosity converge to each other and finally align to the decay rates of pure model. This again confirms that fluid-based energy damping becomes dominant with the increase in the dimensionless wavelength.

The last row of Fig. 2 compares the decay rates of materials with different elasticity numbers. As apparent, when a KV model is used, an increase in elasticity number reduces the decay rate. The most significant effects of elasticity on decay rate emerge at short wavelengths, where the solid-based energy damping is expected to be dominant. The decay rates of KV materials with different elasticity numbers converge to those of the pure elastic material. The decay rate of a Maxwell material is proportional to its elastic modulus. Similar to decay rate plots of KV materials with different elasticity numbers, decay rates of Maxwell materials with different elasticity numbers converge to what the pure elastic model predicts. The $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves of pure elastic material with a larger flexural rigidity, ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ peaks at a longer ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$. This means that, in case a RP model highly under-predicts the decay rate with setting a realistic elastic modulus, a much lower elastic modulus may be used for fitting the experimental data with the dispersion relationship. This has been observed by Mosig et al. (2015). They demonstrated that a shear modulus of $\approx \mathrm{3.2}\times {\mathrm{10}}^{\mathrm{7}}$ Pa gives the best fitting, though the shear modulus of ice is expected to be greater than 1 GPa.

Figure 3 Effects of fluid damping coefficient on $\widehat{\mathit{\alpha }}$ vs. the ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves viscoelastic models give.

Figure 3 shows how the consideration of fluid damping can affect the decay rates. Three different plots are presented, each of which shows the decay rates found by setting different values for the fluid damping coefficient. The solid blue curves show the data found by setting a zero fluid damping coefficient, i.e., this plot represents the decay rates calculated in the absence of fluid damping. The two other plots show the decay rates predicted by setting two different values for the damping coefficient.

The left panel of Fig. 3 shows the data related to a KV material. As apparent when the fluid damping coefficient is set to be zero, the decay rate decreases with a very high rate at ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}>\mathrm{180}$. By assuming a non-zero value for the damping coefficient, one abrupt reduction in the $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curve occurs, though the rate of the reduction of the decay rate as a function of the dimensionless wavelength decreases with the increase in the wavelength. Eventually the decay rate decreases with a mild rate over the range of ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}>\mathrm{400}$. This confirms that the fluid-based energy damping becomes dominant over this range. For a large fluid-damping coefficient, as seen, the sudden reduction in the decay rate does not occur and decreases at a low rate after reaching its peak value.

The decay rates of a Maxwell material considering different damping coefficients are plotted in the right panel. Trends of the presented curves are consistent with the ones observed in the left panel of Fig. 3 (KV material). This demonstrates that the effects of fluid-damping on decay rates are insensitive to the nature considered for material, which matches with the presented formulations for the dispersion relationships.

Figure 4 Log–log plots of α vs. ω curves constructed using different models. Curves are constructed by setting shear modulus to be 1 GPa.

To understand the role of fluid-based damping and the effects of viscosity on the decay rate more deeply, α vs. ω curves are plotted on a log–log scale. Dissimilar to previous figures, the dimensional data are presented as they help to understand the dependency of α as a function of ω. The α vs. ω curve corresponding to pure elastic material (RP model with added mass) peaks at a frequency of ≈1 rad s⁻¹, though curves constructed using the Maxwell and KV models do not peak. Note that the present data are dimensional, and if the dimensionless data (i.e., $\widehat{\mathit{\alpha }}$) were presented, all curves would peak as was observed before (see Fig. 2).

α values predicted using the pure are proportional to ω³ in the long-wave regime. This matches with the analysis presented by Meylan et al. (2018). In the short-wave range, decay rate of pure elastic material is proportional to ω^−0.35. We recall that it was discussed that one solution to match predictions of RP model with experimental data over short-range wave is to decrease the shear modulus of material, which may significantly widen the zone over which α is proportional to ω³.

The α vs. ω curves of viscoelastic materials are approximately proportional to ω³ over a wider range of frequency as compared to the curves constructed using the pure elastic model. In the short-wave zone, α values predicted using the KV model are proportional to ωⁿ, where n is slightly lower than 0.5. Overall, it can be concluded that consideration of solid-based damping widens the frequency range over which α∝ω³ and cancels out the decrease of decay over short-wave range.

The dispersion relationships are also employed to calculate wavenumbers of waves propagating into viscoelastic covers. Results are depicted in Fig. 5. The density of covers is set to be ${\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{w}}=\mathrm{0.9}$, which is very close to that of the sea ice. The two upper panels of Fig. 5 show the dimensionless wavenumbers found for KV materials with different elasticity numbers. Zero added mass coefficients (solid curve) and two non-zero added mass coefficients are considered.

Figure 5 Effects of the added mass force (a, b) and the dynamic viscosity (c, d) on the dispersion curves of waves propagating in an ice cover.

As apparent, when the elasticity number of a KV material is greater (left) and the added mass is nil, the dimensionless wavenumber is below 1.0 at small values of ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$. This means that, compared to an open-sea condition, gravity waves traveling into a solid cover become longer over the range of a short open-water wavelength. ${\widehat{k}}_{\mathrm{i}}$ increases with the increase in ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ and eventually converges into 1.0, where plotted curves are flattened out, meeting the reference line (the dashed blue line). Interestingly, added mass can affect the wavenumber. Under the action of added mass force, wavenumber becomes greater. The effects of added mass on the wavenumber are more significant over $\mathrm{100}<{\widehat{\mathit{\lambda }}}_{\mathrm{o}}<\mathrm{600}$. When the added mass coefficient is set to be very large, the dimensionless wavenumber becomes greater than 1.0, reaching a maximum value, and then decreases, converging into 1.0. The same conclusion can be made for a pure elastic model as the dispersion process of pure elastic material and KV model is similar.

For a KV material with low density and a zero added mass coefficient, the wavenumber is below 1.0 when the dimensionless open-water wavelength is small. Compared to a KV material with a greater elasticity number (left panel), ${\widehat{k}}_{\mathrm{i}}$ becomes lower than 1.0 over a narrower range of ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$. This implies that waves propagating into a solid cover are lengthened over a wider range of dimensionless open-water wavelength when flexural rigidity increases. ${\widehat{k}}_{\mathrm{i}}$ becomes greater than 1.0 and reaches a peak value over the range of short the open-water wavelength, which is in contrast with what was observed for material with a larger elasticity number and zero added mass (left panel).

The added mass force can affect the dispersion process of waves propagating into cover with lower elasticity number by increasing the wavenumber, though its influences on the waves advancing in the cover with lower elasticity number are more noticeable compared to cover with larger elasticity number.

The two lower panels of Fig. 4 show the calculated dimensionless wavenumbers of gravity waves propagating into a solid cover hypothesized to behave in the same way Maxwell materials do. Left and right panels respectively show the curves corresponding to covers with relatively large and low elasticity numbers. The trends of curves plotted in the left panel confirm that the dynamic viscosity can slightly affect the wavenumber of waves advancing in the cover with a large elasticity number. The influence of dynamic viscosity on wavenumber is not significant and emerges over a relatively narrow range, (the close-up view provides the evidence). Dimensionless wavenumber is seen to be insensitive to the dynamic viscosity of Maxwell materials when Elasticity number is low (right panel).

3.2 Ability of models in prediction of the decay rates

This section presents comparisons between predictions of the models and decay rates found through field and flume measurements. In all runs, the $\widehat{A}$ is set to be 1.0. Decay rates of two recent field tests are presented in Fig. 6 (circle markers). Upper and lower rows respectively display the data corresponding to field tests that took place in the Arctic and Antarctic. The first, second and third columns of Fig. 6 respectively show the decay rates predicted by prescribing Kelvin–Voigt, Maxwell and pure elastic materials. The decay rates found by setting a zero damping coefficient are also plotted to demonstrate how the inclusion of the fluid damping can improve the accuracy of models in the prediction of the decay rates.

Viscoelastic models cannot follow the field data when fluid damping is set to be zero, and tails of curves diverge from those of the field data. This can be seen in both upper and lower rows of Fig. 6. In this condition, the increase of the dynamic viscosity cannot affect the trend of $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$, and may only shift the curve upward or downward, depending on the nature considered for the material. But, when the fluid damping coefficient is set to be non-zero, the $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ follow the field data. This confirms that fluid damping can contribute to energy damping occurring under a landfast ice cover. Interestingly, fluid damping coefficients used to predict the decay rates of viscoelastic are similar.

For the KV material, a dynamic viscosity of $\approx \mathrm{1.3}\times {\mathrm{10}}^{\mathrm{8}}$ Pa s gives the best fitting for both ice covers. But for a Maxwell material, a dynamic viscosity of $\approx \mathrm{4.4}\times {\mathrm{10}}^{\mathrm{10}}$ Pa s gives the best fitting. Finally, the pure elastic material assumption can also be used to match curves with the field experiments. But, compared to KV and Maxwell models, a greater value of fluid damping coefficient gives the best fitting for an elastic model.

Figure 6 Comparisons between $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves predicted by different models and the data measured by Voermans et al. (2021). Upper and lower rows respectively show data measured in the Arctic and Antarctic. In the first and second columns, decay rates predicted by setting a zero fluid damping are also plotted. The shear modulus of ice is set to be 1 GPa.

The decay rates of two different broken ice fields are calculated using the wet beam models and are compared against measurements. Field data and calculated decay rates are plotted in Fig. 7. The data presented in the upper and lower panels of this figure are respectively documented in Wadhams et al. (1988) and Meylan et al. (2014).

Wet beam models established for viscoelastic materials can predict decay rate curves of a broken ice field with a good level of accuracy when fluid damping is incorporated into the calculations. Similar to curves plotted in Fig. 6, increase or decrease in the dynamic viscosity of the material can shift the curves vertically, not affecting the tail of the $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curve. The curve given by a pure elastic material is not accurate at all dimensionless wavelengths. This implies that to accurately compute the decay rate of a broken ice field, solid-based energy damping caused by the viscoelastic behavior of the solid needs to be taken into consideration, which is being lacked in a pure elastic model. The dimensional data and curves related to field tests documented in Wadhams et al. (1988) are presented in Appendix A. The aim is to show whether a peak may emerge in α vs. T curves when models are used or not, as the emergence of a peak value in high-frequency range may not be real (Thomson et al., 2021).

Figure 7 Comparisons between $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves predicted by different models and the data measured by Wadhams et al. (1988), upper row, and Meylan et al. (2014), lower row. In the first and second columns, decay rates predicted by setting a zero fluid damping are also plotted. The shear modulus of ice is set to be 1 GPa.

Wet beam models are employed to calculate the decay rates of waves traveling into viscoelastic covers, with an aim to understand whether they can be used to predict attenuation rates measured in flume tests or not. The reconstructed $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves are presented in Fig. 8. The upper and lower rows show the data correspond to covers with low and large rigidity.

The curves reconstructed by viscoelastic models cannot follow the experimental data if the fluid damping is not considered. Dashed and dashed-dotted curves plotted in the first and second columns provide shreds of evidence. Their results only match with experimental data at very short dimensionless open-water wavelengths. Viscoelastic-based models can accurately predict decay rates when the fluid damping is set to be non-zero. The elastic model can calculate decay rates with an acceptable level of accuracy. The fluid damping that gives the best fitting for pure elastic material is slightly greater than those of models built for viscoelastic materials.

Figure 8 Comparison between $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves predicted by models and the data measured by Sree et al. (2018). Upper and lower rows show the data related to a cover with a shear modulus of 20 and 80 KPa respectively.

At the final stage, the decay rates of the freshwater ice formed in the wave flume of the University of Melbourne are predicted by using the presented models. Figure 9 displays the decay rates plots and experimental data. The results presented in the upper and lower rows of Fig. 9 are related to 1 and 1.5 cm thick ice covers. Before conducting flume experiments, Parra et al. (2020) measured the Young's modulus of the dry freshwater ice, reporting that Young's modulus is ≈3 GPa.

The KV model can accurately predict the decay rates of 1 cm thick ice when dynamic viscosity is set to be 2.6×10⁹ Pa s. Inclusion of the fluid damping is seen to be ineffective when the KV model is used. The Maxwell model can also predict the decay rate curve of the 1 cm thick ice cover with a good level of accuracy. But the decay rates that the Maxwell model predict depend on the fluid damping. The model developed for pure elastic materials can construct $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curve if fluid damping is set to be 1300 Pa s m⁻¹, which is 13 times greater than those considered for the KV and Maxwell models. This significant difference between fluid damping of pure elastic models and viscoelastic models, which was never observed in previous Figures, indicates that the solid-based damping highly contributes to energy dissipation over the range of tested waves in the laboratory. To compensate for the absence of the solid-based energy damping of the pure elastic model, a very large fluid damping coefficient needs to be used. This is in contrast with what was observed for the tests of Sree et al. (2018), where the difference between fluid damping of elastic and viscoelastic models was not significant. The big difference between the elasticity numbers of the covers tested by (Sree et al., 2018) and ice covers tested by Yiew et al. (2019) explains this behavior. It was also demonstrated that when the elasticity number is low, the solid-based energy damping contributes to energy dissipation over a narrow range of wavelengths (second row of Fig. 2), compared to a cover with a larger elasticity number (first row of Fig. 2).

The KV model under-predicts the decay rate curve of 1.5 cm thick ice cover regardless of the values of the dynamic viscosity and the fluid damping coefficients. The Maxwell model, however, can construct $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curve when a dynamic viscosity value of 4×10⁶ Pa s is set. This value is lower than that the dynamic viscosity gave the best fitting for the thinner ice. Artificial effects, boundary conditions and the presence of side walls may cause a larger energy damping pattern when fluid interacts with the thicker ice (Sutherland et al., 2016). Thus, different values of dynamic viscosity give the best fitting for 1 and 1.5 cm thick ice covers. It is important to note that dynamic viscosity values giving the best fitting for the landfast ice were not seen to be different (Fig. 6), where artificial effects are less likely to contribute. As seen, fluid damping does not have any significant effect on the curve over the range of tested waves. It can only affect the tail of curve. This again confirms that the solid-based energy damping is dominant over the range of tested waves. If tests covered a wider range of open-water wavelengths, especially longer ones, a proper value for fluid damping could be found through fitting predicted curves with experimental data. When the pure elastic model is used to calculate decay rates, a very large fluid damping coefficient needs to be set, but the curve can never meet the experimental data.

Figure 9 Comparison between $\widehat{\mathit{\alpha }}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves predicted by models and the data measured by Yiew et al. (2019). Upper and lower rows show the data related to freshwater ice covers with thicknesses of 1 and 1.5 cm respectively.

Figure 10 Comparison between ${\widehat{k}}_{\mathrm{i}}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves predicted by wet beam models and the data measured by Voermans et al. (2021). Upper and lower rows respectively show data related to tests performed in Arctic and Antarctic. The shear modulus of ice is set to be 1 GPa.

3.3 Ability of models in the prediction of the dispersion process

The accuracy of models in the calculation of the dispersion process of waves traveling through viscoelastic covers is also evaluated. First, the dispersion curves of waves propagating into landfast ice are constructed, and then the dispersion curves of waves advancing in freshwater ice are plotted.

Figure 10 shows the normalized wavenumber ${\widehat{k}}_{\mathrm{i}}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves of landfast ice. Symbols indicate the field measurements and curves denote the calculated dimensionless wavenumbers. Left and right panels respectively show predictions made for KV and Maxwell materials. As seen, both models can predict the dispersion process of waves traveling into the landfast ice with an acceptable level of accuracy. The inputs that are used for the construction of dispersion curves are similar to what gives the best fitting for decay rates.

Figure 11 Comparison between ${\widehat{k}}_{\mathrm{i}}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves predicted by wet beam models and the data measured by Yiew et al. (2019). Upper and lower rows show the data related to freshwater ice covers with thicknesses of 1 and 1.5 cm respectively.

Figure 11 shows the ${\widehat{k}}_{\mathrm{i}}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ plots of the flume tests performed in the University of Melbourne. It was observed that none of the models can predict the dispersion process of the waves interaction with freshwater ice with the inputs giving the best fitting for the decay rates plots (Fig. 9). This is in contrast with what was observed in Fig. 10, where models accurately predicted the wave dispersion process with similar inputs utilized for the calculation of decay rates. As seen when the dispersion process of waves interacting with freshwater is constructed by using a KV model, the Young's modulus values of 5×10⁵ and 2.5×10⁸ Pa give the best fitting for the dispersion processes of 1 and 1.5 cm ice covers respectively. As explained in the Introduction, some other researchers observed that the effective Young's modulus of a material is different and smaller from what is measured in dry tests. For a Maxwell material, however, different values of dynamic viscosity can be used to fit the curves with the experimental data. It was previously demonstrated that the dynamic viscosity can affect ${\widehat{k}}_{\mathrm{i}}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ plots by increasing the wavenumber over the range of short waves. By setting larger values for the dynamic viscosity, the constructed dispersion curve matches with the experimental data (the solid curve).

Note that wavenumbers of field tests documented in Wadhams et al. (1988) and Meylan et al. (2014) are not presented. Therefore, we were not able to compare the capability of wet beam models in the predictions of wavenumbers of the field with broken ice against those of a broken ice field.

3.4 Other models

In Sect. 3.2 and 3.3, it was observed that proposed models cannot accurately calculate the dispersion process and decay rates of waves advancing into the freshwater ice. Artificial effects may have contributed to the dispersion and dissipation process, though the large difference between inputs giving the best fit for the decay rate and dispersion plots still leaves us with a big question mark about the main reason for such a difference.

As discussed, other researchers have hypothesized that the effective elasticity or dynamic viscosity of the wet ice interacting with water waves can be different from what have been measured in dry tests. While the whole paper looks into ability of KV and Maxwell in the prediction of decay rates and dispersion process and discusses how accurate their results can be, the idea of developing other models by considering various viscoelastic materials arises. We can use the other solid models to evaluate whether the Young's modulus giving the best result for decay rates and dispersion process matches with what was found in dry tests or not.

Figure 12 Standard linear solid models. Panels (a) and (b) respectively show the SLS KV and SLS Maxwell models.

To provide an answer to the above question, the linear combination of KV and Maxwell materials, known as the standard linear solid (SLS) model, is used to describe the mechanical behavior of ice. Two spring elements (G₁ and G₂) are used to formulate these models. Standard linear solid models are also called “Zener”. Two different Zener models have been introduced in this paper. The first one is an SLS model established using a Kelvin approach (a spring element in series with a KV arm, also known as the first order generalized Kelvin–Voigt model). The other model is acquired by employing a Maxwell approach (a spring element in parallel with a Maxwell arm, also known as the first order generalized Maxwell model). The schematic of both of these models is shown in Fig. 12. The storage and loss moduli of these two solid models are formulated as

$$\begin{array}{}\text{(25)}& G_E=G_{\mathrm{0}}(\mathrm{1}-{\displaystyle \frac{p}{\mathrm{1}+{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}}})-i{G}_{\mathrm{0}}{\displaystyle \frac{p\mathit{\tau }\mathit{\omega }}{\mathrm{1}+{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}}}.\end{array}$$

For the SLS KV material, G₀, p and τ are given by

$$\begin{array}{}\text{(26)}& \begin{array}{lll}{G}_{\mathrm{0}}=G_{\mathrm{1}},& p=\frac{G_{\mathrm{1}}}{G_{\mathrm{1}}+G_{\mathrm{2}}},& \mathit{\tau }=\frac{\mathit{\eta }}{G_{\mathrm{1}}+G_{\mathrm{2}}}.\end{array}\end{array}$$

For the SLS Maxwell material, G₀, p and τ are calculated as

$$\begin{array}{}\text{(27)}& \begin{array}{lll}{G}_{\mathrm{0}}=G_{\mathrm{1}}+G_{\mathrm{2}},& p=\frac{G_{\mathrm{2}}}{G_{\mathrm{1}}+G_{\mathrm{2}}},& \mathit{\tau }=\frac{\mathit{\eta }}{G_{\mathrm{2}}}.\end{array}\end{array}$$

Readers who are interested in the above formulations are referred to Serra-Aguila et al. (2019). The same method used to establish dispersion relationships in Sect. 2 is applied, and the related dispersion relationship for SLS materials is formulated as

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{\mathit{\omega }}^{\mathrm{2}}& =({\displaystyle \frac{G_{\mathrm{0}}(\mathrm{1}-\frac{p}{\mathrm{1}+{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}})h^{\mathrm{3}}}{\mathrm{6}(\mathrm{1}-\mathit{\nu }){\mathit{\rho }}_{\mathrm{w}}}}{k}^{\mathrm{4}}-i{\displaystyle \frac{G_{\mathrm{0}}\frac{p\mathit{\tau }\mathit{\omega }}{\mathrm{1}+{\mathit{\tau }}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}}{h}^{\mathrm{3}}}{\mathrm{6}(\mathrm{1}-\mathit{\nu }){\mathit{\rho }}_{\mathrm{w}}}}{k}^{\mathrm{4}}\\ & +g-i\mathit{\omega }{\displaystyle \frac{B}{{\mathit{\rho }}_{\mathrm{w}}}}-{\displaystyle \frac{({\mathit{\rho }}_{\mathrm{i}}h+A)}{{\mathit{\rho }}_{\mathrm{w}}}}{\mathit{\omega }}^{\mathrm{2}}\left)k\tanh \right(kD).\end{array}\end{array}$$

Decay rates and dispersion of waves propagating into the 1.5 cm freshwater ice cover are recalculated by using the two introduced SLS models and plotted in Fig. 13. The left and right panels respectively show the results found using the SLS KV and SLS Maxwell models. Dispersion curves are plotted in the upper panels, and the decay rates are illustrated in the lower panels. SLS KV and SLS Maxwell models both construct ${\widehat{k}}_{\mathrm{i}}$ vs. ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ plots similarly. Interestingly, the equivalent Young's modulus of the SLS KV model is 2.5× 10⁴ Pa and that of the SLS Maxwell model is 3.3× 10⁹ Pa. This confirms that the SLS KV model, which is built using the KV approach, can only predict the predicted process with a very low equivalent Young's modulus, which is much smaller than the Young's modulus found in dry tests (which is around 2.6×10⁹ Pa). But the equivalent Young's modulus of the SLS Maxwell model is close to what is found in dry tests. This signifies that the freshwater ice formed in the flume is more likely to behave similarly to an SLS Maxwell material. This can explain the difference between the effective Young's modulus reported in different experimental research. Researchers who conducted those flume/basin experiments concluded that a lower Young's modulus should be used to calculate the dispersion process of waves propagating into ice when a dispersion equation built on the basis of pure elastic material (or the KV material) is employed. But as observed here, the material is more prone to show a rheological behavior falling in between those of the KV and Maxwell materials, arranged using a Maxwell approach, the equivalent Young's modulus of which is close to what is measured in dry tests.

Figure 13 Comparison between predictions of SLS models and values measured by Yiew et al. (2019). Upper and lower rows show the decay rates and dispersion plots. The equivalent Young's modulus of the SLS KV model is 2.5×10⁴ Pa, and the equivalent Young's modulus of the SLS Maxwell model is 3.3×10⁹ Pa. Dynamic viscosity of the SLS KV model is 4×10⁷ Pa s, and dynamic viscosity of the SLS Maxwell model is 3.5×10⁹ Pa s.

The decay rates are seen to be well predicted by both models. The interesting point is that the KV model was found not to be able to construct the decay rates of 1.5 cm thick. But as seen here, the results of the SLS KV model fairly agree with experimental measurements. The SLS Maxwell model can predict the decay rates very well. The trend of the decay rates given by SLS Maxwell and SLS KV are very similar.

Note that adding more elements may make the dispersion relationship more complex, leaving us with more options (i.e., more inputs), though it can also lead to over-fitting as more parameters are needed to be tuned. In the future, more studies can be carried out to understand the mechanical behavior of different types of viscoelastic ice, and it can be investigated whether it is required to use solid models with more than two elements or not.

3.5 Final discussion

Decay rates and dispersion curves of gravity waves propagating into an ice cover can be constructed using the presented wet beam models. With the same setups (e.g., dynamic viscosity and Young's modulus), different models can give different curves for the decay rates and the dispersion process. Thus, the choice of the model for the prediction of the decay rates and the dispersion process is very important. This needs a better understanding of the ice mechanics and the way it is formed. Simply stated, the mechanics of ice is a very complicated field of research, with lots of open questions that have not been answered yet.

Models formulated in the present research follow the basis of the Euler–Bernoulli beam theory, which describes displacements of a solid layer by assuming small motions. To provide a clear picture of this theory, transverse sections of a beam flexed due to an external/internal load should be assumed. If the displacement field follows the Euler–Bernoulli law, the normal vector of any transverse section is always parallel to the axis of the beam. That is, no local rotational motion occurs. This is different from what happens for a solid body following the Timoshenko–Ehrenfest beam theory, where rotational motions are taken into consideration (Timoshenko and Woinowsky-Krieger, 1959). Furthermore, to establish any model, the mechanical behavior of a substance should be formulated. It left us with two common choices as we are following a beam theory. The ice can be hypothesized to be either elastic or viscoelastic (we may have other choices for a solid ice layer, for example poroelastic ice). When elastic behavior is considered, Hooke's law is utilized to formulate the relationship between the stress and the resulting strain. If viscoelastic behavior is assumed, different linear models can be used. The KV and Maxwell models are commonly used for this aim, the former of which represents a viscoelastic solid body and the latter of which represents a mass of viscoelastic fluid.

Using mechanical behavior prescribed for the ice and employing an Euler–Bernoulli equation, dispersion relationships can be formulated. If the ice layer is assumed to be dry and no fluid–solid interaction is taken into account, dispersion relationships for wave motions in dry beams are acquired. If the multi-physical problem is considered, the ice layer is assumed to settle down on the water surface and fluid forces emerge. If the fluid motion is assumed to be linear, a fluid damping force can also be included in the beam motion equation, as was suggested by Robinson and Palmer (1990). The original models of Greenhill (1886) and Fox and Squire (1991) do not consider fluid damping. Fluid damping was considered in some other models, which were formulated for an elastic ice cover as dissipation was either formulated through introducing a complex term in flexural rigidity of the ice cover (solid-based energy damping) or by considering a linear damping term (fluid-based energy damping, Squire et al., 2021). But, the fluid damping and the solid damping can be considered at the same time, leading to a more general dispersion relationship compared to available ones.

In the present research, it is hypothesized that energy dissipation is caused by solid-based energy damping and fluid-based energy damping. The former is dominant over high frequencies (corresponding to short waves in an open-water condition), and the latter is dominant over small frequencies (corresponding to long waves in an open-water condition). The solid-based energy damping is caused by the viscoelastic resistance emerging in the solid and is absent for an elastic material. Regardless of the viscoelastic model used to treat the mechanical behavior of the material, the solid-based energy damping decreases with an increase in the wavelength. With an increase in elasticity number, the wavelength range over which solid-based energy damping is dominant becomes wider. This matches with physics. The body is more rigid, and solid responses can contribute to energy damping of longer waves compared to bodies with low rigidity. The fluid-based energy dissipation is generated by a velocity dependent force, decreasing with the increase in wavelength. The energy damping triggered by a pure elastic cover is only due to the presence of fluid damping force. Thus, to compensate for the lack of solid-based energy damping of an elastic cover, a larger fluid damping coefficient needs to be considered. Energy decay triggered by the landfast ice covers which have large rigidity can be computed using the viscoelastic models if fluid damping is included. If it is not, models cannot work properly, and the predicted decay curves diverge from field measurements. This well confirms that fluid-based energy damping contributes to total energy dissipation. Both viscoelastic models were seen to construct the decay curves and fairly follow field measurements. But the dynamic viscosity values were seen to be much different. A KV material may be a more realistic indicator of ice behavior as the landfast ice is expected to be solid.

Decay rates of broken ice fields were seen to be calculated by both viscoelastic models when the fluid-based energy damping is considered. This provided us with another piece of evidence for the contribution of fluid damping coefficient force to energy dissipation. Compared to landfast ice covers, different values of the dynamic viscosity were observed to give the best fitting for the decay rates of waves advancing in water partially covered with broken ice floes. The coexistence of open water and ice floes on the upper layer of the fluid domain can explain this. Models are formulated for an integrated layer of solid ice. If water is included, a mixed thin layer represents the cover. A volume fraction model may describe the dynamic viscosity of the layer. The dynamic viscosity of the water is much smaller compared to that of ice, but the water entrapped between ice floes may be turbulent, leading to eddy generation (turbulent flow can cause damping of waves, e.g., Tavakoli et al., 2022). In some mathematical models, an artificial damping term is usually introduced to incorporate the gap effects on fluid pressure (Lu et al., 2010). The KV model gives the best fitting of decay rate when the dynamic viscosity is reduced, compared to the landfast ice. In contrast, the Maxwell model gives the best fitting with a larger dynamic viscosity. Thus, what can be concluded is that if water reduces the dynamic viscosity of the upper layer, the KV model will more likely to be a better indicator of the ice behavior. Otherwise, the Maxwell model prescribes the mechanical behavior of the material. In addition, in a broken ice field, the gap between ice floes may lead to emergence of an extra damping mechanism, which is introduced as an artificial damping in studies highlighting interactions between floating bodies. In the present research, the aim was to develop a model for an integrated viscoelastic ice cover with incorporation of fluid damping, and thus the artificial damping was not employed.

The decay rates of the freshwater ice were also calculated. The pure elastic materials can reconstruct the decay rates only if a very large fluid damping coefficient is used, which may not be realistic. But, viscoelastic models were seen to predict the decay rate plots. Models, however, were not able to capture the dispersion process under freshwater ice covers with the same input observed to give the fitted decay rates. Effective values were seen to construct dispersion plots with an acceptable level of accuracy. This has been observed and reported by other scholars, who measured the wavelength and phase speed of disintegrated elastic/viscoelastic covers. The interesting point is that, when a Maxwell model is used, the dynamic viscosity can affect the dispersion process. This motivates us to build other models which formulate the storage modulus by applying the dynamic viscosity. Two available linear models were introduced. One is SLS KV and the other is SLS Maxwell. Both models include two springs and one dashpot. The equivalent Young's modulus giving the best fitting when SLS KV is used is $\approx \mathrm{2.5}\times {\mathrm{10}}^{\mathrm{4}}$ Pa, which is much smaller compared to that of real freshwater ice. But, for SLS Maxwell, the equivalent Young's modulus was seen to be $\approx \mathrm{3.3}\times {\mathrm{10}}^{\mathrm{9}}$ Pa, which is close to that of freshwater ice. While it is very interesting to add more elements to the solid model, over-fitting may happen as more elements are needed to be tuned. This leaves us with another question. Is it required to consider more elements in the solid model or has the effective shear modulus of the ice formed in the flume been highly reduced because of the ice formation process (for example, porosity may reduce effective shear modulus, Zong, 2022)?

Based on what was observed when SLS models were used, it can be concluded that a two-parameter models such as Maxwell and KV, with a realistic Young's modulus, cannot reproduce mechanical behavior of the freshwater ice formed in the flume, and its behavior is different from that of the landfast ice. There are still doubts about the behavior of broken sea ice. Both the Maxwell and KV models were found to give the best fitting. One reason is that the wave phase is not measured in most of the field tests that took place in the broken ice field, and researchers were mostly concerned with decay rates. Since the phase speed and dispersion plots are not available, the performance of models in the reconstruction of dispersion plots cannot be evaluated. Especially, this could show whether changes in the dynamic viscosity can modify the accuracy of the Maxwell model in the prediction of dispersion or not (the dispersion curves the Maxwell model gives are sensitive to dynamic viscosity). That may be still an open question for the future.

4 Conclusions

The common approach used to predict decay rates and the dispersion of waves penetrating an ice cover is to formulate a dispersion relationship of a continuum medium, the roots of which give the decay rate and the relative wavenumber in an ice-covered sea. The majority of models developed for wave–ice interaction have been developed based on two common approaches. First, ice was assumed to be a viscous fluid or a viscoelastic solid layer. Second, ice was assumed to be elastic, and a fluid damping term is used to calculate the decay rate. This paper aimed to present wet beam models by considering both of fluid-based and solid-based energy dissipation mechanisms by accommodating different rheological for the ice layer (KV and Maxwell). Thus, the models were called “wet beam”, which refers to water-based forces that are taken into consideration.

Predictions of viscoelastic models and field measurements were quantitatively compared against each other. KV and Maxwell models were seen to reconstruct the decay rates and dispersion process of landfast ice with a great level of accuracy. Decay rates were observed to be poorly predicted if the fluid-based energy damping is not taken into account, suggesting that this mechanism has a very important role in ice-induced energy decay over the range of long waves. The decay rates of unconsolidated ice fields were also seen to be accurately predicted by KV and Maxwell models by setting a non-zero value for the fluid damping coefficient.

Decay rates and wave dispersion of freshwater ice were predicted. The pure elastic model was seen to predict the decay rate with an unrealistic fluid damping coefficient. The decay rates predicted by KV and Maxwell models were seen to agree with experiments, though the dispersion plots were observed to diverge from the experimental data.

Two standard linear solid models were used and two other dispersion relationships were formulated. These relationships were seen to predict the attenuation rate and the dispersion plots with a good level of accuracy. But, the SLS model which was fundamentally based on the KV material gave the best fitting with an unrealistic Young's modulus.

Overall, the wet beam dispersion relationships developed by accommodating two-parameter solid models are able to predict the decay rates and dispersion process of ice fields. But, for a freshwater ice flume, a three-parameter solid model may increase the accuracy of the predictions. As wet beam models were observed to be capable of fitting decay rates of different field measurements, they can be employed in wave spectral modeling of polar seas in the future and can also be coupled with ice break-up models to simulate the evolution of marginal ice zones. In the future, nonlinear models can also be developed to consider the effects of wave steepness on dispersion, and other beam theories can be employed. In addition, viscosity and turbulence effects can also be incorporated into wet beam models.

Appendix A: An example of dimensional data

The data presented in Sect. 3.2 are dimensionless. As was observed in Fig. 4, α (dimensional decay rate) found using viscoelastic models may not peak. This is because $\widehat{\mathit{\alpha }}$ versus ${\widehat{\mathit{\lambda }}}_{\mathrm{o}}$ curves of the viscoelastic models are normalized using an open-water wavenumber. Recently, it has been shown that instrument noise and local non-linear wave generation of high-frequency waves may cause a peak in short-wave regimes which may not be real (Thomson et al., 2021). To make it clearer, the comparison between field measurements documented in Wadhams et al. (1988) and predictions of models are presented in a dimensional way (Fig. A1). As shown, field data reach a peak value at wave period of ≈5.15 s, though the decay rates predicted by viscoelastic models decrease under the increase of wave periods (i.e., they never peak). Interestingly, decay rate versus wave period curve predicted using the pure elastic model (RP with an additional added mass term) reaches a peak value in the short-wave regime.

Figure A1 Comparison between α versus T curves constructed using different wet beam models against those of field measurements documented in Wadhams et al. (1988). Curves are constructed using the inputs giving the best fitting, which can be seen in Legends and captions of Fig. 7.