Sum-frequency triad interactions among surface waves propagating through an ice sheet

Abstract

We study nonlinear resonant wave–wave interactions which occur when ocean waves propagate into a thin floating ice sheet. Using multiple-scale perturbation analysis, we obtain theoretical predictions of the wave amplitude evolution as a function of distance travelled past the ice edge for a semi-infinite ice sheet. The theoretical predictions are supported by a high-order spectral (HOS) method capable of simulating nonlinear interactions in both open water and the ice sheet. Using the HOS method, the amplitude evolution predictions are extended to multiple (coupled) triad interactions and a single ice sheet of finite length. We relate the amplitude evolution to mechanisms with strong frequency dependence – ice bending strain, related to ice breakup, as well as wave reflection and transmission. We show that, due to sum-frequency interactions, the maximum strain in the ice sheet can be more than twice that predicted by linearised theory. For an ice sheet of finite length, we show that nonlinear wave reflection and transmission coefficients depend on a parameter in terms of wave steepness and ice length, and can have values significantly different than those from linear theory. In particular, we show that nonlinear sum-frequency interactions can appreciably decrease the total wave energy transmitted past the ice sheet. This work has implications for understanding the occurrence of ice breakup, wave attenuation due to scattering in the marginal ice zone and the resulting ice floe size distribution.

1. Introduction

Ocean surface waves propagating from open water into ice-covered regions can cause mechanical ice breakup, which is a key process affecting the ice floe size distribution in the marginal ice zone (MIZ), the region of fragmented ice at the boundary between open water and consolidated sea ice. For small incident wave steepness, observations coincide with a scattering model based on linear wave theory where the ice is treated as a thin elastic plate (Kohout & Meylan 2008). For large wave events, however, ice breakup has been observed $O(100\ \text {km})$ past the ice edge, which is beyond where linear theory predicts wave-induced ice breakup can typically occur. A number of possible mechanisms to explain these observations have been proposed, such as nonlinear focusing due to lateral pack compression (Liu & Mollo-Christensen 1988) and inverse energy cascades characteristic of open water waves (Kohout et al. 2014). In this work, we propose a simple mechanism based on (leading-order) nonlinear resonant interactions among surface waves on the ice sheet. In particular, the dispersion relationship admits nonlinear resonant triads which generate sum-frequency interactions that significantly modify total wave amplitudes and can lead to an appreciable increase in maximum ice strain that is not predicted by linear theory.

Nonlinear wave–wave interactions among flexural-gravity waves have received less attention than, for example, periodic waves of maximum height (Forbes 1986; Vanden-Broeck & Părău 2011), solitary wave solutions (Guyenne & Părău 2012; Alam 2013) and a moving load on ice (Părău & Dias 2002; Bonnefoy, Meylan & Ferrant 2009) for the flexural-gravity wave problem. The existence of triad resonant interactions between flexural-gravity waves in a floating ice sheet of infinite extent was noted by Marchenko (1999), who also provided temporal amplitude evolution equations for the resonating components. The current work focuses on triad interactions among solely flexural-gravity waves without considering the effects of stratification, lateral ice compression and current, which add further richness to the dispersion relation and the possible triad interactions (Das, Sahoo & Meylan 2018; Bisht et al. 2022).

None of the aforementioned works consider nonlinear wave–wave interactions in a finite ice sheet. This is a problem of significant interest because the MIZ can be modelled as a collection of finite ice sheets (Kohout & Meylan 2008). The majority of work concerning nonlinear effects on wave–ice interactions in the MIZ (e.g. Liu et al. 2020; Waseda et al. 2022) uses phase-averaged models, such as WaveWatch3. To study general wave–ice interactions, we use a phase-resolved high-order spectral (HOS) method which explicitly captures the direct physical mechanisms, including wave reflection and transmission at ice edges. These phase-resolved findings are complementary to phase-averaged models and may inform their ongoing development.

The main objective of this work is to elucidate nonlinear sum-frequency triad interactions in a semi-infinite or finite-length ice sheet under specific incident ocean wave conditions. Of special interest are the spatial wave amplitude evolution and key frequency-dependent mechanisms, such as wave reflection/transmission and the maximum ice bending strain. For simplicity, we consider the well-established model of two-dimensional (2-D) potential flow partially covered by a thin ice sheet, treated as an Euler–Bernoulli beam; see e.g. Fox & Squire (1990) and Marchenko (1999) in two dimensions and Marchenko & Shrira (1991) and Kirby (1992) in three dimensions.

The problem statement and the potential flow nonlinear boundary value problem are specified in § 2. In § 3, we review the triad resonance condition and derive the multiple-scales (MS) perturbation evolution equations for the nonlinear wave amplitudes. The formulation and generalisation of the HOS method for nonlinear wave–ice interactions is outlined in § 4 along with supporting convergence tests. We present our results and findings in § 5 for semi-infinite and finite ice sheets. In § 5.1, we present theoretical (MS) and numerical results for a semi-infinite ice sheet – nonlinear wave amplitude evolutions and the variation of the ice strain, including the magnitude and location of its maximum. The HOS results show agreement with all the theoretical predictions in the domain where MS analysis is valid. The case of a finite-length ice sheet is considered in § 5.2 where, in addition to amplitude evolution and ice strain, the wave transmission and reflection associated (mainly) with the trailing edge are of special interest owing to their implications for wave attenuation in the MIZ (Wadhams et al. 1988; Kohout et al. 2020). We provide a brief summary in § 6.

2. Problem statement

We consider 2-D potential flow of wave propagation in deep water partially covered by a thin ice sheet of length $L$ and uniform thickness $h$, assumed to be small relative to the characteristic (ice) wavelength $\lambda \gg h$. We assume a linear elastic ice sheet with flexural rigidity $D=Eh^3/(12(1-\nu ^2))$, Young's modulus $E$ and Poisson ratio $\nu$. Let $x$ be the horizontal coordinate with ambient waves incident from negative $x$. For the ice sheet located in $0\leq x\leq L$, we consider two cases of interest: a semi-infinite ice sheet ($L\rightarrow \infty$) and a finite ice sheet ($L$ finite). The vertical coordinate is $z$, where $z=0$ corresponds to the mean free surface, $z$ is positive upward, and $z=\eta (x,t)$ denotes the elevation of the water ($x<0$, $x>L$) and ice ($0\leq x\leq L$) surface, which is assumed to be comparable to $h$. See figure 1.

Figure 1. Definition sketch.

The flow can be described by a velocity potential $\varPhi (x,z,t)$ satisfying the nonlinear initial boundary-value problem

(2.1a)$$\begin{gather} \nabla^2\varPhi=0 \quad \mbox{for}\ z\in(-\infty,\eta(x,t)], \quad x\in(-\infty,\infty), \end{gather}$$

(2.1b)$$\begin{gather}\eta_t+\varPhi_x \eta_x=\varPhi_z \quad \mbox{on}\ z=\eta(x,t), \quad x\in(-\infty,0) \cup (0,L) \cup (L,\infty), \end{gather}$$

(2.1c)$$\begin{gather}\varPhi_t+\frac{1}{2}|\nabla \varPhi|^2+g\eta+\frac{D}{\rho_w}\eta_{xxxx}=0 \quad \mbox{on}\ z=\eta(x,t),\quad x\in[0,L], \end{gather}$$

(2.1d)$$\begin{gather}\eta_{xx}=\eta_{xxx}=0 \quad \mbox{on}\ z=\eta(x,t), \quad x=0, L, \end{gather}$$

(2.1e)$$\begin{gather}\varPhi_t+\frac{1}{2}|\nabla \varPhi|^2+g\eta=0 \quad \mbox{on}\ z=\eta(x,t),\quad x\in (-\infty,0)\cup(L,\infty), \end{gather}$$

(2.1f)$$\begin{gather}|\nabla \varPhi| \rightarrow 0 \quad \mbox{as}\ z\rightarrow -\infty, \quad x\in(-\infty,\infty) , \end{gather}$$

where $g$ is gravitational acceleration and $\rho _w$ is water density. The initial boundary-value problem is complete with the inclusion of a radiation condition for outgoing scattered waves in the far field and suitable initial conditions for $\eta (x,0)$ and $\varPhi (x,\eta (x,0),0)$.

Note that in (2.1), we include nonlinearity in the flow equations but, consistent with $h/\lambda \ll 1$ and $\eta /h=O(1)$, we incorporate a linear Euler–Bernoulli beam model (2.1c) which does not consider the ice inertia (cf. e.g. Alam 2013). Past works regarding nonlinear flexural-gravity waves have used more complex beam models such as the nonlinear Kirchoff–Love model (e.g. Vanden-Broeck & Părău 2011) or that of Plotnikov & Toland (2011) (e.g. Guyenne & Părău 2012). To second order in wave steepness, the leading order considered in the current work, these models reduce to the linear Euler–Bernoulli beam model we use.

3. Theoretical results

We consider the conditions for nonlinear resonant triad interactions of waves on the ice surface and obtain the evolution equations of these wave components using MS perturbation analysis. As expected, the results resemble those of three-wave interactions appearing in other physical contexts, e.g. capillary-gravity waves (McGoldrick 1965) and nonlinear optics (Armstrong et al. 1962). We define the wave perturbation parameter $\varepsilon =O(\kappa a) \ll 1$, where $\kappa =2{\rm \pi} /\lambda$ and $a$ are the wavenumber and complex amplitude in ice, respectively ($\varepsilon _w$, $k$ and $A$ in open water). We expand the dependent variables to second order in $\varepsilon$:

(3.1a,b)\begin{equation} \varPhi(x,z,t)=\varepsilon \varPhi^{(1)} + \varepsilon^2 \varPhi^{(2)} + O(\varepsilon^3) ;\quad \eta(x,t)=\varepsilon \eta^{(1)} + \varepsilon^2 \eta^{(2)} + O(\varepsilon^3). \end{equation}

For the general triad interaction case, with wavenumber $\kappa _j$ and frequency $\omega _j$ for $j=1,2,3$, the first-order solution is

(3.2a)$$\begin{gather} \varPhi(x,z,t) = \sum_{j=1}^{3} -\frac{\text{i}}{2} \frac{\omega_j}{|\kappa_j|} a_j \exp{( \text{i}(\kappa_j x - \omega_j t) + |\kappa_j| z )} + \text{c.c.} , \end{gather}$$

(3.2b)$$\begin{gather}\eta(x,t) = \sum_{j=1}^{3} \frac{a_j}{2} \exp{( \text{i}(\kappa_j x - \omega_j t) )} + \text{c.c.} , \end{gather}$$

(3.2c)$$\begin{gather} \omega_j^2=g |\kappa_j| + \frac{D}{\rho_w} |\kappa_j|\kappa_j^4, \end{gather}$$

where c.c. denotes complex conjugate. A special case of (3.2) is the double-frequency case, where $(\omega, \kappa )_j$ are identical for $j=1,2$ and $a_1=a_2=a_3/2$.

The triad resonance conditions for flexural-gravity waves, similar to those for capillary-gravity waves, are given by Marchenko (1999), which we summarise here for completeness. The dispersion curve has an inflection point because the gravitational and flexural restoring exhibit different functional dependences. As a result, flexural-gravity waves admit triad resonant interactions at $O(\varepsilon ^2)$ while, in the limit of either pure gravity or pure flexural waves, resonant interactions do not occur until $O(\varepsilon ^3)$ for the 2-D deep water problem. The general triad resonance conditions are as follows:

(3.3a,b)\begin{equation} \kappa_1 \pm \kappa_2 ={\pm} \kappa_3, \quad \omega_1 \pm \omega_2 ={\pm} \omega_3, \end{equation}

where each wave component satisfies (3.2c) and the same sign combinations must be taken in both equations. Given our focus on sum-frequency resonant interactions, we consider only the plus signs in (3.3a,b) and assume $|\kappa _1| \le |\kappa _2|$.

A particular case of general sum-frequency resonance which is of special interest is the double-frequency case, where $\kappa _1=\kappa _2\equiv \kappa _0$, hereafter referred to as the primary wave, and $\kappa _3=2\kappa _0$, the double-frequency wave, form a resonant triad. Using the dispersion relationship (3.2c), we obtain $\kappa _0=2{\rm \pi} /\lambda _0=( \rho _w g/(14D))^{1/4}$ and $\omega _0=2{\rm \pi} /T_0\simeq 0.744 (g^5 \rho _w/D)^{1/8}$, which are functions of the ice sheet physical parameters. For convenience, we specify a general triad (3.3a,b) relative to the double-frequency condition using the parameter $\gamma =(\kappa _2-\kappa _0)/\kappa _0$. Figure 2(a) displays the wavenumber combinations which form sum-frequency triads.

Figure 2. (a) Wavenumber combinations that form resonant triads with $\kappa _1$ (black), $\kappa _2$ (blue) and $\kappa _3$ (red) as a function of $\gamma =(\kappa _2-\kappa _0)/\kappa _0$; (b) phase speed $C_p$ (——–) and group speed $C_g$ (– – – ) as functions of frequency.

To obtain a uniformly valid solution when the resonance conditions are satisfied, we apply standard MS analysis (e.g. Mei, Stiassnie & Yue 2005) up to $O(\varepsilon ^2)$. We introduce slow variables $\bar {t}=\varepsilon t$ and $\bar {x}=\varepsilon x$ and obtain governing equations for $\varPhi = \varPhi (x,\bar {x},z,t,\bar {t})$ and $\eta = \eta (x,\bar {x},t,\bar {t})$. For the general triad resonance case, the amplitude evolution equations are

(3.4a)$$\begin{gather} \frac{\partial a_1}{\partial \bar{t}} + C_{g1} \frac{\partial a_1}{\partial \bar{x}} =-\text{i} \frac{\omega_1 \omega_2}{2|C_{p1}|} a_2^*a_3 , \end{gather}$$

(3.4b)$$\begin{gather}\frac{\partial a_2}{\partial \bar{t}} + C_{g2} \frac{\partial a_2}{\partial \bar{x}} =-\text{i} \frac{\omega_1 \omega_2}{2|C_{p2}|} a_1^*a_3 , \end{gather}$$

(3.4c)$$\begin{gather}\frac{\partial a_3}{\partial \bar{t}} + C_{g3} \frac{\partial a_3}{\partial \bar{x}} =-\text{i} \frac{\omega_1 \omega_2}{2 |C_{p3}|} a_1 a_2, \end{gather}$$

where an asterisk indicates complex conjugate; here $C_{pj}=\omega _j/\kappa _j$ is the phase speed and $C_{gj}= \omega _j/(2\kappa _j) + 2 D\kappa _j^3 |\kappa _j|/(\rho _w \omega _j)$ the group speed, for $j=1,2,3$, both shown for a range of frequencies in figure 2(b). We remark that for the double-frequency case, $C_{p1,2}=C_{p3}=C_{p0}$, while $C_{g1,2}< C_{p0}< C_{g3}$.

It is known (e.g. Simmons (1969) for capillary-gravity waves) that the double-frequency case is not recovered in the limit as $\omega _1 \rightarrow \omega _2$ in (3.4), but the double-frequency evolution equations can be derived separately to obtain

(3.5a)$$\begin{align} \frac{\partial a_1}{\partial \bar{t}} + C_{g1} \frac{\partial a_1}{\partial \bar{x}} &=-\text{i}\frac{\omega_1 |\kappa_1|}{2} a_1^* a_3, \end{align}$$

(3.5b)$$\begin{align}\frac{\partial a_3}{\partial \bar{t}} + C_{g3} \frac{\partial a_3}{\partial \bar{x}} &=-\text{i}\frac{\omega_1 |\kappa_1|}{4} a_1^2. \end{align}$$

Equations (3.4) and (3.5) show that for resonant triads, significant interactions occur at a propagation distance of $L_{\varepsilon }=O(\lambda _0/\varepsilon )$. It is instructive to see how this compares with the length scale $L_{\nu }$ at which viscous dissipation beneath the ice sheet must be considered. In terms of the double-frequency parameters, the dissipation length scale is (see Liu & Mollo-Christensen 1988)

(3.6)\begin{equation} L_{\nu} = \frac{C_{g0}}{\kappa_0\sqrt{\nu_w \omega_0}} , \end{equation}

where $\nu _w$ denotes the eddy viscosity used to parametrise the turbulent boundary layer. Viscous dissipation effects can be neglected relative to nonlinear triad interactions for $L_\varepsilon /L_{\nu }< O(1)$. In terms of wave and ice parameters, this gives

(3.7)\begin{equation} \varepsilon>O(\sqrt{\nu_w \omega_0}/C_{g0}) . \end{equation}

For realistic sea ice and water parameters, table 1 gives values of (3.7) for different values of ice thickness $h$ and the corresponding double-frequency resonant wave period $T_0$. For all subsequent results (§ 5), the values of $\varepsilon$ satisfy (3.7).

Table 1. Threshold wave steepness values $\varepsilon$ given by (3.7) for different values of ice thickness $h$ and eddy viscosity (corresponding to different ice roughness) $\nu _w\,$$= 4\ \text {cm}^2\ \text {s}^{-1}$ (Liu, Holt & Vachon 1991) and $\nu _w=24\ \text {cm}^2\ \text {s}^{-1}$ (Hunkins 1966), for realistic physical parameters $E=6$ GPa, $\nu =0.3$, $\rho _w=1025\ \text {kg}\ \text {m}^{-3}$ and $g=9.81\ \text {m s}^{-2}$.

4. High-order spectral method for wave–ice interactions

We develop an $O(N)$ numerical solution for wave–ice interactions by extending the HOS method for nonlinear wave–wave interactions (Dommermuth & Yue 1987) through the inclusion of the dynamic boundary condition (2.1c) for ice. We use this method to validate the theoretical results for semi-infinite ice sheets and extend the analysis to more general wave–ice interaction problems. These involve multiple resonant interactions and higher-order effects, for which theoretical analysis becomes increasingly cumbersome, or ice sheets of finite length, where the application of theoretical analysis is not as straightforward given the feedback from the trailing edge. The HOS method is capable of modelling nonlinear interactions of a large number $N$ of spectral modes up to an arbitrary order $M$ in wave steepness. It can be shown through linear stability analysis (cf. Pan 2020) that, due to the flexural term, the maximum time step $\Delta t \sim 1/N^{5/2}$ for sufficiently large $N$, as compared to $\Delta t \sim 1/N^{1/2}$ for open water computations. Time integration is performed using the fourth-order Runge–Kutta method, and the time step is limited by stability rather than accuracy. The operation count is $O(M N \ln N)$ per time step and the method exhibits exponential convergence with respect to $N$.

For future reference, we denote wave amplitudes in open water by $A$, with superscripts $(\,)^{I,R,T}$ respectively indicating incident, reflected and transmitted waves, and denote the reflection and transmission coefficients by $R_j=A_j^R/A_1^I$ and $T_j=A_j^T/A_1^I$, with $j=1,2,3$, for a finite ice sheet. Since $R$ and $T$ capture the cumulative effects of reflection/transmission at the leading and trailing edges as well as the nonlinear amplitude evolution within the ice sheet, it is useful to consider the reflection/transmission at each edge separately. These single-edge reflection/transmission coefficients are equivalent to considering, at the leading edge, waves propagating from open water into a semi-infinite ice sheet or, for right-going waves at the trailing edge (or left-going waves at the leading edge), propagation from a semi-infinite ice sheet into open water (Meylan & Squire 1993). For wave amplitudes, $a$, on ice, we denote right-going and left-going waves with superscripts $+$ and $-$, respectively. We define $\mathcal {R}=(A^R/A^I)_{x=0}$ and $\mathcal {T}=(a^+/A^I)_{x=0}$ as the leading-edge reflection and transmission coefficients for waves propagating from open water to ice. We define the ice-to-water reflection and transmission coefficients as $\tilde {\mathcal {R}}=(a^+/a^-)_{x=0}=(a^-/a^+)_{x=L}$ and $\tilde {\mathcal {T}}=(A^R/a^-)_{x=0}=(A^T/a^+)_{x=L}$, respectively. For the semi-infinite case, we drop the $\pm$ since only right-going waves are present.

The HOS method has been used to study wave–ice interactions in the past (Bonnefoy et al. 2009; Guyenne & Părău 2012), with the majority of the applications considering ice sheets of infinite extent, and only more recently finite ice sheets (Guyenne & Părău 2017a,b; Xu & Guyenne 2023). Similar to these works, in the context of the HOS method, we model a finite ice sheet by smoothly varying the flexural rigidity term $(D(x)\eta _{xx})_{xx}$ in (2.1c) from $0$ to $D$ over a taper length $L_T$. While the numerical results are affected by values of the numerical parameter $L_T$, we remove this uncertainty by performing convergence tests so that the final results are insensitive to further decreases in (small) $L_T$. We do not explicitly impose the ice edge boundary conditions (2.1d), which we find does not generally alter the numerical results. Figure 3 shows a representative converged prediction of the HOS method for the single-edge reflection and transmission coefficients $\mathcal {R}$ and $\mathcal {T}$ as well as $\tilde {\mathcal {R}}$ and $\tilde {\mathcal {T}}$ for small wave steepness $\varepsilon _w$, compared with linearised numerical predictions (Fox & Squire 1990).

Figure 3. Single-edge reflection and transmission coefficients (a) from open water to ice and (b) from ice to open water, denoted by $\mathcal {T}$($\tilde {\mathcal {T}}$) (——–) and $\mathcal {R}$($\tilde {\mathcal {R}}$) (– – – ), obtained from eigenfunction matching (cf. Fox & Squire 1990), HOS ($M=1$) ($\circ$), and HOS ($M = 2$) for $\varepsilon _w = 0.05$ ($\varepsilon = 0.01$ for $\tilde {\mathcal {T}}$, $\tilde {\mathcal {R}}$) ($\times$) and for $\varepsilon _w = 0.1$ ($\varepsilon = 0.02$) ($\Delta$).

We perform systematic convergence studies for the computational parameters of the HOS method for semi-infinite and finite ice sheets for increasing $M$ and increasing number of spectral modes per open water wavelength $N/N_\lambda$. Table 2 shows nonlinear wave component amplitudes for a semi-infinite ice sheet, where $a_1^I = \mathcal {T}_1 A_1^I$, while table 3 displays reflection and transmission coefficients for a finite ice sheet. These results consistently show exponential convergence with respect to $M$ and $N/N_{\lambda }$. We further show convergence to steady-state reflection and transmission coefficients with increasing simulation time $T_S$ (starting from the time incident waves enter a finite ice sheet) in table 3.

Table 2. Convergence of wave component amplitudes at $x/\lambda _0=25$ with the number of spectral modes per open-water wavelength $N/N_{\lambda }$ and the nonlinear order $M$ for a semi-infinite ice sheet using the representative double-frequency case with $\varepsilon =0.04$.

Table 3. Convergence of reflection and transmission coefficients $R_j$ and $T_j$ for double-frequency triad wave components $j=1,3$ for a finite ice sheet of length $L/\lambda _0=20$ and $\varepsilon =0.02$ with (a) the number of spectral modes per open-water wavelength $N/N_{\lambda }$ and the nonlinear order $M$, and (b) increasing total simulation time $T_S/T_0$ with $M=3$.

Based on systematic convergence results such as these, for all the HOS simulations we use the following settings. For semi-infinite cases, we use $N/N_{\lambda }=16$, $L_T/\lambda _0=0.175$, computation domain size $N_{\lambda }=k_0\ell /(2{\rm \pi} )=256$ and integration time step $T_0/\Delta t=704$. For the finite ice sheet, the double-frequency reflection and transmission coefficients of interest require a somewhat higher grid resolution, and we use $N/N_{\lambda }=32$, $L_T/\lambda _0=0.102$, $T_0/\Delta t=3952$, $\ell /L=6.9$ and simulation time to reach steady-state $T_S>O(100T_0)$. With these computational parameters, all the results we present are expected to have converged to within 1 % for semi-infinite cases and less than 5 % for finite ice cases. As an additional check, our computations are shown to conserve total fluid volume to within $O(10^{-5})$ and total energy to within 1 % of their initial values.

5. Results

5.1. Semi-infinite ice sheet: strain increase due to sum-frequency interactions

We first consider the case of the semi-infinite ice sheet, where we obtain theoretical (MS) solutions, and compare and complement these with HOS computations. We consider ambient incident wave components satisfying the double-frequency resonance condition, the general triad conditions or some combination thereof, yielding multiple coupled resonances.

5.1.1. Closed-form solutions for the nonlinear strain

We consider the steady, spatially varying solution pertaining to both the double-frequency case (3.5) and a general triad case (3.4), and decompose the complex amplitude $a_j(\bar {x})=r_j \text {e}^{\text {i} \theta _j}$, where $r_j$ and $\theta _j$ are real and functions of $\bar {x}$. Consider (i) pure amplitude modulation (constant phase) where, for the double-frequency case, $2\theta _1-\theta _3=\pm {\rm \pi}/2$ and, for a general triad resonance case, $\theta _3-\theta _2-\theta _1=\pm {\rm \pi}/2$, and (ii) the situation where the shortest wave has zero initial amplitude, $r_1(x = 0)=r_1^I$ and $r_3(x = 0)=0$ for double frequency, and $r_1(x = 0)=r_1^I$, $r_2(x = 0)=r_2^I$ and $r_3(x = 0)=0$ for a triad. For the double-frequency case, the amplitude evolution in closed form is

(5.1a,b)\begin{align} r_1(x)= r_1^I \text{sech}\left(r_1^I |\kappa_1| \frac{\omega_1}{\sqrt{8 C_{g1} C_{g3}}} x \right);\quad r_3(x)= r_1^I \sqrt{\frac{C_{g1}}{2 C_{g3}}} \text{tanh} \left(r_1^I |\kappa_1| \frac{\omega_1}{\sqrt{8 C_{g1} C_{g3}}} x \right). \end{align}

For the general triad interaction, the amplitude evolution can be expressed using Jacobi elliptic functions:

(5.2a–c) \begin{equation} r_1(x)=r_1^I {\rm dn}(u|m), \quad r_2(x)=r_2^I {\rm cn}(u|m) , \quad r_3(x)=r_2^I \sqrt{\frac{C_{p2}C_{g2}}{C_{p3}C_{g3}}} {\rm sn}(u|m) , \end{equation}

where

(5.3)\begin{equation} u=r_1^I \frac{\omega_1 \omega_2}{\sqrt{C_{p2}C_{g2}C_{p3}C_{g3}}} x \end{equation}

and

(5.4)\begin{equation} m=\frac{(r_2^I)^2 C_{p2}C_{g2}}{(r_1^I)^2 C_{p1}C_{g1}} \end{equation}

is the ratio of the energy transfer rate by the two incident wave components, and, by definition, $m=1$ for the double-frequency case. When $m=1$, each component transfers energy to the short wave at an equal rate such that the solution behaviour is monotonic rather than periodic (McGoldrick 1965).

A main interest here is the (maximum) bending strain $S(x,t)$ in the ice sheet, which governs the breakup of the sheet. For the Euler–Bernoulli beam, this is given by

(5.5)\begin{equation} S(x,t) = \frac{h}{2} \frac{\eta_{xx}}{(1+\eta_x^2)^{3/2}} = \frac{h}{2} \eta_{xx} + O(\varepsilon^3).\end{equation}

We define the ratio $Q_{T,D}(x)$ of the slowly varying strain envelope $S_{env}(x)$ relative to its value at the ice edge $S_{env}^I=S_{env}(0^+)$ (which is approximately equivalent to the linear solution). For the general triad ($T$) ($\,j=1,2,3$) and double-frequency ($D$) ($\,j=1,3$) interactions,

(5.6)\begin{equation} Q_{T,D}(x) = \frac{S_{env}(x)}{S_{env}^I} = \sum_{j=1} q_j(x), \quad \text{where}\ q_j(x)=\frac{\kappa_j^2r_j(x)}{\sum_{i=1}^{2} \kappa_i^2r_i^I} . \end{equation}

For the double-frequency interaction, we obtain the maximum of the envelope to be

(5.7)\begin{equation} \hat{Q}_D={\rm max}_x{Q_D(x)} = \sqrt{1+8 C_{g1}/C_{g3}} \approx 2.06, \end{equation}

where the ratio of the group speeds is a constant when the double-frequency condition is satisfied.

The location of this maximum is an intermediate point where both $r_1$ and $r_3$ are non-zero,

(5.8)\begin{equation} \frac{x}{\lambda_0} = \frac{\sqrt{8 C_{g1}C_{g3}}}{\varepsilon \omega_0 \lambda_0} \ln \left(4\sqrt{\frac{C_{g1}}{2C_{g3}}} + \sqrt{8\frac{C_{g1}}{C_{g3}}+1} \right) \approx \frac{0.605}{\varepsilon}. \end{equation}

Additional parameters become important when considering general triad interactions, namely, the ratio of the initial wave amplitudes, encompassed in the parameter $m$, and the deviation of the triad from the double-frequency resonance condition, $\gamma$. We remark that $\hat {Q}_T={\rm max}_x{Q_T(x)}$ is obtained in the monotonic limit of (5.2a–c) where $m=1$ and that $\hat {Q}_T$ decreases from approximately $1.6$ near $\gamma =0$ to an asymptotic value of $\hat {Q}_T\approx 1.4$ for large $\gamma$, reflecting the difference in conversion efficiency between the double-frequency and triad interaction cases.

5.1.2. Amplitude and strain evolution due to double-frequency interactions

We study the double-frequency interaction with incident wave frequency $\omega _0$ for two initial wave steepnesses ($\varepsilon =0.04,0.08$) using perturbation analysis and HOS computations with $M = 2,3$. Figure 4(a,c) show the steady-state amplitude evolution for the two $\varepsilon$ values considered. At the leading ice edge, most of the primary incident wave is transmitted with single-edge transmission coefficient $|\mathcal {T}_1|=0.83$. However, since an ice edge acts as a low-pass filter, much less of the incident double-frequency wave (if present) would be transmitted ($|\mathcal {T}_3|=0.25$). Nonlinear sum-frequency interactions provide a mechanism for introducing significant high-frequency wave energy in the ice sheet that is minimally transmitted past the leading edge. Considering the amplitude evolution in the ice sheet, the MS predictions and HOS results agree up to a nonlinear interaction length $\varepsilon x/\lambda _0\leq O(1)$ but deviate for greater distances, where the HOS results exhibit recurrence/oscillatory features in contrast to the monotonic behaviour of the analytic solution (5.1a,b). This is expected as the HOS method contains many (non-resonant) interacting modes, free and locked, not represented in the perturbation theory, and it is in principle valid (up to the included order $M$) for long times/distances beyond $\bar {x}=O(1)$ of the MS analysis. In this case, the presence of the locked wave component manifested in a non-zero $a_3^I$ in the HOS incident waves is found to be the main source of the deviation in figure 4. This latter point can be verified by adjusting the initial condition of the perturbation evolution equations (3.5) to account for $a_3^I$ (in the HOS method), thus obtaining oscillatory predictions in close qualitative agreement with the HOS results for longer distances.

Figure 4. Double-frequency interaction in semi-infinite ice sheet with (a,b) $\varepsilon =0.04$ and (c,d) $\varepsilon =0.08$. (a,c) Spatial amplitude evolution: MS (5.1a,b) (– – – ), HOS ($M=2$) (——–) and HOS ($M=3$) (– - – - –); (b,d) spatial strain evolution for HOS ($M=3$), where $a_j$, $q_j$ for $j=1$ (black) and 3 (red), $Q_D$ (green), $Q_M$ (cyan) and MS prediction ($\times$) from (5.7) and (5.8).

To show the effect in the HOS results when higher-order interactions are included, figure 4 plots results also for $M = 3$, which shows that the present double-frequency resonance mechanism is well captured at leading order $M = 2$.

The accompanying strain evolution from HOS computations at $M = 3$ is displayed in figure 4(b,d). It is clear that the strain from the $\omega _3$ component exceeds that of the $\omega _1$ component as energy is transferred to the higher wavenumber, causing an increase in $Q_D$. Both the sum of strain magnitudes $Q_D$ and the complex summation of strain from $N$ wave components $Q_M={\rm max}_t (\tfrac{1}{2}h \eta _{xx} )/S_{env}^I$ are plotted to show that the bulk of the total wave energy is captured in the two components of the nonlinear interaction and that discarding phase does not significantly alter the result. We consider only the sum of strain magnitudes (e.g. $Q_D$) for the remaining results. The MS predictions ((5.7) and (5.8)) are also plotted in figure 4(b,d), showing that MS provides a reasonable estimate of $\hat {Q}_D$ for this case.

In summary, we show that double-frequency interactions can result in approximately double the maximum strain, irrespective of initial wave steepness, relative to that predicted by linear theory. Furthermore, because this is a result of nonlinear resonant interactions, this maximum strain occurs at a distance $O(\lambda _0/\varepsilon )$ from the ice edge, in contrast to linear theory, which predicts that the maximum bending strain occurs within $O(\lambda )$ of the ice edge (Fox & Squire 1991).

The amplification of the maximum strain by an approximate factor of 2, predicted computationally (see e.g. figure 4b,d) and theoretically in (5.7), can be argued heuristically. From figure 4(a) for example (also cf. figure 4c), it is seen that there is some point in the interaction where most of the energy in the incident wave is transferred into the double-frequency $\kappa _3=2\kappa _0$ component. The amplitude of $a_3$ can be obtained by equating the energy flux associated with this component, $\mathcal {F}_3$, to that of the initial energy flux of the incident component, $\mathcal {F}_1^I$. For a flexural-gravity wave, $\mathcal {F}=\tfrac{1}{2} \rho _w (g+D\kappa ^4/\rho _w)a^2 C_g$, which gives $a_3\approx 0.45 a_1^I$. From (5.5), $S\sim h \kappa ^2 a$, so that ${\rm max}_x(S_3)/S_1^I\approx 2$.

To the authors’ knowledge, no experimental work specifically addressing nonlinear interactions among flexural-gravity waves has been conducted. However, the above results provide a plausible explanation for the ice breakup observed by Herman, Evers & Reimer (2018), who conducted wave tank tests for an ice sheet under varying (monochromatic) incident conditions. The intent of these experiments was not to study nonlinear effects, but, coincidentally, the frequency of the incident waves in some cases is within $1\,\%$ of $\omega _0$ (runs 2010–2050 of test group A). Using the ice and wave parameters matching the experimental conditions of run 2030, where ice breakup was observed, we calculate, using linear theory, a maximum strain of $S=2.9\times 10^{-3}$, which is below their reported critical breaking strain $S_C=5.3\times 10^{-3}$. However, (5.7) predicts a maximum strain of $S_{env}=5.9\times 10^{-3}$, which is above the breaking threshold. The spatial extent of ice breakup observed during the experiment reasonably aligns with the theory as well. Equation (5.8) predicts the maximum strain to occur at 24 m past the ice edge, which is within the region of observed ice breakup extending 34 m from the ice edge.

5.1.3. Amplitude and strain evolution due to general triad interactions

We now consider general triad resonance with a bichromatic incident wave consisting of components $\omega _1$ and $\omega _2$. Figure 5 plots representative results for $\kappa _2=2\kappa _1$ ($\kappa _3=3\kappa _1$, $\gamma =0.357$) and with (initial) amplitudes chosen to correspond to $m=1$. Because $m=1$, the transfer of energy from $a_1$ and $a_2$ occurs at a comparable rate (despite the difference in the initial steepnesses $\varepsilon _1$ and $\varepsilon _2$) and their evolutions are almost identical. As before, HOS and MS agree well for $\bar {x}\leq O(1)$. Figure 5(b) plots the corresponding strain evolution. Numerical integration of (3.4) yields that $\hat {Q}_T \approx 1.6$ for small $\gamma$ with $m=1$, which is corroborated reasonably well by direct HOS simulation. For broadband incident waves, mechanisms such as the above are expected to produce significant strain amplification over a broad range of frequencies.

Figure 5. Triad interaction in semi-infinite ice sheet where $\kappa _2=2\kappa _1$ and $\kappa _3=3\kappa _1$ ($\gamma =0.357$) with $\varepsilon _1=0.02$ and $\varepsilon _2=0.04$. (a) Spatial amplitude evolution: MS (5.2a–c) (– – – ) and HOS ($M=2$) (——–); (b) spatial strain evolution for HOS ($M=2$), where $a_j$, $q_j$ for $j=1$ (black), 2 (blue) and 3 (red), and $Q_T$ (green).

5.1.4. Features of multiple resonant interactions

So far, we have considered incident waves with specific frequency components that satisfy a triad resonance. Realistic broadband seas contain many wave components which may satisfy multiple (coupled) resonant triad conditions and, for greater interaction times/distances, lead to generalised multiple-resonance (MR) interactions. HOS is capable of capturing the evolution of these general MR interactions involving many interacting components (e.g. Alam, Liu & Yue 2009; Xiao et al. 2013). To elucidate the key mechanisms and features of such MR interactions, we consider the simplest case involving two interacting triads, which is shown to lead to a further increase in the total strain. Consider a semi-infinite ice sheet and representative bichromatic incident waves consisting of a primary wave of frequency $\omega _1=\omega _0$ and a long wave of frequency $\omega _L$. The double-frequency triad component associated with $\omega _1$ is $\omega _3=2\omega _0$. The long-wave component is chosen such that $\omega _L$ and $\omega _3$ form another resonant triad ($\gamma =1$) that generates a new short wave $\omega _S ({>}\omega _3)$. See figure 6(a). As before, we perform HOS simulations (with $M=2$) to obtain steady-state results.

Figure 6. Multiple resonances in a semi-infinite ice sheet with initial steepness $\varepsilon _1=0.05$ and $\varepsilon _L\approx 0.01$. (a) Wavenumber combinations for multiple resonances; (b) spatial amplitude evolution for HOS ($M=2$); (c) strain evolution, where $a_j,q_j$ for $j=1$ (purple), $L$ (black), $3$ (blue) and $S$ (red), and $Q_{MR}$ (green).

Figure 6(b) shows that the initial evolution of $a_1$ and $a_3$, up to $x \approx 10\lambda _0$, is similar to that of § 5.1.2 (see figure 4a), with $a_1$ and $a_3$ respectively decreasing and increasing monotonically with propagation distance. As $x/\lambda _0$ increases and the magnitude of $a_3$ becomes appreciable, however, the $a_L \unicode{x2013} a_3 \unicode{x2013} a_S$ resonant triad comes into play, generating and growing the short $a_S$ component, with some decrease in $a_3$ and $a_L$.

We define the MR total strain amplification relative to linear predictions: $Q_{MR} =\sum _j \kappa _j^2 |a_j|/\sum _i \kappa _i^2 |a_i^I|$, where $j = 1,3,S,L$ and $i = 1,L$. Figure 6(c) plots $Q_{MR}$ as a function of evolution distance, showing that it attains $\hat {Q}_{MR}={\rm max}_x Q_{MR}\approx 2.73$ (at $x/\lambda _0 \approx 26.1$), which, as expected, is primarily a combination of the strain $q_3$ and $q_S$ due to the short-wave components. The MR maximum amplification $\hat {Q}_{MR}\approx 2.73$ in this case is appreciably greater than the theoretical limit of $\hat {Q}_D\approx 2.06$ (cf. (5.7)) and $\hat {Q}_T\approx 1.57$ from (5.6) (also for $\gamma =1$). This can be expected for general MR interactions involving broadband incidence waves as multiple resonant interaction mechanisms are obtained which could transfer energy from lower-frequency to higher sum-frequency resonant components.

5.2. Finite ice sheet: strain and reflection/transmission modified by nonlinear interactions

We consider nonlinear wave propagation over a finite-length ice sheet. We show that in this case, sum-frequency interactions can increase the maximum strain to more than double the linear result and that the maximum strain generally occurs near the trailing edge of the ice sheet.

Wave transmission and reflection by a finite ice sheet are of special interest in the context of wave propagation through aggregate ice fields. In contrast to the semi-infinite case, finite ice sheet wave reflection and transmission are significantly affected by nonlinear interactions which transfer energy to frequencies at which the (trailing edge) reflection/transmission coefficients can be very different to those of the incident waves. Because of sum-frequency nonlinear interactions, the scattered waves surrounding the ice sheet are generally at higher frequencies and steeper than those predicted by linear theory, with a greater portion of the incident wave field reflected. In figure 7, we illustrate a double-frequency interaction in an ice sheet of moderate length ($L/\lambda _0=20 = O(\varepsilon ^{-1})$). The amplitude evolution of the right-going wave components is qualitatively similar to that in the semi-infinite case (e.g. figure 4a).

Figure 7. Double-frequency interaction in a finite ice sheet with $L = 20\lambda _0$ for $\varepsilon = 0.019$. (a) Spatial amplitude evolution with right-going waves $a_j^+$ (——–) and left-going waves $a_j^-$ ($\cdots \cdots \cdots$); (b) spatial strain evolution, where $a_j^{\pm }$, $q_j^{\pm }$ for $j = 1$ (black) and $3$ (red), and $Q_D'$ (green).

From linear theory, the single-edge reflection and transmission coefficients (see figure 3) are as follows: at $\omega _1 = \omega _0$, $|\mathcal {T}_1|=0.829$, $|\tilde {\mathcal {T}}_1|=1.21$ and $|\mathcal {R}_1| = |\tilde {\mathcal {R}}_1| = 0.034$; at $\omega _3 = 2\omega _0$, $|\mathcal {T}_3| = 0.246$, $|\tilde {\mathcal {T}}_3| = 3.53$ and $|\mathcal {R}_3| = |\tilde {\mathcal {R}}_3| = 0.364$. At the trailing edge ($x = L$), since $|\tilde {\mathcal {R}}_3|\gg |\tilde {\mathcal {R}}_1|$, the left-going $\omega _3$ component, $a_3^-$, is greater than the left-going $\omega _1$ component, $a_1^-$. While the left-going components satisfy the double-frequency resonance condition, they are of much lower wave steepness than the right-going components, so the nonlinear interaction is not significant over the length of the ice sheet. The reflected wave components effectively modify the initial conditions of the nonlinear interaction at $x = 0$ slightly, leading to the small initial decrease in the $a_3^+$ amplitude. We note that it can be shown that counter-propagating flexural-gravity waves do not resonantly interact at second order.

Compared with the semi-infinite case, the strain in a finite ice sheet is further amplified (see table 4) since $a_3^-$ significantly contributes to the total strain amplification (including the right-going ($+$) and left-going ($-$) waves) $Q_D' = \sum _j(q_j^+ + q_j^-)$ for $j = 1,3$, where $q_j^{\pm } = \kappa _j^2|a_j^{\pm }|/(\kappa _1^2 |a_1^I|)$. Unlike in the linear case where the strain is approximately uniform within the ice sheet, figure 7(b) shows the strain increasing with distance from the leading edge, with the maximum at the trailing edge.

Table 4. Finite ice sheet reflection and transmission coefficients and nonlinear strain amplification $\hat {Q}_D'$ corresponding to figure 7.

The wave field surrounding the ice sheet is modified due to the double-frequency interaction. Compared with the linear case, less wave energy is transmitted since a greater portion of energy transferred to $a_3^+$ is reflected, evident from table 4. Note that, as the $\omega _3$ component propagates from ice into open water, in both the transmitted and the reflected wave fields it satisfies the open-water dispersion relationship so that $k_3=4 k_1$. From table 4, the nonlinear double-frequency transmission and reflection coefficients are $O(1)$, implying that the double-frequency transmitted and reflected wave amplitudes are of comparable order of magnitude to the amplitude of the incident wave. Consequently, the steepnesses of the transmitted wave ($4k_1A_3^T=0.08$) and reflected wave ($4k_1A_3^R=0.03$) are much greater than those predicted by linear theory ($k_1A_1^T=0.02,k_1A_1^R=0.001$) for the case considered.

We remark that, while the current work considers only a single ice sheet, in the context of an aggregate ice field with multiple ice sheets, the energy transferred to the $\omega _3$ component is reflected to a greater extent than the $\omega _1$ component by each ice edge it encounters. The net result is that nonlinear interactions could appreciably increase wave attenuation due to scattering relative to the linear prediction.

According to linear theory, the reflection and transmission of a wave with frequency $\omega$ and wavenumber $\kappa$ over a finite ice sheet are determined by two parameters: $\omega /\omega _0$ and $\kappa L$. When double-frequency resonance conditions are satisfied, however, the triad interaction length scale parameter $\mathcal {L}_{NL} = L/L_\varepsilon = \varepsilon L/\lambda _0$ becomes significant. We demonstrate the dependence of the reflection and transmission coefficients and the maximum strain amplification $\hat {Q}_D'$ on $\mathcal {L}_{NL}$ by performing HOS simulations of the double-frequency case over a broad range of this parameter for different $\varepsilon$.

Figure 8(a,b) plot the primary and double-frequency transmission and reflection coefficients $T_1$, $R_1$, $T_3$ and $R_3$, as well as the maximum strain amplification $\hat {Q}_D'$, all as functions of $\mathcal {L}_{NL}$. The results for different values of $\varepsilon$ all collapse as a function of $\mathcal {L}_{NL}$, confirming that triad nonlinear interactions are the dominant underlying mechanisms. As expected, with increasing interaction distance, the transmission of the incident component $|T_1|$ decreases, while that of the resonant double-frequency component $|T_3|$ (as well as $|R_3|$) increases. Here, $|R_1|$ (figure 8a) is a combination of reflection from the leading edge ($|\mathcal {R}_1|$) and from the trailing edge ($|\tilde {\mathcal {R}}_1|$). From linear theory, figure 3, $|\mathcal {R}_1| = |\tilde {\mathcal {R}}_1| \approx 0.034$. With increasing $\mathcal {L}_{NL}$, the contribution from the leading edge is fixed while that from the trailing edge decreases with $a_1$ at $x = L$. The net result is an almost constant $|R_1|$ with a slight decrease for greater $\mathcal {L}_{NL}$. Similar to the semi-infinite case (e.g. figure 4b), the total strain amplification $Q_D'$ increases with interaction distance. For finite $L$, $\hat {Q}_D'$ is further enhanced by a factor (comparing figures 7b and 4b) associated with the strong reflection at the trailing edge.

Figure 8. Dependence on $\mathcal {L}_{NL}$ for (a) $|T_1|$ and $|R_1|$, (b) $|T_3|$ and $|R_3|$, and (c) $\hat {Q}_D'$. Results are for HOS ($M=3$) with $\varepsilon = 0.012$ ($\circ$), $\varepsilon = 0.015$ ($\square$) and $\varepsilon = 0.019$ ($\Delta$), and SPA (5.9) for $|T_j|$ ($\hat {Q}_D'$ in panel c) (——–) and $|R_j|$ (– – – ).

Because the reflection of waves propagating from ice to open water is somewhat weak ($|\tilde {\mathcal {R}}_1| = 0.034$, $|\tilde {\mathcal {R}}_3| = 0.364$), the multiple reflection results in figure 8 can be approximated, say, by using a single-pass (one right-going pass and one left-going pass with subsequent reflections ignored) approximation (SPA). This gives

(5.9a)$$\begin{gather} |T_1|=| \tilde{\mathcal{T}}_1 \mathcal{T}_1{a_1^+(L)}/{a_1^I}|,\quad |R_1|=|\tilde{\mathcal{R}}_1\tilde{\mathcal{T}}_1 \mathcal{T}_1{a_1^+(L)}/{a_1^I}|+|\mathcal{R}_1|, \end{gather}$$

(5.9b)$$\begin{gather}|T_3|=|\tilde{\mathcal{T}}_3\mathcal{T}_1{a_3^+(L)}/{a_1^I}|, \quad |R_3|=|\tilde{\mathcal{R}}_3\tilde{\mathcal{T}}_3 \mathcal{T}_1{a_3^+(L)}/{a_1^I}|, \end{gather}$$

(5.9c)$$\begin{gather}\hat{Q}_D'=\left| \frac{a_1^+(L)}{a_1^I} \right| (1+|\tilde{\mathcal{R}}_1| ) + 4 \left| \frac{a_3^+(L)}{a_1^I} \right| (1+|\tilde{\mathcal{R}}_3|), \end{gather}$$

where $|a_j^+(L)|/|a_1^I|$ is determined using (5.1a,b). Comparing (5.9) with HOS results, in this case the error is shown to be generally $O(\tilde {\mathcal {R}}^2)$.

The maximum amplification of the strain, represented by $\hat {Q}'_D$ for the finite ice sheet in figure 8(c), can be estimated from (5.9c) using the semi-infinite ice analytic result (5.1a,b) for $|a_3^+(L)|$. This yields the SPA estimate ${\rm max}_{\mathcal {L}_{NL}} \hat {Q}_D' \approx 2.66$ which is obtained at $\mathcal {L}_{NL} \approx 0.717$. The error is less than 2 % compared to the HOS prediction.

Primarily because of the higher reflection coefficient of shorter waves (at the trailing edge), nonlinear sum-frequency interactions generally reduce the energy flux $\mathcal {F}_{NL}$ into open water in the lee of a finite ice sheet. The magnitude of this reduction is of interest for the attenuation of waves through an aggregate ice field. For double-frequency triad interactions, we can obtain a first-principle argument for the lower bound for the ratio of nonlinear to linear energy fluxes, $\mathcal {F}_{NL}/\mathcal {F}_{L}$, where $\mathcal {F}_{L}=\tfrac{1}{2} \rho g |A_1^I|^2 C_{gw1}$ is the maximum linear energy flux assuming perfect transmission. Using SPA (5.9), the total energy flux into open water is, to leading order,

(5.10)\begin{align} \mathcal{F}_{NL}\equiv \mathcal{F}_1 + \mathcal{F}_3 & = \tfrac{1}{2} \rho g ( |A_1^T|^2 C_{gw1} + |A_3^T|^2 C_{gw3}) \nonumber\\ & = \tfrac{1}{2} \rho g |A_1^I|^2 ( |T_1|^2 C_{gw1} + |T_3|^2 C_{gw3}) . \end{align}

For (sufficiently) long $L$, there will be some points $x = x_k$ near which $a_1(x_k) \approx 0$, and $a_3(x_k)$ is maximum with $a_{3,max}/a_1^I \approx 0.45$ (cf. § 5.1.2). The lower bound of (5.10) is obtained for ice sheet lengths close to any $x_k$, yielding, with $C_{gw3}=C_{gw1}/2$ and using asymptotic values of $|T_1| = 0$ and $|T_3| \approx 1.32$ for large $L$,

(5.11)\begin{equation} {\mathcal{F}_{NL}}/{\mathcal{F}_{L}}=|T_1|^2 + |T_3|^2 {C_{gw3}}/{C_{gw1}} \approx 0.867 . \end{equation}

6. Conclusion

We study sum-frequency nonlinear resonant wave–wave interactions which occur when surface waves propagate through semi-infinite and finite-length thin floating ice sheets. For semi-infinite ice sheets, using multiple-scale (MS) perturbation analysis, we obtain theoretical predictions of the spatial evolution for the resonantly interacting wave components. We modify a numerical approach based on the high-order spectral (HOS) method of Dommermuth & Yue (1987) which accounts for an arbitrary number $N$ of interacting modes and nonlinear order $M$. The theoretical (MS) and HOS results agree well in the domain where MS is valid, while HOS provides a powerful method for more general cases involving finite-length ice sheets, multiple resonant interactions and higher-order effects.

For the semi-infinite ice sheet, in addition to nonlinear amplitude evolutions of the interacting wave components, a main interest is the maximum ice bending strain which may lead to ice breakup. Sum-frequency triad interactions generally amplify the maximum strain relative to linear predictions by some factor $\hat {Q} >$1. For the double-frequency ($D$) special case, MS theory, confirmed by direct HOS simulations, predicts $\hat {Q}_D \approx 2.06$ at a distance $O(\lambda _0/\varepsilon )$ from the ice edge, where $\lambda _0$ is the incident primary wavelength. By considering the total wave energy flux, we provide a heuristic argument for $\hat {Q}_D \approx 2$. For the general triad ($T$) case, $1 < \hat {Q}_T < \hat {Q}_D \approx 2$ in general. Multiple (sum-frequency) resonances can increase $\hat {Q}$ considerably. This has implications for general broadband waves incident on an ice sheet. We provide an example of two interacting triads involving four wave components, with $\hat {Q}_{MR} \approx 2.73$.

When the ice sheet length $L$ is finite, we emphasise the significance of both strain amplification and wave energy transmission past the ice sheet, the latter being of special interest for wave propagation and attenuation through aggregate ice fields. For the double-frequency special case, the maximum strain amplification that is $\hat {Q}_D' \approx 2.7 > \hat {Q}_D \approx 2$ occurs at the trailing edge with $L = O(\lambda _0/\varepsilon )$. For $L \leq O(\lambda _0/\varepsilon )$, we offer a single-pass approximation (SPA) accounting primarily for right-going amplitude evolution and reflection from the trailing edge. Using SPA, we obtain an estimate of the maximum strain amplification that is in good agreement with $\hat {Q}_D' \approx 2.7$ from direct HOS computation. Due to sum-frequency resonant interactions, wave energy is transferred to higher-frequency components which generally have much greater reflection coefficients (at the trailing ice edge). Consequently, the total nonlinear wave energy flux $\mathcal {F}_{NL}$ into the lee of the ice sheet is generally less than the predicted linear wave flux $\mathcal {F}_{L}$. For the double-frequency special case, we use SPA to obtain a lower bound $\mathcal {F}_{NL}/\mathcal {F}_{L} \approx 0.87$, useful as an estimate of nonlinear effects on wave attenuation. We show that for a range of $L$ and $\varepsilon$, the nonlinear quantities $\hat {Q}_D'$ and the finite ice reflection and transmission coefficients each collapse to a single curve as a function of the nonlinear interaction parameter $\mathcal {L}_{NL} = \varepsilon L/\lambda _0$.

While the MS equations, (3.4) and (3.5), and HOS are applicable to difference-frequency triad interactions, we have focused on sum-frequency interactions, which are primarily responsible for nonlinear amplification of the maximum ice strain and a reduction in transmitted wave energy flux. Difference-frequency interactions which transfer energy to lower wavenumber components are not expected to increase the total strain. For wave transmission over a finite ice sheet, the property of transmission/reflection coefficients at smaller wavenumbers should allow such long-wave components to be transmitted with minimum reflection at the trailing ice edge. In this case, the heuristic arguments leading to (5.11) should still hold to provide an estimate for the wave energy flux due to nonlinear difference-frequency interactions.

As an initial effort, we have considered only the 2-D problem. For the three-dimensional problem, which can account for oblique incidence and general ice sheet geometry, we expect sum- and difference-frequency interactions to lead to behaviour such as longshore wave generation (e.g. Alam, Liu & Yue (2010) for nonlinear wave–bottom interactions) not present in this work.