3.1. Flexural-gravity triad resonance conditions for waves generated by a moving load

It is well known that the leading-order resonance for nonlinear interactions of flexural-gravity waves in deep water can occur at second order in wave steepness, involving three wave components. Triad resonant interactions occur among combinations of wave components with wavenumbers $\{k_1,k_2,k_3\}$ and corresponding frequencies $\{\omega _1, \omega _2, \omega _3\}$ that satisfy the conditions

(3.1a) \begin{equation} k_1 \pm k_2 \pm k_3 = 0 , \quad \omega _1 \pm \omega _2 \pm \omega _3 = 0 , \end{equation}

(3.1b) \begin{equation} \omega _j^2 = g k_j + D k_j^5/\rho , \quad j=1,2,3 , \end{equation}

with the same combination of signs taken throughout (3.1a). Without loss of generality, we consider combinations where $k_3\gt k_2\gt k_1\gt 0$ and $\omega _3\gt \omega _2\gt \omega _1\gt 0$ such that all wave components propagate in the direction of travel of the moving load.

Focusing on triad interactions that affect the response of the moving load, we reduce the set of all possible combinations $\{k_1,k_2,k_3\}$ satisfying (3.1) to the subset that contains $k_L$ and/or $k_T$ . The qualitative behaviour of these interactions depends on the sign of $C_{gj} - U$ for $j=1,2,3$ . Of the many possible interactions, the combination $\{k_1, k_2=k_m, k_3\}$ at $U=c_m$ is interesting because resonant energy transfer can affect the rate at which energy accumulates in the LP $k_m$ component since $U=C_{g2}$ . Cases also exist for $U\gt c_m$ , where an RP wave has a group speed equal to $U$ such that energy transferred to this wave cannot radiate away from the load. Despite their interesting features, we do not pursue further analysis of these cases, but instead focus on the double-frequency case as justified below.

In general, energy transfer through triad interactions can occur if at least two of the three wave components have non-zero amplitude. Yet we observe in (2.2) that only one wave component has non-zero amplitude in either the leading or trailing wave field. Therefore, an additional wave component, present in the ambient wave field or due to body oscillations, is necessary for nonlinear energy transfer in general. However, at the double-frequency condition, a special case where $k_1 = k_2 = k_T$ and $\omega _1 = \omega _2 = \omega _T$ , nonlinear energy transfer occurs without an additional wave since the trailing wave interacts with itself. We pursue further analysis of the double-frequency case in this work since a nonlinear interaction occurs with or without an ambient wave, and we do not address other general triad interactions that generally require an additional wave.

Figure 2.

Phase speed $c_j=U$ (solid line) and group speed $C_{gj}$ (dashed line) as functions of wavenumber $k_j$ for $j=L,T$ , denoting leading and trailing waves, respectively. This is overlaid with the double-frequency condition at $U=c_d$ , which intersects at $k_T=k_d$ and $k_L=2k_d$ .

The double-frequency condition is the main interest of this work for this reason, and because the RP waves coherently interfere with the LP waves. In this case, (3.1a) simplifies to $2k_1 = k_3$ and $2\omega _1 = \omega _3$ , and the waves satisfy the dispersion relation (3.1b) when $c_1 = c_3=c_d$ as shown in figure 2. The double-frequency condition is obtained at a specific speed $U = c_d \approx 1.44 ( D g^3/\rho )^{1/8} \approx 1.09 c_m$ with $\omega _1 = (2 \pi /T_d) \approx 0.74 (g^5 \rho /D)^{1/8}$ and, likewise, $k_1 = k_d = (2 \pi /\lambda _d) = (\rho g/(14 D))^{1/4}\approx 0.68 k_m$ . Significantly, $C_{g1} \approx 0.69 c_m \lt U$ and $C_{g3} \approx 1.70 c_m \gt U$ .

Before discussing the amplitude evolution, we briefly explain how the double-frequency interaction associated with a moving load differs from the well-known behaviour in the context of ambient waves. The difference occurs since the LP component $k_T = k_1$ ( $C_{g1}\lt U$ ) trails behind the load, while the RP $k_3$ component ( $C_{g3}\gt U$ ) radiates towards the load. As it passes the load, the RP wave with $k_3 = 2k_d$ coherently interferes with the LP leading wave with $k_L = 2k_d$ . The behaviour is further complicated when a small disturbance with wavenumber $k_1$ is present in addition to the LP waves. Such a disturbance could plausibly result from the wakes of other vehicles moving nearby (van der Sanden & Short Reference van der Sanden and Short2017), ambient waves propagating from open water, or body oscillations (Duffy Reference Duffy1991). The $k_1$ disturbance instigates a difference–frequency interaction in the leading wave field where the $k_3=2k_d$ wave radiating ahead of the load transfers energy to its half-harmonic $k_1 = k_d$ radiating back towards the load. This creates a feedback loop with the sum–frequency interaction in the trailing wave field.

3.2. Forced double-frequency amplitude evolution equations

To describe the exchange of energy between wave components, we derive the amplitude evolution equations, in the presence of a moving pressure forcing, following the standard method of multiple scales (cf. Mei et al. Reference Mei, Stiassnie and Yue2005) for the double-frequency case. We take perturbation expansions in powers of wave steepness $\varepsilon$ ( $\ll 1$ ). In terms of the long length scale $X=\varepsilon x$ and long time scale $T=\varepsilon t$ , the expansions of the pressure forcing, velocity potential and ice sheet elevation are $P(x) = \varepsilon\, p(x) + O(\varepsilon ^2)$ , $\Phi (x,z,t;X,T) = \varepsilon\, \Phi ^{(1)}(x,z,t;X,T) + \varepsilon ^2\, \Phi ^{(2)}(x,z,t;X,T) + O(\varepsilon ^3)$ and $\eta (x,t,X,T) = \varepsilon\, \eta ^{(1)}(x,t;X,T) + \varepsilon ^2\, \eta ^{(2)}(x,t;X,T) + O(\varepsilon ^3)$ .

The forced amplitude evolution equations describing the second-order, double-frequency interaction in the moving frame are obtained as

(3.2a) \begin{align} \frac {\partial a_1}{\partial T} + C^{m}_{g1}\, \frac {\partial a_1}{\partial X} &= -\text {i}\, \frac {c_d k_1^2}{2 }\, a_1^*a_3 - \frac {1}{\varepsilon }\,\frac {\mathrm {d} p}{\mathrm {d} x}\, \frac {\exp (-\text {i} k_1 x)}{\rho c_d k_1} , \end{align}

(3.2b) \begin{align} \frac {\partial a_3}{\partial T} + C^m_{g3}\, \frac {\partial a_3}{\partial X} &= -\text {i}\, \frac {c_d k_1^2}{4 }\, a_1^2 - \frac {1}{\varepsilon }\,\frac {\mathrm {d} p}{\mathrm {d} x}\, \frac {\exp (-\text {i} k_3 x)}{\rho c_d k_3} , \end{align}

where $a_1$ and $a_3$ are the complex amplitudes of the $k_1=k_d$ and $k_3=2k_d$ components, and an asterisk indicates complex conjugate. Here, $C^m_{g1} = C_{g1} - c_d \lt 0$ and $C^m_{g3} = C_{g3} - c_d \gt 0$ are the group speeds relative to the moving frame. For definiteness, we assume a symmetric forcing given by a Gaussian distribution, $p(x) = F\exp (-x^2/2w^2)/\sqrt {2 \pi w^2}$ , where $w$ measures the length of the forcing, and $F$ is the forcing magnitude per unit width.

Our interest is the case of a relatively narrow forcing, $w \ll \lambda _d/\varepsilon$ . In this case, in the manner of Akylas (Reference Akylas1987), we replace the pressure distribution terms in (3.2) with the limiting form as $\varepsilon \to 0$ , capturing the leading-order contribution

(3.3) \begin{equation} \frac {1}{\varepsilon }\,\frac {\mathrm {d} p}{\mathrm {d} x} \exp (-\text {i} k_j x) \to {2 \pi }\, \mathcal {F}\left (\frac {\mathrm {d} p}{\mathrm {d} x}\right )(k_j)\, \delta (X) , \quad j=1,3 , \end{equation}

where $\mathcal {F} ({\mathrm {d} p}/{\mathrm {d} x} )(k)$ is the Fourier transform of $\mathrm {d} p/\mathrm {d} x$ :

(3.4) \begin{equation} \mathcal {F}\left (\frac {\mathrm {d} p}{\mathrm {d} x}\right )(k) = \frac {1}{2 \pi }\int _{-\infty }^\infty \frac {\mathrm {d} p}{\mathrm {d} x} \exp (-\text {i} k x)\,\mathrm {d}x . \end{equation}

This limiting form represents the effect of the pressure distribution collapsing to a point when viewed from the far field.

Using (3.3) in (3.2), we obtain

(3.5a) \begin{align} \frac {\partial a_1}{\partial T} + C^{m}_{g1}\, \frac {\partial a_1}{\partial X} &= -\text {i}\, \frac {c_d k_1^2}{2 }\, a_1^*a_3 - \text {i}\, \frac {F }{\rho c_d }\exp \left (-\frac {1}{2}\, w^2 k_1^2\right )\,\delta (X) , \end{align}

(3.5b) \begin{align} \frac {\partial a_3}{\partial T} + C^m_{g3} \frac {\partial a_3}{\partial X} &= -\text {i} \frac {c_d k_1^2}{4 } a_1^2 - \text {i} \frac {F }{\rho c_d }\exp \left (-\frac {1}{2}\, w^2 k_3^2\right )\,\delta (X) .\end{align}

3.3. Steady amplitude evolution

We consider steady solutions of (3.5) relative to the frame $X=\varepsilon (x^{\prime}-Ut)$ moving with the pressure forcing. We solve separately for the solutions corresponding respectively to the trailing $X\lt 0$ and leading $X\gt 0$ waves, with complex amplitudes $a_j^-$ , $a_j^+$ , $j=1,3$ . These are described by

(3.6a) \begin{align} C^{m}_{g1}\, \frac {\mathrm {d} a_1^\pm }{\mathrm {d} X} &= -\text {i}\, \frac {c_d k_1^2}{2 } (a_1^\pm )^*a_3^\pm , \end{align}

(3.6b) \begin{align} C^m_{g3}\, \frac {\mathrm {d} a_3^\pm }{\mathrm {d} X} &= -\text {i}\, \frac {c_d k_1^2}{4 } (a_1^\pm )^2 , \end{align}

with boundary conditions $a_1^+(X\to \infty )=a_1^{+\infty }$ and $a_3^+(X=0^+) = a_3^{0^+}$ for $X\gt 0$ , and $a_1^-(X=0^-) = a_1^{0^-}$ and $a_3^-(X\to -\infty ) = a_3^{-\infty }$ for $X\lt 0$ .

Our goal is to determine closed-form solutions of (3.6) that decay to zero as $X\to \pm \infty$ in the domain where a nonlinear interaction occurs. We decompose $a_j^\pm = r_j^\pm \exp (\mathrm {i}\theta _j^\pm )$ into its real-valued magnitude $r_j^\pm$ and phase $\theta _j^\pm \in [- \pi , \pi ]$ . To aid in solving (3.6), we consider its two leading-order spatially invariant quantities, the Hamiltonian $\mathcal {H}$ and the wave energy flux $\mathcal {F}$ , which decay to zero as $X\to \pm \infty$ . We invoke the spatial invariance of $\mathcal {H}$ to show that $\theta _j^\pm$ are constant, which reduces (3.6) from four coupled equations in terms of $r_j^\pm$ and $\theta _j^\pm$ to two coupled equations in terms of $r_j^\pm$ . In a similar manner, we use conservation of  $\mathcal {F}$ to decouple the remaining two equations, which can then be solved easily.

We express the Hamiltonian of (3.6) $\mathcal {H}^\pm (X) = 2\mathrm {i}\,\mathrm{Re} \{a_1^\pm (X)\, a_1^\pm (X)\, (a_3^\pm (X))^* \}$ , where $\mathrm{Re}$ indicates the real part. Invoking invariance at $X\gtrless 0$ relative to points at $X\to \pm \infty$ in the far field respectively for $\mathcal {H}^\pm$ , we write

(3.7) \begin{equation} \mathcal {H}^\pm (X)= \mathcal {H}^\pm (X\to \pm \infty ) = 0 , \end{equation}

with the same signs taken throughout. Expressing in terms of magnitude and phase, and assuming $r_j^{\pm \infty } \ll r_j^\pm (X=O(\pm 1))$ , the right-hand side of (3.7) is zero, and we are left with the relative phase relationship $2\theta _1^\pm -\theta _3^\pm = \pm \pi /2 + 2 \pi n$ , $n\in \mathbb {Z}$ . This is well known to correspond to pure amplitude modulation with constant phases (e.g. McGoldrick (Reference McGoldrick1970) for capillary-gravity waves).

In a similar manner, we express the wave energy flux $\mathcal {F}^\pm _j(X) = \rho \omega _j c_d (r_j^\pm (X))^2 C_{gj}^m/2$ for $j=1,3$ , where the group speed is relative to the moving frame. Invoking wave energy flux invariance in a similar manner, we write

(3.8) \begin{equation} \mathcal {F}_1^\pm (X) + \mathcal {F}_3^\pm (X) = \mathcal {F}_1^\pm (X\to \pm \infty ) + \mathcal {F}_3^\pm (X\to \pm \infty ) = 0 , \end{equation}

where the same sign is taken throughout. We have assumed that $r_j^{\pm \infty } \ll r_j^\pm (X=O(\pm 1))$ so that the right-hand side of (3.8) is zero. This follows from the fact that for both leading and trailing waves, as $|X|$ increases, energy is converted from the wave component radiating away from $X=0$ to the wave component radiating towards $X=0$ until the remaining outgoing energy is negligible. In a typical scenario where wave propagation and energy flow are in the same direction, the sum of the wave energy flux of each component is conserved. However, in this case, energy flow and wave propagation are in the same direction for the $k_3$ component but in opposite directions relative to the moving frame for the $k_1$ component. Consequently, the energy flux of the $k_1$ component is equal and opposite to that of the $k_3$ component. We note that the right-hand sides of (3.7) and (3.8) are only approximately zero in the presence of a small ambient $k_1$ wave with amplitude $r_1^a$ as considered in § 5.4, but the subsequent results are accurate as long as $r_1^a \ll r_j^\pm (X=O(\pm 1))$ .

We use (3.8) to decouple the remaining two equations from (3.6) as

(3.9a) \begin{align} \frac {\mathrm {d} r_1^\pm }{\mathrm {d}X} &= -\frac {c_d k_1^2}{2C^{m}_{g1}}r_1^\pm r_3^\pm = \frac {c_d k_1^2}{\sqrt {-8C^{m}_{g1}C^m_{g3}}}(r_1^\pm )^2 \; , \end{align}

(3.9b) \begin{align} \frac {\mathrm {d} r_3^\pm }{\mathrm {d}X} &= -\frac {c_d k_1^2}{4C^m_{g3}}(r_1^\pm )^2 = \frac {c_d k_1^2}{2C^{m}_{g1}}(r_3^\pm )^2 \; , \end{align}

with boundary conditions $r_1^-(X=0^-) = r_1^{0^-}$ for $X\lt 0$ , and $r_3^+(X=0^+) = r_3^{0^+}$ for $X\gt 0$ , and we consider $r_j^\pm \gt 0$ without loss of generality. The solution of (3.9) is obtained to be

(3.10a) \begin{align} r_1^-(x) &= {r_1^{0^-}}/{(1-\alpha r_1^{0^-} x)} , \quad r_3^-(x) = \chi\, r_1^-(x), \quad X\lt 0 , \end{align}

(3.10b) \begin{align} r_3^+(x) &= r_3^{0^+}/(1-\beta r_3^{0^+}x) , \quad r_1^+(x) = {r_3^+(x)}/{\chi }, \quad X\gt 0 , \end{align}

where $\alpha = c_dk_1^2/\sqrt {-8C^{m}_{g1}C^m_{g3}}\gt 0$ , $\beta = c_d k_1^2/(2C^{m}_{g1})\lt 0$ and $\chi = ({-C^{m}_{g1}/(2C^m_{g3})})^{1/2}\gt 0$ . We remark that the temporally invariant spatial amplitude evolution described by (3.10) is similar to the spatially-invariant temporal evolution associated with the three-wave explosive instability (Segur Reference Segur2008). While the explosive instability results in ‘blow up’ in finite time, this is not the case for (3.10) since the singularity is located outside the solution domain. Similar spatial amplitude evolution is known to arise in nonlinear optics in the context of second-harmonic generation in a metamaterial where the indices of refraction of the fundamental ( $k_d$ ) and second-harmonic ( $2k_d$ ) components are positive and negative, respectively (Popov & Shalaev Reference Popov and Shalaev2006).

3.4. Connecting the leading and trailing wave solutions

To determine the remaining unknowns $r_1^{0^-}$ and $r_3^{0^+}$ , (3.10a) and (3.10b) are connected by jump conditions at $X=0$ resulting from the delta functions in (3.5) as

(3.11a) \begin{align} a_1^{0^+} - a_1^{0^-} &= \hat {r}_1 \exp (-\textrm {i} \hat {\theta }_1) , \end{align}

(3.11b) \begin{align} a_3^{0^+} - a_3^{0^-} &= \hat {r}_3 \exp (\textrm {i} \hat {\theta }_3) , \end{align}

where $\hat {r}_j = |F| \exp (-w^2k_j^2/2)/(\rho c_d |C_{gj}^m|)$ and $\hat {\theta }_j = -\mathrm {sgn}(F)\, \pi /2$ , $j=1,3$ , denote the real amplitude and phase, respectively, of the LP wave components, which are equivalent to the linear solution. For a concentrated line load ( $w=0$ ) moving with speed $U=c_d$ , (3.11) agrees with (2.2), the linear prediction of Schulkes & Sneyd (Reference Schulkes and Sneyd1988) for far-field amplitudes and phases. We express (3.11) in terms of magnitude and phase:

(3.12a) \begin{align} &r_1^{0^+}\cos (\theta _1^{0^+}) - r_1^{0^-}\cos (\theta _1^{0^-}) = 0 , \quad r_1^{0^+}\sin (\theta _1^{0^+}) - r_1^{0^-}\sin (\theta _1^{0^-}) = \mathrm {sgn}(F)\, \hat {r}_1(w) , \end{align}

(3.12b) \begin{align} &r_3^{0^+}\cos (\theta _3^{0^+}) - r_3^{0^-}\cos (\theta _3^{0^-}) = 0 , \quad r_3^{0^+}\sin (\theta _3^{0^+}) - r_3^{0^-}\sin (\theta _3^{0^-}) = -\mathrm {sgn}(F)\, \hat {r}_3(w) . \end{align}

Even with the restriction to double-frequency triad interactions, there are different possible cases of interest determined by the sign and length $w$ of the forcing, and the possible presence of an ambient wave at wavenumber $k_1$ .

In general, a sum–frequency interaction always occurs in the trailing wave field. This is clear from (3.9b) (middle term), where the nonlinear term proportional to $(r_1^-)^2$ is non-zero due to the LP trailing wave. On the other hand, difference–frequency interaction in the leading wave field occurs only in the presence of a co-propagating ambient $k_1$ wave. This is evident from (3.9a) (middle term) where the nonlinear term proportional to $r_1^+r_3^+$ is non-zero if both the LP leading wave $r_3^+$ and an ambient $k_1$ component $r_1^a$ are non-zero.

In the absence of an ambient $k_1$ wave ( $r_1^a=0$ ), and hence a difference–frequency interaction, the leading wave amplitude evolution is not governed by (3.9) and is instead constant with $r_1^+=r_1^{0^+}=0$ and $r_3^+=r_3^{0^+}$ . Upon selecting $2\theta _1^- -\theta _3^- = \pi /2+2 \pi n$ with $n\in \mathbb {Z}$ , we can now determine from (3.10) and (3.12) that

(3.13) \begin{align} r_1^{0^-} = \hat {r}_1 , \quad \theta _1^{0^-} = -\mathrm {sgn}(F)\, \pi /2, & \quad r_3^{0^-} = \chi \hat {r}_1 , \quad \theta _3^{0^-} = \pi /2 , \nonumber \\ r_3^{0^+} = \hat {r}_3 - \mathrm {sgn}(F)\,\chi \hat {r}_1, & \quad \theta _3^{0^+} = -\mathrm {sgn}(F)\, \pi /2 . \end{align}

We observe that the sign of $F$ determines whether the LP leading wave with amplitude  $\hat {r}_3$ constructively or destructively interferes with the RP wave with amplitude $\chi \hat {r}_1$ . These cases can be further categorised based on the load length $w$ , which changes $\hat {r}_3/\hat {r}_1$ . For $F\gt 0$ , $r_3^{0^+}=0$ when $k_1w=k_1W=( \ln (-2C^{m}_{g1}/C^{m}_{g3} )/3)^{1/2}$ . We refer to $W$ as the critical load length since it corresponds to complete destructive interference in the leading wave field. For $w\gt W$ and $F\gt 0$ , (3.11) and (3.12) must be modified to be consistent with $r_j^\pm \gt 0$ , and further increase in $w$ results in an increase in $r_3^{0^+}$ .

We now consider the presence of a small ambient $k_1$ wave with amplitude $r_1^a \ll r_1^+$ , consistent with the assumption in (3.7) and (3.8) of vanishing far-field wave amplitudes. For simplicity, we restrict our analysis to the case where $w=0$ and $F\gt 0$ . Upon selecting $2\theta _1^- -\theta _3^- = \pi /2+2 \pi n$ and $2\theta _1^+ -\theta _3^+ = - \pi /2+2 \pi n$ with $n\in \mathbb {Z}$ , we find the following solution from (3.10) and (3.12):

(3.14) \begin{align} r_1^{0^-} = \left (\hat {r}_3/\chi - \hat {r}_1 \right )/2 ,\quad \theta _1^{0^-}= \pi /2 ,\quad & r_3^{0^-} = \left (\hat {r}_3 - \chi \hat {r}_1 \right )/2 ,\quad \theta _3^{0^-} = \pi /2 , \nonumber \\ r_1^{0^+} = \left (\hat {r}_3/\chi + \hat {r}_1 \right )/2 ,\quad \theta _1^{0^+} = \pi /2 ,\quad & r_3^{0^+} = \left (\hat {r}_3 + \chi \hat {r}_1 \right )/2 ,\quad \theta _3^{0^+} = - \pi /2 . \end{align}

5.1. Positive line load resulting in substantial reduction in wave drag

We consider the effect of the double-frequency interaction on the representative case of a concentrated line load moving on top of the ice sheet by considering $F\gt 0$ and $w=0$ (see figure 3(a) for a schematic), which exhibits a significant reduction in wave drag compared to linear theory.

Figure 4.

For $F\gt 0$ and $w=0$ , spatial variation of steady-state (a) amplitude and (b) phase: MS (3.10) and (3.13) (solid line) and HOS for $\varepsilon _1 = 0.05$ (dashed line) and $\varepsilon_1=0.1$ (dotted line), with $j=1$ (black) and $j=3$ (red). Also, HOS temporal evolution for $\varepsilon _1=0.05$ of (c) amplitude at $X=0^\pm$ , for $r_j^{0^+}$ ( solid line) and $r_j^{0^-}$ (dot-dashed line) with $j=1$ (black) and $j=3$ (red), and (d) wave drag, for $D_w^{\ell }$ (dot-dashed line) and $D_w^{\textit{n}\ell }$ (solid line).

Figures 4(a) and 4(b) show the steady-state amplitude and phase evolution as functions of the long spatial scale $\varepsilon _1 x/\lambda _d$ for a positive line load ( $F\gt 0$ , $w=0$ ) for two values of wave steepness ( $\varepsilon _1=0.05,0.1$ ). In the trailing wave field, $r_1^-(x=0^-) = \hat {r}_1$ in agreement with linear theory. Whereas linear theory predicts that $r_1^-(x\leqslant 0^-)= \hat {r}_1$ , figure 4(a) shows that $r_1^-$ decays in the $-x$ direction as it resonantly transfers energy to $r_3^-$ . The $k_3$ component grows as it radiates in the $+x$ direction towards $x=0$ . Figure 4(b) shows that both trailing wave components have constant phase, with the $k_3$ wave setting up out of phase with the $k_1$ wave. The HOS results for both the amplitude and phase evolution in the trailing wave field compare very well with the MS predictions over the range of wave steepness considered.

Upon reaching $x=0$ , the RP $k_3$ component from the trailing wave field with $r_3^-(x=0^-) = \chi \hat {r}_1$ destructively interferes with the LP $k_3$ component with amplitude $\hat {r}_3$ . Destructive interference occurs because the RP wave with $\theta _3^-= \pi /2$ is out of phase with the LP wave with $\hat {\theta }_3 = - \pi /2$ , as seen in figure 4(b). The resulting leading wave has amplitude $r_3^+ = r_3^{0^+}\approx \hat {r}_3 - \chi \hat {r}_1 \approx 0.15\hat {r}_3$ and phase $\theta _3^+ = \theta _3^{0^+} \approx - \pi /2$ . We observe that the leading wave amplitude and phase differ from the MS prediction as wave steepness increases. From figure 4(b), we see that $\theta _3^-(x=0^-)$ , the phase of the RP $k_3$ wave, agrees well with theory. However, $\theta _3^+$ has a constant offset from the MS prediction since $\hat {\theta }_3$ differs from theory as the $\varepsilon \to 0$ assumption of (3.3) becomes less applicable.

Figure 4(c) shows the temporal evolution of $r_{1,3}^{0^\pm }$ to a steady state that occurs over the nonlinear time scale $T_{NL}=T_d/\varepsilon _1=20T_d$ for $\varepsilon _1=0.05$ . We observe that $r_3^+$ initially grows before decaying simultaneously with the growth of $r_3^-$ , providing further evidence for destructive interference.

Figure 4(d) shows the temporal evolution of wave drag for the $\varepsilon _1=0.05$ case, including both linear ( $M=1$ ) and nonlinear ( $M=4$ ) simulations. We observe that wave drag from the linear simulation $D_w^\ell$ grows monotonically to its steady-state value, while the nonlinear wave drag $D_w^{\textit{n}\ell }$ grows to a peak before decreasing to its steady-state value. For the $\varepsilon _1 = 0.05$ case $\mathcal {D} = 0.02$ , and for $\varepsilon _1 = 0.1$ , $\mathcal {D} = 0.07$ . The $O(100)$ reduction in wave drag compared to the linear prediction can be explained by equating wave drag to outgoing wave energy flux relative to the moving frame. It is well known that $D_w = (-\mathcal {F}^- + \mathcal {F}^+)/U$ . For the linear case, $\mathcal {F}^- = \mathcal {F}_1^-$ with $r_1^- = \hat {r}_1$ , and $\mathcal {F}^+ = \mathcal {F}_3^+$ with $r_3^+ = \hat {r}_3$ . Accounting for the nonlinear interaction in the trailing wave field, $\mathcal {F}^- = \mathcal {F}_1^- + \mathcal {F}_3^- = 0$ from (3.8) and $\mathcal {F}^+ = \mathcal {F}_3^+$ where $r_3^+ \approx 0.15\hat {r}_3$ . Heuristically, since $\mathcal {F}_1^-\approx \mathcal {F}_3^+$ for the linear solution, $\mathcal {D} \approx (0.15\hat {r}_3)^2/(2\hat {r}_3^2) \approx 0.01$ , in agreement with the values measured with the near-field approach. For all the cases that we consider with non-zero wave drag, we find that the near-field measurement of $D_w^{\textit{n}\ell }$ agrees within $5\%$ with the far-field approach $D_w^{\textit{n}\ell } = \mathcal {F}_3^+/U$ with $r_3^+$ measured from HOS simulations. This suggests that the bulk of the wave drag reduction results from the double-frequency mechanism, as we describe.

For this case, the nonlinear interaction only slightly increases the maximum bending strain relative to the linear prediction. The maximum bending strain found from HOS simulation for $\varepsilon _1 = 0.1$ is $\max _x(S^{\textit{n}\ell }) = 1.1 \max _x(S^\ell )$ , which occurs at $x=0^-$ where $r_3^-$ is greatest.

The results above can be obtained in realistic conditions. To illustrate this, we provide a specific example from the experimental work of Squire et al. (Reference Squire, Robinson, Langhorne and Haskell1988) that was performed on sea ice floating on deep water with parameters $h=1.6\,\rm m$ and $D=1.8\times 10^9\,\rm Nm$ , a typical flexural strength of sea ice $\sigma _b = 1\,\rm MPa$ from Timco & Weeks (Reference Timco and Weeks2010), and standard constants $\rho =1025\,\rm kg\, m^{-3}$ , $g=9.81\,\rm m\,s^{-1}$ and $\nu = 0.3$ . With these parameters, the double-frequency interaction occurs for a load speed $U=74\,\rm km\,h^{-1}$ , which is within the range of speeds tested by Squire et al. (Reference Squire, Robinson, Langhorne and Haskell1988) and is plausible for normal vehicle travel over sea ice. This corresponds to a wavelength $\lambda _d = 250\,\rm m$ . We then consider the maximum load and corresponding maximum wave steepness at which the ice fails in bending to provide an upper bound on the physical relevance of the double-frequency interaction. We find that the maximum bending stress exceeds $\sigma _b$ for $F\gt 0.3\,\rm MN\, m^{-1}$ and $\varepsilon _1\gt 0.05$ . While representing three-dimensional (3-D) point loads as 2-D line loads offers only an approximation, a typical cargo truck with mass $M=22000\,\rm kg$ and width $b=2\,\rm m$ exerts $F=0.07\,\rm MN\, m^{-1}$ , or a typical heavy-lift cargo plane with mass $3.5\times 10^5\,\rm kg$ and width $6\,\rm m$ exerts $F=0.5\,\rm MN\, m^{-1}$ . Therefore, significant nonlinear effects are plausible in realistic settings. However, the highest steepness $\varepsilon _1=0.1$ cases that we present are likely unrealistic since ice break-up would occur. We include these high steepness cases to clearly show the effect of increased nonlinearity.

Figure 5.

For $F\lt 0$ and $w=0$ , spatial variation of steady-state (a) amplitude and (b) phase: MS (3.10) and (3.13) (solid line) and HOS for $\varepsilon _1 = 0.025$ (dashed line) and $\varepsilon_1=0.05$ (dotted line), where $j=1$ (black) and $j=3$ (red). (c) Envelope of free-surface curvature ( $S\propto \eta _{xx}$ ) from linear (dot-dashed line) and nonlinear (solid line) HOS simulations for $\varepsilon _1 = 0.025$ .

We remark that Squire et al. (Reference Squire, Robinson, Langhorne and Haskell1988) do not report any observations deviating from linear theory at $U=c_d$ during their physical experiments. However, most of their experiments are performed with a small truck with mass $2100\,\rm kg$ and width $2\,\rm m$ such that $F=0.01\,\rm MN\, m^{-1}$ , and any nonlinear effects would likely be small. Furthermore, they take measurements at least 30 m from the centreline of the road, whereas the 2-D predictions that we provide are likely qualitatively accurate only in the vicinity of the centreline. While we have not found evidence of nonlinear effects at $U=c_d$ elsewhere in the experimental record (e.g. Takizawa Reference Takizawa1985), new experimental techniques that capture displacement fields (e.g. van der Sanden & Short (Reference van der Sanden and Short2017), using synthetic aperature radar) rather than point measurements should make this possible.

5.2. Negative line load resulting in significantly increased bending strain

By considering a negative line load ( $F\lt 0$ , $w=0$ ) (see figure 3(b) for a schematic), figure 5(a) shows that the steady-state amplitude evolution for $\varepsilon _1 = 0.025, 0.05$ is similar to the positive line load case in the trailing wave field. However, in the leading wave field, $r_3^+ \approx 2 \hat {r}_3$ rather than being reduced as before. Figure 5(b) shows that switching from $F\gt 0$ to $F\lt 0$ results in $\hat {\theta }_1= \pi /2$ , a shift by $ \pi $ , as expected since the linear solution behaves as an odd function of $F$ . However, $\theta _3^-= \pi /2$ , unchanged when the forcing changes sign since the leading-order nonlinear effects are quadratic and behave as an even function of $F$ . As a result, we observe constructive interference in the leading wave field with $r_3^+ = r_3^{0^+}= \hat {r}_3 + \chi \hat {r}_1 \approx 2\hat {r}_3$ since $\hat {\theta }_3 = \pi /2$ as well. The HOS results and MS predictions of amplitude evolution compare well both fore and aft of the load. However, the HOS results show a linear decrease in $\theta _3^+$ as a function of $\varepsilon _1 x/\lambda _d$ . This corresponds to a wavenumber downshift due to finite amplitude effects, not predicted by the leading-order MS analysis. This wavenumber downshift is consistent with what we observe at non-resonant supercritical load speeds, which is also reported by Bonnefoy et al. (Reference Bonnefoy, Meylan and Ferrant2009).

By nearly doubling the leading wave amplitude, the maximum bending strain induced by the moving load is also doubled compared to the linear prediction. Figure 5(c) shows that the maximum bending strain according to linear theory occurs in the leading wave field with $\max _x(S^\ell ) = S^+\approx 2.5 S^-$ as in the $F\gt 0$ case. Since the nonlinear interaction results in $r_3^+ \approx 2 \hat {r}_3$ , the strain induced in the leading wave field is nearly doubled, and $\max _x(S^{\textit{n}\ell }) \approx 2 \max _x(S^{\ell })$ , as seen in figure 5(c).

We measure $\mathcal {D}=1.97$ for $\varepsilon _1=0.025$ , and $\mathcal {D}=1.96$ for $\varepsilon _1=0.05$ . Using the same heuristic argument as before, $\mathcal {D} \approx (2 \hat {r}_3)^2/(2\hat {r}_3^2) \approx 2$ , supporting that the increase in the leading wave amplitude is enough to nearly double the wave drag compared to the linear prediction despite no outward flux in the trailing wave field.

5.3. Positive forcing with critical forcing length resulting in vanishing wave drag

While line load ( $w=0$ ) forcing is considered in the majority of the results shown, we present results for a finite-length forcing to demonstrate that similar wave drag reduction occurs for a range of $w\gt 0$ . We specifically highlight that wave drag can become vanishingly small at a critical forcing length $w=W$ (see figure 3(c) for a schematic). From our theoretical results for a representative symmetric Gaussian forcing distribution, we observe from (3.11) that the far-field phases of the LP waves are not affected by the forcing length. However, as $w$ increases, both $\hat {r}_1$ and $\hat {r}_3$ decrease, but $\hat {r}_3$ falls off at double the rate. We identify a transition point at $w=W$ where the amplitude of the RP $r_3^-$ and LP leading wave $\hat {r}_3$ are equal but out of phase with each other, leading to $r_3^+=0$ . As in the $w=0$ case, $r_1^-$ and $r_3^-$ approach zero as $x \to -\infty$ due to the resonant interaction. As a result, neither leading nor trailing waves radiate energy to the far field, leading to vanishing wave drag. We remark that for a typical truck length $L=22$ m corresponding to $w=L/2=11$ m, we have $w = 1.6W$ , such that the vanishing drag condition is plausible to achieve. Figure 6 shows that the HOS results agree well with the theoretical prediction in the trailing wave field for a range of wave steepness. In the leading wave field, we observe that $r_3^+$ decreases linearly as a function of $\varepsilon _1$ towards the theoretical prediction $r_3^+=0$ . We measure $\mathcal {D} =0.003$ , $0.003$ and $0.007$ for the $\varepsilon _1 = 0.025$ , $0.05$ and $0.1$ cases, respectively. These values are reduced by an order of magnitude from the $w=0$ cases.

We have performed a systematic study of the parameter space to determine the effects of a finite-length load on wave drag and maximum bending strain more generally. For $F\gt 0$ , we find that as $w$ increases from $0$ to $W$ , $D_w^{\textit{n}\ell }$ decreases to nearly zero. As $w$ increases from $W$ , $D_w^{\textit{n}\ell }$ increases because $\hat {r}_3\lt r_3^{0^-}$ and the destructive interference becomes less complete. Maximum bending strain $\max _x(S^{\textit{n}\ell })$ remains constant as a function of $w$ since it occurs in the trailing wave field and is unaffected by the changing interference pattern in the leading wave field. For $F\lt 0$ , both $D_w^{\textit{n}\ell }$ and $\max _x(S^{\textit{n}\ell })$ decrease monotonically with increasing $w$ because $r_3^+$ is reduced as $\hat {r}_3$ decreases. The critical forcing length $w=W$ bears no special significance for $F\lt 0$ .

Figure 6.

For $F\gt 0$ and $w=W$ , spatial variation of steady-state amplitude for $j=1$ (black) and $j=3$ (red): MS (3.10) and (3.13) (solid line) and HOS for $\varepsilon _1=0.025$ (dashed line), $\varepsilon _1=0.05$ (dot-dashed line), and $\varepsilon _1=0.1$ (dotted line).

5.4. Resonant interactions in both leading and trailing wave fields

The previous cases consider a nonlinear interaction occurring only in the trailing wave field, which redirects energy flow to the leading wave field with wavenumber $2k_d$ . However, this leading wave is known to be unstable in the presence of its sub-harmonic $k_d$ , which can result in significant energy transfer even if the amplitude of the ambient $k_d$ wave satisfies $r_1^a\ll r_3^+$ . We revisit the positive line load case ( $F\gt 0$ , $w=0$ ) for $\varepsilon _1 = 0.05,0.1$ described in § 5.1, but with the addition of a small ambient noise component with $k^a = k_d$ and $\varepsilon _1^a = k^a r_1^a = 0.01 \varepsilon _1$ (see figure 3(d) for a schematic). Figures 7(a) and 7(b) display that the steady-state amplitude and phase evolution differ markedly from the no-noise ( $r_1^a = 0$ ) case in figures 4(a) and 4(b).

Figure 7.

For $F\gt 0$ and $w=0$ with an ambient $k_d$ wave with $\varepsilon _1^a = 0.01\varepsilon _1$ , spatial variation of steady-state (a) amplitude and (b) phase: MS (3.10) and (3.14) (solid line) and HOS for $\varepsilon _1 = 0.05$ (dashed line) and $\varepsilon_1=0.1$ (dotted line) for $j=1$ (black) and $j=3$ (red). Also, HOS temporal evolution for $\varepsilon _1=0.05$ of (c) amplitude at $X=0^\pm$ , for $r_j^{0^+}$ (solid line) and $r_j^{0^-}$ (dot-dashed line) for $j=1$ (black) and $j=3$ (red), and (d) wave drag, for $D_w^{\ell }$ (dot-dashed line) and $D_w^{\textit{n}\ell }$ (solid line).

In figure 7(a), we observe that the wave component amplitudes in the trailing wave field are much smaller than in figure 4( a ), but still decay algebraically. In the leading wave field, wave amplitude decay of both the $k_1$ and $k_3$ components now occurs since the ambient $k_1$ wave instigates the difference–frequency interaction. The result is a locally confined wave packet containing mutually trapped $k_1$ and $k_3$ components, which is steady in the moving frame. We observe that HOS results for both wave steepness values agree well with the amplitude and phase evolution in the leading wave field predicted by MS analysis. However, we see in figure 7(b) that there is an apparent phase offset in the LP trailing wave due to finite amplitude effects similar to those described in § 5.1.

Figure 7(c) shows that the temporal evolution of the wave component amplitudes at $x=0^\pm$ reach a steady state at $t\approx 700 T_d$ . However, the time to reach a steady state is an order of magnitude longer than the case without ambient noise (see figure 4 c). The time to reach steady state can be decreased by increasing $\varepsilon _1^a$ . We also observe that $r_1^-$ decays to zero at $t\approx 500 T_d$ , at which point $\theta _1^-$ jumps by $ \pi $ , which is equivalent to a sign change in $r_1^-$ . This sign change leads to the steady-state solution by arresting the growth of $r_3^{0^+}$ .

Figure 7(d) shows the temporal evolution of wave drag. The evolution behaves similarly to the $r_1^a = 0$ case (see figure 4 d) until $t\approx 400 T_d$ , where $D_w^{\textit{n}\ell }$ briefly becomes negative before settling to a steady-state value $D_w^{\textit{n}\ell }=O(10^{-4})$ . A wave drag coefficient of nearly zero implies that the forcing is not performing work on the system. A localised wave packet exists despite zero power input because the triad interactions in the leading and trailing wave fields prevent energy from radiating to the far field.

Although we focus on the positive line load case, we have systematically examined the effect of a small $k_d$ ambient wave component on a range of $w$ with $F \gtrless 0$ . For $0 \leqslant w \lt W$ , we find steady-state solutions qualitatively similar to figure 7. At $w=W$ , $r_3^+=0$ so no difference–frequency interaction occurs in the leading wave field, and the solution resembles figure 6. For $w\gt W$ , we find that the solutions are generally unsteady given sufficient time. These cases exhibit $D_w^{\textit{n}\ell }$ oscillating in time and energy radiated to the far field. Similar behaviour is also observed for cases where $F\lt 0$ .