Nonlinear flexural-gravity waves for flows over bottom topography

Chuang Liang; Baoyu Ni; Yuriy Semenov

doi:10.1017/jfm.2025.10436

Nonlinear flexural-gravity waves for flows over bottom topography

Published online by Cambridge University Press: 22 September 2025

and

Chuang Liang: Affiliation:
College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, PR China
Baoyu Ni*: Affiliation:
College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, PR China
Yuriy Semenov*: Affiliation:
College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, PR China
*: Corresponding authors: Baoyu Ni, nibaoyu@hrbeu.edu.cn; Yuriy Semenov, yuriy.a.semenov@gmail.com
Corresponding authors: Baoyu Ni, nibaoyu@hrbeu.edu.cn; Yuriy Semenov, yuriy.a.semenov@gmail.com

Article contents

Abstract
Introduction
Formulation of the problem and analysis
Results and discussion
Conclusions
Funding
Declaration of interests
References

Rights & Permissions

Abstract

The interaction between the flow in a channel with multiple obstructions on the bottom and an elastic ice sheet covering the liquid is studied for the case of steady flow. The mathematical model employs velocity potential theory and fully accounts for the nonlinear boundary conditions at the ice/liquid interface and on the channel bottom. The integral hodograph method is used to derive the complex velocity potential of the flow, explicitly containing the velocity magnitude at the interface. This allows the boundary-value problem to be reduced to a system of nonlinear equations for the unknown velocity magnitude at the ice/liquid interface, which is solved using the collocation method. Case studies are carried out for a widened rectangular obstruction, whose width exceeds the wavelength of the interface, and for arrays of triangular ripples forming the undulating bottom shape. The influence of the bottom shape on the interface is investigated for three flow regimes: the subcritical regime, $F \lt F_{{cr}}$, for which the depth-based Froude number is less than the critical Froude number, and the interface perturbation decays upstream and downstream of the obstruction; the ice-supercritical and channel-subcritical regime, $F_{cr} \lt F \lt 1$, for which two waves of different wavelengths extend upstream and downstream to infinity; and the channel-supercritical regime, $F \gt 1$, for which the hydroelastic wave extends downstream to infinity. The results revealed a trapped interface wave above the rectangular obstruction and the ripple patch. The resonance behaviour of the interface over the undulating bottom occurs when the period of ripples approaches the wavelength of the ice/liquid interface.

JFM classification

Interfacial Flows (free surface): Interfacial Flows (free surface)Waves/Free-surface Flows: Channel flow Waves/Free-surface Flows: Wave-structure interactions

Information

Type: JFM Papers
Information: Journal of Fluid Mechanics , Volume 1019 , 25 September 2025, A40
DOI: https://doi.org/10.1017/jfm.2025.10436 [Opens in a new window]
Copyright: © The Author(s), 2025. Published by Cambridge University Press

1. Introduction

The classic problem of interaction between a moving fluid with an elastic ice sheet and the bottom topography is widely studied in mathematics and engineering due to its extensive physical applications (Kheisin Reference Kheisin1967; Korobkin, Părău & Vanden-Broeck Reference Korobkin, Părău and Vanden-Broeck2011; Squire Reference Squire2020). The interaction among ice cover, moving fluid and bottom topography plays a crucial role in river and oceanic flows, influencing ice stability, wave propagation and the operation of coastal structures. Recently, this subject has garnered renewed interest due to climate change and the accelerated melting of Arctic ice. The opening of new shipping routes and the expansion of resource exploration in these regions requires better understanding ice-fluid-structure interactions (IWSI) in Arctic waters (Squire et al. Reference Squire, Dugan, Wadhams, Rottier and Liu1995; Părău & Dias Reference Părău and Dias2002; Blyth, Părău & Vanden-Broeck Reference Blyth, Părău and Vanden-Broeck2011).

Most studies related to bottom topography phenomena rely on linear free-surface water wave theories (Davies Reference Davies1982; Mei Reference Mei1985; Chamberlain & Porter Reference Chamberlain and Porter1995; Alam, Liu & Yue Reference Alam, Liu and Yue2009; Korobkin et al. Reference Korobkin, Malenica and Khabakhpasheva2018). Linear theories enable the derivation of dispersion relations and determination of two critical wave speeds: one at the Froude number $F=1$ , associated with the finite-depth channel, and another at the critical Froude number $F=F_{{cr}}$ , corresponding to the minimum phase speed of flexural-gravity waves (Kheisin Reference Kheisin1963). The linear theories enable us to predict the qualitative behaviour of ice/liquid interfaces and evaluate transmission and reflection coefficients for water waves or an ice/liquid interface propagating over various bottom profiles, including right steps, backward steps and undulating shapes (Karmakar, Bhattacharjee & Sahoo Reference Karmakar, Bhattacharjee and Sahoo2010; Boral et al. Reference Boral, Meylan, Sahoo and Ni2023; Barman et al. Reference Barman, Chanda, Tsai and Mondale2024a , Reference Barman, Chanda, Tsai and Hsub ).

Linear models have also been applied to studying the response of ice sheets to external disturbances, such as loads moving along the ice surface or bodies moving beneath the ice. The time-dependent response of ice in shallow water was investigated by Sturova (Reference Sturova2009), while Porter & Porter (Reference Porter and Porter2004) examined the effects of varying ice thickness and bottom topography using a variational approach. Further studies have analysed the impact of moving loads in frozen channels and submerged moving bodies (Shishmarev et al. Reference Shishmarev, Khabakhpasheva and Korobkin2016, Reference Shishmarev, Khabakhpasheva and Korobkin2019; Xiong et al. Reference Xiong, Ni, Semenov and Korobkin2025) on the ice/liquid interface. In addition, linear theory has been used to explore hydroelastic interactions in the presence of ice leads (Xue et al. Reference Xue, Zeng, Ni, Korobkin and Khabakhpasheva2021). Ni et al. (Reference Ni, Han, Di and Xue2020) provided a comprehensive review of both theoretical and experimental studies on IWSI.

Although linear theories provide a fundamental framework for analysing wave propagation in ice-covered waters, they rely on small-amplitude approximations and assume linearised boundary conditions. One major limitation of linear wave theory is its inability to describe wave behaviour near the critical wave speed, or the Froude number $F=F_{{cr}}$ , at which the wavelengths of the elastic upstream and gravity downstream waves coincide (Semenov Reference Semenov2021), resulting in complex superposition akin to resonance. Near the critical value, the linear theory predicts infinite displacement of the interface, leading to singularities in the solution. This unphysical result indicates that nonlinear effects become dominant and must be accounted for to obtain a physically meaningful results.

In real-world scenarios, wave interactions are governed by nonlinear mechanisms such as wave steepening, energy transfer and modulation instability. Linear theories assume non-interacting wave components, which is valid only in the small-amplitude limit. However, for finite-amplitude waves, nonlinear effects such as wave–wave coupling and amplitude-dependent dispersion emerge, significantly altering the wave structure (Milewski, Vanden-Broeck & Wang Reference Milewski, Vanden-Broeck and Wang2011; Guyenne & Părău Reference Guyenne and Părău2012). These higher-order behaviours necessitate the use of nonlinear corrections or fully nonlinear models to accurately capture interface deformation and the associated ice bending phenomenon (Bonnefoy, Meylan & Ferrant Reference Bonnefoy, Meylan and Ferrant2009; Gao, Wang & Milewski Reference Gao, Wang and Milewski2019). A proper understanding of nonlinear effects is essential for improving predictive models of ice-covered waterways and evaluating the influence of submerged obstacles or bottom topography on the fluid motion.

The nonlinear study of IWSI began in the 1980s with the work of Forbes (Forbes Reference Forbes1986, Reference Forbes1988). He employed a complex function theory to solve the nonlinear boundary-value problem for the water flow and used the Kirchhoff–Love plate theory to model the elastic ice cover. His analysis predicted periodic waves in a channel with an obstruction. Later, Părău & Dias (Reference Părău and Dias2002) investigated nonlinear effects slightly below the critical wave speed, or $F\lt F_{{cr}}$ , deriving a nonlinear Schrödinger equation and identifying elevation and depression solitary waves with decaying oscillatory tails. Further advancement in the nonlinear theory was made by Page & Părău (Reference Page and Părău2014), who examined the steady hydraulic fall problem in the presence of an ice sheet and bottom geometry. They applied Cosserat theory of hyper elastic shells to model the ice sheet and employed a boundary-integral equation methods to solve problem for the liquid domain. Their solutions revealed a wave train trapped between two bottom obstructions. Most recently, Ni et al. (Reference Ni, Semenov, Khabakhpasheva, Părău and Korobkin2024) showed that a hydro-elastic wave may extend to infinity upstream and downstream at a Froude number in the range $F_{{cr}}\lt F\lt 1$ . For $F\gt 1$ the interface may exhibit a hydro-elastic wave extending to infinity downstream. Near the critical Froude number $F_{{cr}}$ , steady flow may cease to exist due to a resonance interference between the elastic and gravity waves.

In this paper, we consider steady-state nonlinear flexural-gravity waves generated by bottom topography in a finite-depth channel. We consider various bottom profiles, including rectangular obstructions and trenches, as well as arrays of periodic obstructions modelling undulating bottoms. The primary objective is to investigate previously unexplored nonlinear effects of the bottom shape on the ice/fluid interface. We employ velocity potential theory in conjunction with the integral hodograph method, which is based on an integral formula that determines a complex function from the values of its argument and its magnitude given on the real and imaginary axes of the first quadrant. Accordingly, we choose the first quadrant as the domain of the parameter variable. Using this integral formula, we derive the complex velocity potential, which includes the velocity magnitude at the interface in explicit form. This allows the boundary-value problem to be reduced to a system of nonlinear equations for the unknown velocity magnitude at the interface, which is solved using the collocation method.

Case studies are presented for a rectangular obstruction/trench which width exceeds the wavelength of the ice/liquid interface. Such a large width also enables us to model both the forward-facing and backward-facing bottom steps. For the array of triangular ripples forming the undulating bottom, the obtained results revealed the length of ripples for which the interface exhibits both the resonance and localised behaviour. The resonance occurs when the period of the ripples approaches the wavelength of the interface.

The structure of this paper is as follows. Section 2 presents the mathematical formulation of the problem, detailing the governing equations and boundary conditions. Section 3 presents results for different obstruction geometries and flow regimes, analysing the effects of key parameters such as the Froude number and obstruction height and width. Finally, § 4 summarises the main findings. The details of the method are presented in the appendices.

2. Formulation of the problem and analysis

We analyse a two-dimensional steady flow through a channel with complex bottom topography covered by a thin ice sheet. Far upstream and downstream, the bottom and interface are flat, and the flow is uniform, with a constant velocity $U$ and a fixed channel depth $H$ . The ice sheet has a constant thickness, denoted by $\varDelta$ . A Cartesian coordinate system $XY$ is introduced, with its origin positioned on the channel bed. The $X$ -axis is aligned with the flow direction, while the $Y$ -axis extends vertically upward. A schematic of this configuration is shown in figure 1(a).

Figure 1.

( $a$ ) Physical plane, ( $b$ ) parameter, or $\zeta$ -plane, and compex potential, or $w$ -plane.

The liquid is assumed to be inviscid and incompressible, subjected to downward gravitational acceleration $g$ . The flow is irrotational, allowing the use of potential flow theory. The channel bottom profile consists of straight segments $A_i A_{i+1}$ , each forming an angle $\beta _i$ with the $X$ -axis.

We introduce the complex velocity potential of the flow, $W(Z) = \varPhi (X,Y) + i\varPsi (X,Y)$ , where $\varPhi (X,Y)$ is the velocity potential and $\varPsi (X,Y)$ is the streamfunction, with $Z = X + iY$ . The governing boundary-value problem for the velocity potential is formulated as

(2.1)

\begin{equation} {\nabla} ^2 \varPhi = 0, \quad {\nabla} ^2 \varPsi = 0 ,\end{equation}

in the liquid domain; the boundary condition on the channel bottom $Y_b = Y_b(X)$ enforces impermeability

(2.2)

\begin{equation} \frac {\partial \varPhi }{\partial Y} = \frac {\partial \varPhi }{\partial X} \frac {\partial Y_b}{\partial X}, \quad \varPsi = 0; \end{equation}

the dynamic boundary condition on the ice–liquid interface $Y = Y(X)$ is given by Bernoulli’s equation

(2.3)

\begin{equation} \frac {\rho V^2}{2} + \rho g Y + p_{\textit{ice}} + p_{{ext}} = \frac {\rho U^2}{2} + \rho g H + p_\infty , \end{equation}

where $V = |\boldsymbol{\nabla }\varPhi |$ is the velocity magnitude, $p_{\textit{ice}}(X)$ is the pressure along the ice/liquid interface (representing both the elastic and hydrodynamic pressure) and $p_\infty = p_a + \rho _i g \varDelta$ is its value at infinity. Here, $p_a$ is atmospheric pressure and $\rho _i$ is the ice density. The external pressure $p_{{ext}}(X)$ is applied in the intervals $X_{P_2} \lt X \lt X_{P_1}$ and $X_{T_1} \lt X \lt X_{T_2}$ to ensure a waveless interface far upstream and downstream.

We seek a solution satisfying the asymptotic condition $Y(X)_{X \to -\infty } = H$ , where $H$ is the channel depth. Since the flow is steady, the streamfunction at the ice–liquid interface remains constant and equals the flow rate across the channel:

(2.4)

\begin{equation} \varPsi = \textit{UH}. \end{equation}

The far-field condition is

(2.5)

\begin{equation} \boldsymbol{\nabla }\varPhi \to U, \quad X \to -\infty , \quad 0 \leqslant Y \leqslant H. \end{equation}

The elastic sheet is modelled under the principal assumptions that it is thin and that both its inertia and stretching can be neglected. Three models are commonly used in the literature to describe the behaviour of the elastic sheet: the Euler–Bernoulli beam, the Kirchhoff–Love plate and the Cosserat shell model. Milewski & Wang (Reference Milewski and Wang2013) provided a comprehensive analysis of these models and their application to three-dimensional flexural-gravity (FG) waves. They found that, although all three models yield qualitatively similar results, there are substantial quantitative differences in the modelling of three-dimensional FG waves. The Cosserat model offers an advantage as it includes a nonlinear term absent in the other models. Ni et al. (Reference Ni, Semenov, Khabakhpasheva, Părău and Korobkin2024) compared the Kirchhoff–Love and Cosserat models in the context of two-dimensional flow in a channel, specifically analysing the critical bending moment at which an ice sheet breaks. They observed that the difference between the Cosserat theory and the Kirchhoff–Love plate model, where the cubic curvature term is omitted, is small due to the relatively low curvature of the ice sheet before failure begins. In the present study, we employ the Cosserat model of hyperelastic shells (Plotnikov & Toland Reference Plotnikov and Toland2011), which is better justified from a theoretical point of view (Milewski & Wang Reference Milewski and Wang2013)

(2.6)

\begin{equation} p_{\textit{ice}} - p_a = D'\left (\frac {d^2 \kappa }{{\rm d}s^2} + \frac {1}{2}\kappa ^3\right ), \end{equation}

where $\kappa$ is the curvature of the interface, $D' = (({E \varDelta ^3})/{12(1 - \nu ^2)})$ is the flexural rigidity of the ice sheet, $E$ is Young’s modulus and $\nu$ is Poisson’s ratio.

The interaction between the obstruction, the flow and the elastic sheet may lead to generation of the interface waves propagating upstream and downstream, for example at Froude number $F_{{cr}}\lt F\lt 1$ . To confine these waves within the computational domain, the external pressure is applied in the intervals $P_1 P_2$ and $T_1 T_2$ (see figure 1 $a$ ). These intervals are placed sufficiently far from the bottom irregularities to avoid interference with the solution. The applied external pressure is given by

(2.7)

\begin{equation} p_{{ext}} = C_d V \frac {{\rm d}V}{{\rm d}X}, \end{equation}

where $C_d$ is a coefficient characterising wave attenuation. This coefficient increases linearly from zero at $P_1$ to a maximum value $C_{{up}} \gt 0$ at $P_2$ , after which it remains constant. Similarly, $C_d$ increases from zero at $T_1$ to $C_{{dw}} \gt 0$ at $T_2$ , and then remains constant; for Froude $F\gt 1$ , the damping coefficient must be negative, $C_{{dw}} \lt 0$ . Within the interval $P_1 T_1$ , or in the solution domain, $C_d = 0$ . This wave attenuation technique was originally proposed by Semenov (Reference Semenov2021) for a similar problem, but with an infinite water depth.

To facilitate the solution, we introduce a non-dimensionalisation of the variables using velocity $U$ and channel depth $H$ as reference quantities. The dimensionless coordinates are defined as $x = X/H$ , $y = Y/H$ and $s = S/H$ , while the dimensionless ice thickness is $\delta = \varDelta /H$ . The bottom and interface profiles are given by $y_b(x) = Y_b(X)/H$ and $y(x) = Y(X)/H$ , respectively. Similarly, the velocity potential $\varPhi$ and streamfunction $\varPsi$ are normalised by the characteristic quantity $\textit{UH}$ , yielding the dimensionless variables $\phi = \varPhi /\textit{UH}$ and $\psi = \varPsi /\textit{UH}$ . Under this normalisation, the streamfunction takes the values $\psi = 0$ along the channel bottom and $\psi = 1$ at the ice–liquid interface.

The non-dimensionalised dynamic boundary condition (2.3) becomes

(2.8)

\begin{equation} v^2 = 1 - \frac {2(y - 1)}{F^2} - 2D \left (\frac {d^2 \kappa }{{\rm d}s^2} + \frac {1}{2} \kappa ^3\right ) - \frac {2C_d}{H} v \frac {{\rm d}v}{{\rm d}s}, \end{equation}

where

(2.9)

\begin{align} v=|\boldsymbol{\nabla }\phi |=V/U, \qquad E_b=\frac {D^\prime }{\rho gH^4 }, \qquad D=\frac {E_b}{ F^2}, \end{align}

and $\kappa$ represents the curvature of the interface. The Froude number, defined as

(2.10)

\begin{align} F = \frac {U}{\sqrt {gH}}, \end{align}

characterises the flow speed. Equation (2.8) expresses the velocity magnitude in terms of the wave elevation and its derivatives. These variables will be further related through the derived expression for the complex potential in the following sections.

2.1. Complex potential

The classical hodograph method was developed by Joukowskii (Reference Joukowskii1890) and Michell (Reference Michell1890). They introduced an auxiliary parameter plane, or $\zeta $ -plane, typically chosen as an upper half-plane. They then defined two functions: the complex potential $w(\zeta )$ and the function omega

(2.11)

\begin{equation} \omega (\zeta ) = -\ln \left (\frac {{\rm d}w}{v_0 {\rm d}z}\right ) = -\ln \left (\frac {v}{v_0}\right ) + i\beta , \end{equation}

where $v$ and $\beta$ are the velocity magnitude and direction, respectively, and $v_0$ is the velocity magnitude along the free surface, assumed to be constant. Once these functions are determined, the velocity and the flow region can be obtained in parametric form as follows:

(2.12)

\begin{equation} \frac {{\rm d}w}{{\rm d}z} = v_0\exp [-\omega (\zeta )], \quad z(\zeta ) = z_0 + \int _0^{\zeta } \frac {{\rm d}w}{{\rm d}\zeta '}/\frac {{\rm d}w}{{\rm d}z} {\rm d}\zeta ', \end{equation}

where the function $z(\zeta )$ is referred to as the mapping function.

The Joukowskii–Michell approach provides a powerful framework for solving free-surface flow problems involving bodies with polygonal geometries under the condition of constant velocity along the free surface or interface (i.e. in the absence of gravity, surface tension and other external forces). In such cases, the functions $\omega (\zeta )$ and $w(\zeta )$ define polygonal domains. The conformal mapping of the upper half-plane onto these polygonal domains can be obtained using the Schwarz–Christoffel integral to complete the solution.

However, significant complexity arises when the slope of the body contour varies or when the velocity magnitude on the free surface or interface is influenced by gravity, surface tension or other forces. In these cases, the boundary conditions for the complex velocity function differ: along the solid boundary, the velocity direction is dictated by the slope of the body, whereas on the free surface or interface, the velocity magnitude derived from Bernoulli’s equation is varying. This results in a mixed boundary-value problem for the complex velocity function.

In this paper, we adopt an integral formula (Semenov & Cummings Reference Semenov and Cummings2007; Semenov & Yoon Reference Semenov and Yoon2009) that allows the determination of a complex function from its argument along the real axis and its magnitude along the imaginary axis in the first quadrant. Accordingly, we choose the first quadrant as the region for the parameter variable $\zeta = \xi + i\eta$ instead of a half-plane, as illustrated in figure 1(b). In the physical plane $z = x + iy$ (figure 1 $a$ ), this parameter region corresponds to the liquid domain: the real axis represents the channel bottom, while the imaginary axis corresponds to the interface.

The conformal mapping theorem permits us to fix the positions of three points: $O(O')$ at $\zeta = 0$ , $C$ at $\zeta = 1$ and $B(B')$ at $\zeta = \infty$ , as shown in figure 1(b). However, the locations of the points $A_i$ , at $\zeta = a_i$ $i = 1, \ldots , n$ are unknown and must be found using additional physical constraints. The complex velocity function derived in Appendix B takes the form

(2.13)

\begin{equation} \frac {{\rm d}w}{{\rm d}z} = \prod _{i=1}^{n} \left ( \frac {a_i - \zeta }{a_i + \zeta } \right )^{\alpha _i / \pi } \exp \left [ -\frac {i}{\pi } \int _0^{\infty } \frac {d \ln v}{{\rm d}\eta } \ln \left ( \frac {i\eta - \zeta }{i\eta + \zeta } \right ) {\rm d}\eta \right ], \end{equation}

where $\alpha _i = \beta _i - \beta _{i-1}$ is the change in the velocity direction across the corner point $A_i$ , $i = 1, \ldots , n$ .

On the bottom of the channel, the streamfunction $\psi \equiv 0$ , and at the interface, $\psi \equiv 1$ , as follows from the boundary conditions (2.2) and (2.4), while the potential $\phi$ varies from $-\infty$ to $+\infty$ . Thus, the domain of the complex potential $w = \phi + i\psi$ is an infinite strip of unit width: $-\infty \lt \phi \lt \infty$ , $0 \leqslant \psi \leqslant 1$ . Due to the simplicity of the complex potential domain, $w$ , we can use a conformal mapping to directly express the complex potential $w$ as a function of the parameter variable $\zeta$

(2.14)

\begin{equation} w(\zeta ) = \frac {2}{\pi } \ln \zeta , \qquad \frac {{\rm d}w}{{\rm d}\zeta } = \frac {2}{\pi \zeta }. \end{equation}

The complex potential in (2.14) is a logarithmic function of $\zeta$ , implying that $\zeta$ depends exponentially on the complex potential $w = \phi + i\psi$ . The arc-length coordinates $s_b \sim \phi$ and $s \sim \phi$ correspond to the bottom and the interface, respectively. This means that $\zeta$ depends exponentially on $s$ , which causes computational difficulties when the length of the computational domain exceeds approximately $5H$ .

We can eliminate the logarithmic singularity by removing the parameter variables $\zeta$ , $\xi$ and $\eta$ from (2.13). Substituting the real part of (2.14) into (2.13), we obtain the complex velocity as a function of the complex potential $w$ , whose inverse function is the derivative of the mapping function $z = z(w)$ (see Appendix B)

(2.15)

\begin{equation} \frac {{\rm d}z}{{\rm d}w} = \prod _{i=1}^{n} \left ( \frac {a_i + e^{w'}}{a_i - e^{w'}} \right )^{\alpha _i/\pi } \exp \left [ \frac {i}{\pi } \int _{-\infty }^{\infty } \frac {d \ln v}{{\rm d}\phi '} \ln \left ( \frac {i e^{\phi '} - e^{w'}}{i e^{\phi '} + e^{w'}} \right ) {\rm d}\phi ' \right ], \end{equation}

where $w' = \pi w/2$ , $\phi ' = \pi \phi /2$ , $a_i = e^{\phi _{A_i}}$ and $\phi _{A_i}$ is the potential at the corner point $A_i$ . The kernel of the integral in (2.15) decays exponentially as the difference $|\phi ' - w'|$ increases. Integration of (2.15) along $-\infty \lt \phi \lt \infty$ with $\psi = 1$ in the $w$ -plane yields the interface $OB$ ; integration along the same interval with $\psi = 0$ yields the bottom surface.

2.2. System of equations for the unknowns $\phi _{A_i}$ and the function $v(\eta )$

The potentials $\phi _{A_i}$ are unknown and must be determined from the lengths of the segments $|A_{i-1} A_i| = s_{bi}$ , given by

(2.16)

\begin{equation} s_{bi} = \int _{\phi _{A_{i-1}}}^{\phi _{A_i}} \frac {{\rm d}s_b}{{\rm d}\phi '} \, {\rm d}\phi ', \end{equation}

where

(2.17)

\begin{equation} \frac {{\rm d}s_b}{{\rm d}\phi '} = \left |\frac {{\rm d}z}{{\rm d}w}\right |_{w = \phi '} \!= \prod _{i=1}^n \left | \frac {a_i + e^{\phi '}}{a_i - e^{\phi '}} \right |^{\alpha _i/\pi }\! \exp \!\left [ \!\frac {i}{\pi } \!\int _{-\infty }^{\infty } \frac {d\ln v}{{\rm d}\phi ''} \big( \pi - 2 \arctan \big(e^{\phi ''} - e^{\phi '}\big)\big) {\rm d}\phi '' \right ]\!. \end{equation}

The velocity magnitude along the interface is obtained using the dynamic boundary condition (2.8), which depends on the interface shape $y(s)$ , its curvature and higher-order derivatives of the curvature. The shape of the ice/liquid interface is determined by integrating the derivative of the mapping function (2.15) along the upper boundary of the strip in the $w$ -plane

(2.18)

\begin{equation} \begin{aligned} y(\phi ) &= y(-\phi ^*) + \operatorname {Im} \left \{ \int _{-\phi ^* + i}^{\phi + i} \left [ \prod _{i=1}^n \left ( \frac {a_i + e^{w'}}{a_i - e^{w'}} \right )^{\alpha _i/\pi } \right .\right . \\&\quad \left .\left . \times \exp \left ( \frac {i}{\pi } \int _{-\infty }^{\infty } \frac {d \ln v}{{\rm d}\phi '} \ln \left ( \frac {i e^{\phi '} - e^{w'}}{i e^{\phi '} + e^{w'}} \right ) {\rm d}\phi ' \right ) \right ] {\rm d}w' \right \}, \end{aligned} \end{equation}

where $y(-\phi ^*) = H$ . Here, $-\phi ^*$ and $\phi ^*$ are the lower and upper boundaries of the computational region. The channel in the physical plane is truncated, and the flow outside the computational region, $|\phi | \gt \phi ^*$ , is uniform.

The arc-length coordinate along the interface is defined by

(2.19)

\begin{equation} s(\phi ) = s_0 + \int _0^{\phi } \frac {{\rm d}\phi '}{v(\phi ')}, \end{equation}

where $s_0 = 0$ is chosen at $\phi = 0$ . Then, the curvature and its higher derivatives are obtained using a spline approximation of the interface $ y(s)$ .

The dynamic boundary condition (2.8) can be rewritten in terms of pressure coefficients. The hydrodynamic pressure coefficient is

(2.20)

\begin{equation} c_p[\phi , v(\phi )] = 1 - v^2(\phi ) - \frac {2 \left [ y(\phi , v(\phi )) - 1 \right ]}{F^2} - C_d[s(\phi )] v^2(\phi ) \frac {{\rm d}v}{{\rm d}\phi }, \end{equation}

and the pressure coefficient due to the bending of the ice sheet is

(2.21)

\begin{equation} c_p^{\textit{ice}}[\phi , v(\phi )] = 2D \left ( \frac {d^2 \kappa }{{\rm d}s^2} + \frac {1}{2} \kappa ^3 \right ). \end{equation}

By equating the pressure coefficients $c_p$ and $c_p^{\textit{ice}}$ at a set of nodes $\phi _m$ , $m = 1, \ldots , M$ , distributed along the interface, a system of nonlinear equations with respect to the unknown velocities $v_m$ is obtained. The details of the numerical method are presented in Appendix C.

2.3. Analysis of the dispersion relation

The dispersion relation can be obtained by assuming that the velocity magnitude along the ice/liquid interface behaves same way as for water waves with a free surface (see Appendix A)

(2.22)

\begin{equation} F^2 \frac {kH}{\tanh (kH)} = 1 + E_{b} (kH)^4. \end{equation}

Figure 2.

Wavenumber vs Froude number for different thicknesses of the ice sheet, $\delta$ .

Figure 2 illustrates the relationship between the dimensionless wavenumber $kH$ and the Froude number, obtained from the solution of (2.22) for various ice sheet thicknesses. It can be seen from figure 2 that, without an ice sheet ( $E_b = 0$ ), each Froude number $F \lt 1$ corresponds to a single wavenumber $kH$ , which tends to zero as $F \to 1$ (thick black line).

In the presence of an ice sheet, there exists a minimum, or critical, Froude number $F_{{cr}}$ for which (2.22) has only one real root, denoted $k_{{cr}}H$ . In the range $F_{{cr}} \lt F \lt 1$ , the dispersion relation yields two real roots for each Froude number, corresponding to two distinct wavenumbers. One, denoted $k_{{gr}}$ , is close to the wavenumber associated with the flow without ice sheet (solid thick line), and the other, denoted $k_{\textit{ice}}$ , arises due to the elastic ice sheet. For clarity, we refer to these two waves as the gravity wave – driven by gravity – and the elastic wave – arising from the elasticity of the ice sheet.

For Froude numbers $F \lt F_{{cr}}$ , (2.22) has a complex root $k = k_r + i k_{{im}}$ , where the real part $k_r$ characterises the wavenumber and the imaginary part $k_{{im}}$ represents wave attenuation. The real part is shown in figure 2 by the dashed line. It is observed that $k_r$ varies only slightly near the critical value $k_{{cr}}$ .

The elastic wave is shorter than the gravity wave since $k_{\textit{ice}} \gt k_{{gr}}$ . The range of Froude numbers for which the dispersion relation admits two real roots is referred to as the ice-supercritical and channel-subcritical flow regime. As the ice thickness increases, the wavenumber $k_{\textit{ice}}$ decreases, and the critical Froude number $F_{{cr}}$ approaches the critical value for the channel, $F_{{cr}} \to 1$ . As a result, the interval $F_{{cr}} \lt F \lt 1$ , in which both gravity and elastic waves can coexist, becomes narrower.

For $ F \gt 1$ , corresponding to a channel-supercritical flow, there exists a single root of the dispersion relation, which continues the branch associated with the elastic sheet. Since perturbations cannot propagate upstream when $ F \gt 1$ , the elastic wave propagates infinitely downstream. We note that for flows with a free surface (Yoon & Semenov Reference Yoon and Semenov2011), as well as for flows with a liquid surface covered by broken ice (Ni, Khabakhpasheva & Semenov Reference Ni, Khabakhpasheva and Semenov2023), the corresponding dispersion equations have no real roots for Froude numbers $ F \gt 1$ . This leads to qualitatively different behaviour between the case of an intact elastic ice cover and that of a free surface or covered by broken ice for $ F \gt 1$ , as will be demonstrated in the next section.

3. Results and discussion

This section presents numerical results for the interface response at different flow regimes and bottom geometries. The flow regimes include the subcritical, ice-supercritical and channel-subcritical and channel-supercritical regimes. The analysis includes two types of bottom geometries: a rectangular obstruction and a periodic array of triangular ripples, which are used to simulate undulating bottom profiles. The influence of geometrical feature of the bottom (e.g. height, width and periodicity) is systematically investigated for all flow regimes. For each configuration, the interface behaviour is analysed, including the occurrence of localised deformation, wave trapping and resonance phenomena. The resonance behaviour is of particular interest when the period of the bottom ripples approaches the natural wavelength of the interface wave, resulting in amplified responses. The numerical methodology is presented in Appendix C.

3.1. Subcritical flows

3.1.1. Rectangular obstruction

In order to model a rectangular obstruction, we set the corner angles as $\alpha _1 = \alpha _4 = 0.5\pi$ and $\alpha _2 = \alpha _3 = -0.5\pi$ . The parameters $\phi _{A_i}$ , $i = 1, \ldots , 4$ , are determined from the system of (2.16) using the height $a$ and width $b$ of the rectangle.

Figure 3.

Nonlinear (solid lines) and linear (dash-dotted lines) approximations of the ice/liquid interface, and the local Froude number above the rectangle (dashed lines), for various rectangle heights shown by the solid thick line at the bottom. (a) Froude number $ F = 0.5$ and obstruction heights $ a = 0.04$ (black), $ 0.06$ (red), $ 0.08$ (blue); (b) Froude number $ F = 0.4$ and obstruction heights $ a = 0.08$ (black), $ 0.15$ (red), $ 0.21$ (blue). The rectangle width is $ b = 12$ , and the ice thickness is $ \delta = 0.003$ .

Figure 3 illustrates the influence of the rectangle height at a fixed width $ b = 12$ on the interface shape for the ice thickness of $ \delta = 0.003$ in the subcritical flow regime. The corresponding critical Froude number is $ F_{{cr}} = 0.5748$ . In figure 3(a), results for $ F = 0.5$ and obstruction heights $ a = 0.04$ , $ 0.06$ and $ 0.08$ are shown in black, red and blue, respectively. The solid lines represent the results of the present nonlinear theory, while the dash-dotted lines correspond to the linear approximation, derived from the present solution using two fundamental properties of linear theory. The first property is that the linear-theory results converge to the nonlinear-theory results as the disturbance height (i.e. the rectangle height) tends to zero. The second property is that the perturbed quantities (such as velocity, pressure or surface shape) are proportional to the perturbation parameter. In this context, it means that the interface’s response varies linearly with the rectangle height; thus, if the rectangle height is doubled, all resulting quantities also double. Based on these properties, the linear approximation of the interface can be expressed as

(3.1)

\begin{equation} y_{{lin}}(a,x) = \left [ y(a^*,x) - 1 \right ] \frac {a}{a^*} + 1, \end{equation}

where $y(a^*,x)$ is the interface shape predicted by the present solution with the height of rectangle $ a^* = 0.005$ .

The rectangular obstruction induces two types of disturbances in the ice–liquid interface. The first arises from the increased flow velocity between the rectangle and the interface, which leads to a pressure drop that deforms the interface downward, forming a trough directly above the obstruction. The second type of disturbance involves wave excitation near the edges of the rectangle. It is not apparent at the small height $ a = 0.04$ (where the local Froude number $ F_{{loc}} \lt F_{{cr}}$ ) and, as a result, is not captured by the linear approximation. This leads to a large discrepancy between the nonlinear results and the linear prediction at the larger height $ a = 0.08$ .

The agreement between the linear approximation and the nonlinear results is relatively good for $ a = 0.04$ and $ a = 0.06$ : the interface perturbations near the edges of the rectangle decay rapidly, and the interface becomes nearly flat above the obstruction. However, at $ a = 0.08$ , which is close to the maximum height for which a steady subcritical solution can be obtained (since $ F_{{loc}} \approx F_{{cr}}$ ), the edge perturbations become more pronounced and decay more gradually toward the centre. The interface is qualitatively different from that predicted by the linear approximation, which fails to capture this second type of response, indicating stronger nonlinear ice/liquid/bottom interactions.

The local Froude number is defined as

(3.2)

\begin{equation} F_{{loc}} = F \frac {v(x)}{\sqrt {y(x)}}, \end{equation}

where $ v(x)$ is the velocity magnitude along the interface and $ y(x)$ is the vertical displacement of the interface. The corresponding local Froude number, shown by the dashed line in figure 3(a) (right axis), reveals that, for $a = 0.08$ , oscillations occur near the critical Froude number, beyond which completely different behaviour of the ice/liquid interface can be expected. Previous studies by Semenov (Reference Semenov2021) and Ni et al. (Reference Ni, Semenov, Khabakhpasheva, Părău and Korobkin2024) demonstrated that, as the system approaches the critical Froude number, the wavelengths of the upstream elastic wave and the downstream gravity wave converge. This convergence leads to strong nonlinear interactions between the two wave modes. Such interactions can prevent solution convergence, even for relatively small obstruction heights.

Similar results are shown in figure 3(b) for a Froude number of $ F = 0.4$ . Since the ratio $ F / F_{{cr}}$ is smaller than in the case shown in figure 3(a), solutions can be obtained for larger obstruction heights – up to $ a = 0.21$ , at which the local Froude number approaches the critical value. A comparison of the interface behaviour in figures 3(a) and 3(b) shows that increasing the obstruction height leads to stronger nonlinear effects, as evidenced by the growing discrepancy between the nonlinear results and the linear approximation.

Figure 4.

Effect of rectangle width on the ice/liquid interface at the onset of convergence of the solution for Froude number $F = 0.5$ and ice thickness $\delta = 0.003$ . (a) Interface profiles for varying rectangle widths and corresponding heights: $b = 12$ , $a = 0.085$ (black); $b = 8$ , $a = 0.085$ (red); $b = 4$ , $a = 0.09$ (blue); $b = 2$ , $a = 0.140$ (olive); $b = 1.5$ , $a = 0.095$ (wine); $b = 1.0$ , $a = 0.10$ (cyan); $b = 0.5$ , $a = 0.130$ (magenta); $b = 0.1$ , $a = 0.205$ (orange). (b) Close-up view: $b = 4.0$ , $a = 0.09$ (blue); $b = 3.0$ , $a = 0.09$ (pink); $b = 2.3$ , $a = 0.095$ (cyan); $b = 2.1$ , $a = 0.135$ (orange); $b = 2.0$ , $a = 0.140$ (olive); $b = 1.9$ , $a = 0.140$ (magenta); $b = 1.8$ , $a = 0.095$ (green).

Figure 4 illustrates how the width of the rectangular obstruction affects the interface shape at the maximum height corresponding to the onset of convergence of the steady flow solution. The rectangular obstruction is shown by solid lines at the bottom, while the resulting interface deformation is represented by solid lines of the same colour at the top.

The results reveal that increasing the length of the obstruction in the range $0.1 \lt b \lt 1.5$ leads to a reduction in the maximum possible height. However, near $b \approx 2$ , there is an anomalous increase which is discussed below. This trend suggests that shorter obstructions induce more localised effects, causing larger surface depressions due to the concentrated influence of the obstacle. However, once the obstruction length exceeds $b \gt 2$ , the maximum height of the rectangle reaches a plateau. This indicates that extending the obstruction beyond this point no longer significantly influences the maximal allowable height, suggesting that the system reaches a geometric saturation.

The wavelength above the rectangle is approximately $\lambda = 2\pi /(kH) \approx 1.5$ for $F = 0.5$ (as seen in figure 2). For cases with $b \lt 2$ , represented by the orange, magenta, cyan and wine-coloured curves, the interface exhibits a localised depression extending approximately one wavelength. This behaviour suggests that shorter obstructions primarily generate sub-wavelength perturbations. For obstruction widths exceeding $b = 2$ , more complex oscillatory patterns emerge. At $b = 2$ (olive line), the interface deformation corresponds to approximately one full wavelength. When $b = 4$ (red line), the deformation shows two full oscillations, indicating that the interface exhibits multiple wave cycles as the obstruction width increases beyond $b = 4$ .

An anomalously large height of the rectangle and corresponding interface deformation is observed for $b \approx 2$ (olive line), which is magnified in figure 4(b) along with additional cases. Within the narrow range $1.9 \lt b \lt 2.3$ , the interface undergoes a notable transformation – from a single trough at $b = 1.9$ (green line) to two troughs at $b = 2.3$ (cyan line). This deviation from the expected trend suggests the presence of resonance-like interactions between the rectangular obstruction and the ice/liquid interface, or the emergence of nonlinear effects that significantly alter the typical interface response.

Figure 5.

(a) Effect of trench depth on the interface shape for a trench width $b = 12$ , Froude number $F = 0.5$ and ice thickness $\delta = 0.003$ . Results are shown for varying trench depths $a$ : 0.2 (orange), 0.5 (blue), 1.0 (olive), 2.0 (magenta), 4.0 (red) and 10 (black). (b) Maximum deformation at infinite trench depth for different trench widths: $b = 12$ (black), 6 (red), 4 (magenta), 2 (olive) and 1 (blue).

For a trench located at the bottom, a steady flow solution can be computed for arbitrarily large trench depths. Figure 5(a) illustrates the interface shapes corresponding to a rectangular trench of length $b = 12$ for varying depths. The trench profiles are represented by thick solid lines at the bottom of the figure, while the resulting interface deformations are shown by the solid lines at the top.

As the trench depth increases, the flow velocity in the region between the trench bottom and the interface decreases, leading to a pressure increase exerted on the interface. This pressure rise causes the interface to deform upward, forming a heel directly above the trench. The deformation remains symmetric about the trench centreline. Within the trench, the local Froude number $F$ decreases, moving further away from the critical value $F_{{cr}} = 0.5748$ . As a result, the interface becomes nearly flat across most of the trench. When the trench depth satisfies $a \gt b/2$ , the interface approaches a limiting configuration. Beyond this threshold, further increases in trench depth do not significantly alter the interface shape, indicating saturation.

The effect of trench width on the interface configuration at the infinite depths is shown in figure 5(b). As the trench width approaches zero ( $b \to 0$ ), the interface remains close to its undisturbed position. However, as the trench width increases, the maximum height of interface deformation also increases, reaching its saturation near $b \approx 12$ . Further increases in trench width do not significantly raise the height of the heel.

3.1.2. Sawtooth patch on the bottom

In this subsection, we investigate the influence of a periodic bottom topography composed of an array of triangles on the interface shape. The results highlight the conditions under which wave trapping or amplification occur, and demonstrate how bottom periodicity can affect the interface.

We consider a patch on the bottom composed of identical ripples in the shape of triangles. We choose the number of triangles to be $m = 5$ ; then, the total number of corner points on the bottom is $n = 2m + 1$ . Each triangle has a height $a$ and a base $b$ , which also represents the period of the bottom profile. The corner angles are defined as follows:

(3.3)

\begin{equation} \alpha _1 = \alpha _n = \gamma , \quad \alpha _{2i} = -2\gamma \ (i = 1, \ldots , m), \quad \alpha _{2j+1} = 2\gamma , \ (j = 1, \ldots , m-1), \end{equation}

where $\gamma = \arctan (2a / b)$ is the angle at the base of the triangle. The parameters $\phi _{A_i}$ , $i = 1, \ldots , n$ , are determined from the system of (2.14), using the side length of the triangle.

Figure 6.

(a) Effect of the sawtooth bottom profile on the interface shape for Froude number $F = 0.5$ , ice thickness $\delta = 0.003$ and two different heights of the ripples: (a) $a = 0.125$ , (b) $a = 0.115$ . The ripple base is indicated in the figure.

Figure 6 illustrates the interaction between the ice/liquid interface and the sawtooth patch on the bottom. The interface is influenced by both the height and the base/wavelength of the ripples. The results correspond to a Froude number of $F = 0.5$ and an ice thickness of $\delta = 0.003$ , for which the real part of the complex wavenumber determines a subcritical wavelength of the decaying wave (see figure 2, vertical dashed line).

When the base of the ripples is shorter than the interface wavelength ( $b \lt \lambda = 1.60$ ) – as in the cases $b = 0.5$ and $b = 1.0$ shown in figure 6(a) – the interface responds to the array of ripples similarly to how it would respond to a single widened obstruction. However, as the bottom wavelength approaches $b = 1.25$ , the interface begins to react to each patch element individually, generating a wave with relatively small amplitude. This wave consists of four crests and five troughs, with the troughs positioned precisely above the peaks of the ripples.

As the base length of the triangular elements increases, the amplitude of the interface wave also increases, reaching a maximum at $b = 1.75$ , which is close to the interface wavelength $\lambda = 1.60$ . Beyond this point, as the bottom wave period continues to increase, the amplitude of the interface wave slightly decreases, as observed for $b = 2.25$ .

For larger base lengths ( $b \gt 2.25$ , figure 6 $b$ ), the interface troughs reach a saturation level, with the minimum $y$ -coordinate stabilising around $y = 0.93$ . In the range $2.25 \lt b \lt 4.25$ , the interface crests are absent, and the profile between the troughs remains nearly flat. At $b = 4.75$ , a small crest begins to emerge above the trough of the bottom profile. With further increases in ripple base, the interface response above each ripple is expected to resemble the response caused by an isolated ripple, as the influence of neighbouring elements diminishes.

3.2. Ice-supercritical–channel-subcritical flows, $F_{{cr}} \lt F \lt 1$

This subsection focuses on the flow regime where the Froude number lies between the critical value for the ice sheet, $F_{{cr}}$ , and the channel critical value, $F = 1$ . In this range, the dispersion relation admits two real roots $k_{{gr}}$ and $k_{\textit{ice}}$ corresponding to the coexistence of gravity and elastic waves, which are extended to infinity downstream and upstream, respectively. This dual-wave interaction leads to a rich variety of nonlinear behaviours, including wave resonance and localisation of the interface deformation. In order to satisfy the radiation condition in numerical computations, damping regions are introduced (the regions $P_1P_2$ and $T_1T_2$ in figure 1 a) in which the waves are attenuated. These damping zones are placed sufficiently far upstream and downstream from the rectangular obstruction so that their effect is excluded from the region of interest. In the damping regions, the downstream-propagating gravity wave and the upstream-propagating elastic wave gradually decay, ensuring that the flow approaches a uniform state at infinity. For further computational details, we refer the reader to Semenov (Reference Semenov2021) and Ni et al. (Reference Ni, Semenov, Khabakhpasheva, Părău and Korobkin2024), where the methodology was originally proposed.

The results presented in this section illustrate how the interface shape responds under varying bottom geometries and how the presence of both wave modes affects steady-state solution existence and structure.

3.2.1. Rectangular obstruction

The effect of the Froude number on the interface profile and the bending moment is shown in figures 7(a) and 7(b), respectively, for a rectangular obstruction with height $a = 0.03$ and width $b = 16$ . The Froude number decreases from top to bottom, approaching the critical value $F_{{cr}} = 0.5748$ .

Figure 7.

(a) Ice–water interface and (b) bending moment, $m=-D\kappa$ , along the interface for ice thickness $\delta = 0.003$ and various Froude numbers in the range $F_{{cr}} \lt F \lt 1$ , as indicated in the figure. The rectangular obstruction has height $a = 0.03$ and width $b = 16$ ; the critical Froude number is $F_{{cr}} = 0.5748$ .

It can be observed that the interface above the obstruction and further downstream exhibits a near-sinusoidal shape, with a longer wavelength above the obstruction than downstream. This behaviour is consistent with the fact that the local Froude number in the region above the rectangle is higher than that downstream, due to a locally reduced channel depth which results in increased flow velocity and, consequently, a higher local Froude number.

For $F = 0.75$ , the upstream-propagating wave generated by the rectangular obstruction is nearly imperceptible in the interface shape, indicating a very low amplitude. However, this wave is clearly visible in the bending moment distribution. The elastic wavelength $\lambda _{\textit{ice}}$ is significantly shorter than the gravity wavelength $\lambda _{{gr}}$ . Moreover, the bending moment exhibits a clear superposition of elastic and gravity waves both above the rectangle and downstream. For this case, the wavelength ratio is $\lambda _{{gr}}/\lambda _{\textit{ice}} = k_{\textit{ice}} / k_{{gr}} = 4.1$ (see figure 2). This superposition is clearly visible in figure 7(b) for $F = 0.75$ .

At $F = 0.70$ , the ratio $\lambda _{{gr}}/\lambda _{\textit{ice}} = 3.2$ is slightly reduced, allowing three full waves to be accommodated above the rectangular instead of two for $F = 0.75$ . The amplitude of the upstream elastic wave is lower than for $F = 0.75$ , and its contribution to the interface deformation above the rectangle is negligible. This observation will be further explained in the following discussion.

At $F = 0.65$ , the wavelength ratio decreases to $\lambda _{{gr}}/\lambda _{\textit{ice}} = 2.4$ . The superposition of elastic and gravity waves remains evident, particularly in the bending moment distribution (figure 7 $b$ ), which displays non-sinusoidal behaviour due to the larger contribution of the elastic wave. The elastic wave re-emerges at the peaks of the gravity wave.

For $F = 0.60$ , the wavelength ratio is $\lambda _{{gr}}/\lambda _{\textit{ice}} = 1.64$ . The interface profile in figure 7(a) clearly reflects the contribution of the gravity wave upstream and the influence of the elastic wave downstream. As the wavelengths and amplitudes of the two wave types become more similar, their interference becomes stronger.

The non-monotonic effect of the Froude number on the amplitudes of the elastic and gravity waves observed in figure 7 arises from their interaction with the widened rectangular obstruction. Figure 8 presents interface profiles for a Froude number $F = 0.7$ and various obstruction widths. As the width of the obstruction changes, the amplitude of the interface wave above the obstruction remains nearly constant, while the amplitudes of the upstream and downstream waves exhibit oscillatory behaviour. The blue and red lines in figure 8(a) correspond to cases where both the upstream elastic wave and the downstream gravity wave attain their maximum amplitudes. These profiles differ from a sinusoidal shape due to the superposition of gravity waves (with larger amplitude) and elastic waves (with smaller amplitude). Between these cases, there exists an intermediate width $b = 15.58$ for which both the upstream and downstream waves vanish. In this configuration, the wave becomes trapped above the obstruction and exhibits a nearly sinusoidal shape. This behaviour persists as the width of the rectangle is further increased or decreased, as shown by the solid lines in (b) for $ b = 13.5$ , (c) for $ b = 11.68$ and (d) for $ b = 9.5$ . The dash-dotted lines represent the linear approximation. At the left end of the rectangle, the agreement between the solid and dash-dotted lines is quite good, and the wave amplitudes are very similar. However, the discrepancy grows downstream due to a mismatch in the wavelengths of the linear and nonlinear predictions. As a result, at the right end of the rectangle, the linear and nonlinear waves are out of phase, which significantly affects the wave amplitude downstream. Thus, nonlinearity affects the wavelength of the interface above the obstruction, which in turn has a significant impact on the interface profile downstream. As the rectangle width increases, the phase difference downstream becomes more pronounced until it completes a full cycle, or the number of waves above the rectangle changes by one.

Figure 8.

Effect of rectangle width $ b$ on the interface shape for Froude number $ F = 0.7$ and rectangle height $ a = 0.03$ . (a) Nonlinear results for widths $ b_{{ab}} = 17.9$ (blue), $ b_{\textit{ag}} = 15.58$ (green) and $ b_{{ar}} = 13.5$ (red); (b), (c) and (d) show comparisons between the nonlinear results (solid lines) and the corresponding linear approximations for widths $ b_{{bb}} = 13.5$ (blue), $ b_{{bg}} = 11.68$ (green) and $ b_{{br}} = 9.5$ (red), respectively.

The difference $b_{\textit{ag}} - b_{{bg}} = 3.9$ is close to the wavelength of the interface wave above the obstruction. This wavelength is larger than the gravity wavelength $\lambda _{{gr}} = 3.12$ downstream of the rectangle due to the increased local Froude number above the obstruction. Consequently, the amplitudes of the upstream and downstream waves as a function of the obstruction width demonstrate the oscillatory dependence.

Figure 9.

Effect of rectangle width $ b$ on the free-surface shape for Froude number $ F = 0.5$ and rectangle height $ a = 0.1$ . (a) Nonlinear results for widths $ b_{{ab}} = 9.4$ (red), $ b_{\textit{ag}} = 8.8$ (green) and $ b_{{ar}} = 8.4$ (blue); (b) shows comparisons between the nonlinear results from ( $a$ ) (solid lines of the same colour) and the corresponding linear approximations (shot-dashed lines).

Figure 9(a) presents the same results as figure 8(a), but for the case of a free surface. In contrast to the ice–liquid interfaces shown in figure 8(a), the free-surface profiles appear purely sinusoidal. The effect of the rectangular obstruction’s width on the amplitude of the downstream free-surface wave is similar to that observed for the ice–liquid interface in figure 8(a). This suggests that the presence of the ice sheet has only a weak influence on the gravity-driven wave. Indeed, examining the dispersion relation in figure 2 reveals that, for a fixed Froude number, the wavenumbers of the FG wave and the pure gravity wave ( $\delta = 0$ ) are very close. Figure 9(b) shows a comparison between the nonlinear free-surface profile and its linear approximation. The nonlinear effects are significant, as both the amplitude and wavelength above the rectangular obstruction differ markedly from the linear prediction. Similar to the results in figure 8(b), the linear approximation fails to capture the increase in wavelength above the rectangle, which results in a phase shift between the linear and nonlinear waves at the downstream edge of the rectangle. Beyond the rectangle, this phase shift remains constant, as the wavelengths of the linear and nonlinear waves match in the downstream region.

Figure 10.

Ice/liquid interface for Froude number $F = 0.7$ , ice thickness $\delta = 0.003$ and ripple height $a = 0.015$ , shown for various ripple periods $b$ indicated in the figure. Resonance is observed at $b = 3.25$ , which is close to the FG wavelength $\lambda _{{gr}} = 3.12$ .

3.2.2. Sawtooth patch on the bottom

The objective of this section is to investigate the nonlinear resonance interactions involving FG waves propagating over a sawtooth bottom. Figure 10 presents the ice/liquid interface profiles for a Froude number $F = 0.7$ and ice thickness $\delta = 0.003$ , corresponding to various periods of the ripples in the range $0.5 \lt b \lt 4$ with the height $a = 0.015$ . This relatively small ripple height is chosen to obtain a converged solution close to the resonance bottom period, $b = 3.25$ , for which the interface amplitude reaches its maximum value. For smaller or larger bottom wavelengths, the limiting ripple height allowing for solution convergence varies within the range $0.015 \lt a \lt 0.07$ .

Figure 10 highlights the resonance behaviour of the interface, which occurs when the period of the bottom undulations approaches the wavenumber of the FG wave. Specifically, resonance is observed for the bottom period $b = 3.25$ , which is close to the wavelength of the ice/liquid interface, $\lambda _{{gr}} = 3.12$ . At this period, the amplitude of the interface increases gradually over the sawtooth patch and then remains constant downstream over the flat bottom.

For smaller bottom periods, the interface response resembles that of a single obstruction, with length equal to the length of the sawtooth patch, similar to that shown in figure 8. The interface response exhibits an oscillatory behaviour as the total patch length $l_p = M \boldsymbol{\cdot }b$ , with $M = 5$ , increases. The amplitude of the interface wave first reaches a maximum and then decreases. When the wave downstream of the patch disappears entirely, the interface over the patch accommodates an integer number of the wave cycles, as seen for $b = 0.66$ , $1.31$ , $1.97$ , $2.625$ and $3.925$ , corresponding to patch lengths $l_p = 3.3$ , $6.55$ , $9.85$ , $13.12$ and $19.62$ , respectively. The period of this sequence is approximately 3.27. The case $l_p = 16.25$ with $b = 3.25$ also fits into this sequence between $b = 2.625$ and $3.925$ . However, instead of leading to wave cancellation, this configuration results in resonance behaviour due to the close match between the bottom period and the wavelength of the FG wave $\lambda _{{gr}} = 3.12$ obtained from the dispersion relation.

3.3. Channel-supercritical flows, $F \gt 1$

It is well known from studies of free-surface channel flows (Dias & Vanden-Broeck Reference Dias and Vanden-Broeck1989) that, in the supercritical regime ( $F \gt 1$ ), two distinct types of steady solutions may exist: one with a smaller wave height, known as the perturbed wave, and another with a larger wave height, termed the soliton wave. Both the perturbed and soliton waves exhibit a wave-free free surface. In contrast to these classical free-surface solutions, the recent research by Ni et al. (Reference Ni, Semenov, Khabakhpasheva, Părău and Korobkin2024) has shown that a wavy ice/liquid interface may exist in a channel covered by an ice sheet, which extends to infinity downstream. The wave interface exhibits oscillations around the perturbed-type free-surface solution with the wavelength corresponding to that of the elastic wave at the supercritical regime ( $F \gt 1$ ) (see figure 2).

Two types of free-surface profiles are shown in figure 11 for supercritical channel flow with $F = 1.3$ , over a rectangular obstruction of height $a = 0.05$ and two widths: $b = 9.0$ (blue) and $b = 5.0$ (red). Attempts to find a soliton-type solution in the presence of an elastic ice sheet have been unsuccessful, in contrast to the cases of free-surface flows (Yoon & Semenov Reference Yoon and Semenov2011) and flows with a broken ice cover (Ni et al. Reference Ni, Khabakhpasheva and Semenov2023). In these earlier studies, both perturbed-type and soliton-type solutions were obtained for supercritical channel flows using the same theoretical framework and numerical methods. However, in the present case involving an elastic sheet, the same approach fails to produce a soliton-type solution for $F \gt 1$ . The key difference in the elastic sheet case lies in the nature of the dispersion relation: for $F \gt 1$ , it possesses a real root, whereas in the free-surface and broken ice scenarios, the dispersion relations have no real roots in the supercritical regime. Consequently, both perturbed and soliton-type solutions in those cases feature a waveless interface in the far field. At present, there is no mathematical proof confirming either the existence or non-existence of soliton-type solutions in the case of an elastic ice sheet.

3.3.1. Rectangular obstruction

In the following, we focus on analysing perturbed-type supercritical channel flows. Figure 12 demonstrates the influence of the height of the rectangular obstruction on the ice/liquid interface at the Froude number, $F = 1.3$ . The nonlinear results are shown by the solid lines, and the linear approximations are shown by the dot-dashed lines. The maximum deformation occurs directly above the rectangle, indicating strong localised interaction between the elastic interface and the rectangle. In this flow regime, the interface deforms upward above the rectangle in contrast to the downward deflection observed for the Froude numbers $F\lt 1$ in figures 3 and 8.

Figure 11.

Soliton-type (solid lines) and perturbed-type (dashed lines) free-surface profiles without an ice sheet for Froude number $F = 1.3$ , with a rectangular obstruction of height $a = 0.05$ and two widths: $b = 9.0$ (blue) and $b = 5.0$ (red).

Figure 12.

Effect of the rectangular obstruction on the ice–water interface for Froude number $F = 1.3$ , ice thickness $\delta = 0.02$ and rectangle length $b = 16$ . The nonlinear results (solid lines) and the linear approximations (dot-dashed lines) are shown for various obstruction heights: $a = 0.01$ (red line), $0.03$ (blue line) and $0.05$ (black line).

As the height of the rectangle increases, both the amplitude and the wavelength of the interface wave above the rectangle also increase. For example, four wave crests are clearly observed when the height is $a = 0.01$ (red line), and the nonlinear solution closely matches the linear approximation (red solid and dot-dashed lines, respectively). For $a = 0.03$ , a noticeable discrepancy emerges between the nonlinear and linear interfaces (blue solid and dot-dashed lines). In particular, the nonlinear interface exhibits a slightly longer wavelength than the linear approximation predicts, leading to a gradual phase lag between the linear and the nonlinear waves. This phase lag accumulates over the rectangle and then remains constant in the downstream wavy region. At $a = 0.05$ , which represents the maximum height considered, the nonlinear wavelength increase becomes so pronounced that it reduces the number of wave crests above the rectangle. The resulting phase shift between the nonlinear and linear interfaces approaches a full wave cycle and remains constant downstream. The limitation on the rectangle height arises from the significant reduction in flow velocity along the interface above the rectangle, which causes a local decrease in the Froude number rendering the flow subcritical in that region.

Figure 13.

Effect of the rectangular width on the ice–water interface for Froude number $F = 1.3$ and ice thickness $\delta = 0.02$ . The values of the rectangle width are indicated in the figure.

Figure 13 illustrates the FG wave propagating downstream of the obstruction, and its amplitude depends on both the height and width of the rectangle. Although the direction of deformation is opposite to that observed in the channel-subcritical flows, the obstruction width influences the interface similarly. At certain widths, where the rectangle accommodates an integer number of wavelengths, the interface becomes symmetric with respect to the rectangle’s midpoint, and the wave amplitude downstream vanishes. In this case, the interface supports a trapped wave localised above the rectangle. For other widths, elastic waves propagate downstream of the obstruction. Thus, increasing the width leads to an oscillatory behaviour in the downstream wave amplitude, indicating a strong dependence of the interface pattern on both the height and width of the rectangular obstruction.

3.3.2. Sawtooth patch on the bottom

Figure 14 shows the impact of the sawtooth patch on the bottom on the interface shape, demonstrating resonance-like behaviour when the bottom period $b$ approaches the natural wavelength of the interface wave, $\lambda _{ice} = 4.25$ . This behaviour is particularly evident for bottom wavelengths $b = 4.0$ and $b = 4.5$ , which are closest to $\lambda _{ice}$ ; in these cases, the interface exhibits amplified oscillations as typically occurs near the resonance.

Figure 14.

Effect of the undulating bottom on the ice/liquid interface for Froude number $F = 1.3$ and ice thickness $\delta = 0.02$ . The ripple periods are indicated in the figure.

For shorter bottom period ( $b \lt 4.0$ ), the interface response resembles that of a single obstruction with a width equal to the total length of the sawtooth patch. In contrast, for the longer bottom periods, such as $b = 6.5$ and $b = 7.5$ , the downstream wave amplitude rapidly decays, and the deformation above the patch approaches its limiting configuration corresponding to the amplitude above an isolated ripple. This suggests that the influence of adjacent ripples decay, and the interface shape is governed predominantly by local interactions.

4. Conclusions

This work presented results on nonlinear steady-state hydroelastic waves generated by fluid flow beneath an ice cover in a channel with complex bottom topography. The mathematical formulation is based on the integral hodograph method within the framework of the theory of complex functions. The boundary-value problem has been transformed into a system of nonlinear equations for the velocity magnitudes along the interface. A robust numerical scheme, based on Newton’s method, was developed to solve this system, for a wide range of bottom geometries and flow regimes.

The paper presents a comprehensive investigation of nonlinear effects in channel flows beneath an elastic ice sheet, where the bottom profile has an arbitrary polygonal shape where the width can be much larger than the channel depth. Specifically, numerical results are presented for two bottom geometries: an elongated rectangular obstruction/trench and a periodic sawtooth-shaped bottom profile. These configurations are studied across all relevant flow regimes from ice-subcritical to channel-supercritical flows.

A key finding is the identification of a constraint on the maximum height of bottom obstructions: the local flow regime above the obstruction must be the same as it is upstream and downstream. It was also shown that the local Froude number defined based on the local magnitude velocity on the interface and the elevation/depression serves as a reliable indicator of the local flow regime.

Another novel result concerns the strong nonlinear effects that emerge as the local Froude number above the obstruction approaches one of the critical values either the ice-critical or channel-critical Froude number. In the subcritical regime, as the local Froude number nears the ice-critical threshold, the interface not only depresses but also develops a wave whose wavelength corresponds to the critical Froude number. In the ice-supercritical/channel-subcritical regime, nonlinear effects manifest not only in increased wave amplitude but also in an elongation of the wavelength due to the local Froude number increase. This results in a cumulative phase shift along the obstruction, leading to a growing discrepancy between the nonlinear solution and the linear approximation, which then remains constant downstream. In the channel-supercritical regime, both the amplitude and wavelength of the interface wave increase significantly along the obstruction, indicating strong nonlinear behaviour.

An important additional finding is the existence of specific obstruction lengths both for the rectangular and sawtooth profiles at which the interface wave becomes trapped above the obstruction. In these cases, the interface remains essentially flat both upstream and downstream, resulting in spatially confined wave structures.

For the case of a rectangular trench, the study reveals a limiting interface shape as the trench depth tends to infinity. Finally, a resonance-like response of the interface is observed when the length of the sawtooth ripples approaches the natural wavelength of the interface. This phenomenon occurs across all flow regimes, from subcritical to channel supercritical.

Funding

This work is supported by the National Natural Science Foundation of China (Nos. 52521009, W2541018, 52192693, 52192690, 52371270), to which the authors are most grateful.

Declaration of interests

The authors report no conflict of interest.

Appendix A. Dispersion relation

We examine a steady, sinusoidal ice/liquid interface of small steepness $\theta _0$ , such that the slope of the interface can be expressed as

(A1)

\begin{equation} \theta (s) = \mathrm{Re}\big[\theta _0 e^{\textit{ikHs}}\big], \end{equation}

where $kH$ is the non-dimensional wavenumber and $s$ is the arc-length coordinate along the interface.

Differentiating the dynamic boundary condition (2.8) with respect to $s$ and using A1), we obtain

(A2)

\begin{equation} v^2 \frac {d \ln v}{{\rm d}s} = -\frac {1}{F^2} \big(1 + E_b(\textit{kH})^4\big)\theta . \end{equation}

In the absence of an ice sheet ( $E_b = 0$ ), (A2) reduces to

(A3)

\begin{equation} v^2 \frac {d \ln v}{{\rm d}s} = -\frac {\theta }{F^2} = -\frac {\textit{kH}}{\tanh (\textit{kH})}\theta , \end{equation}

where we used the classical dispersion relation for free-surface gravity waves in a channel of depth $H$ (Kochin, Kibel & Roze Reference Kochin, Kibel and Roze1964). We assume that the velocity magnitude along the interface follows the same behaviour as it is for the free-surface water waves. Equating the right-hand sides of (A2) and (A3), we derive the dispersion relation for hydroelastic waves

(A4)

\begin{equation} F^2 \frac {\textit{kH}}{\tanh (\textit{kH})} = 1 + E_b(\textit{kH})^4, \end{equation}

which coincides with the classical result for FG waves (Greenhill Reference Greenhill1887).

Appendix B. Complex velocity in the parameter and potential planes

The boundary-value problem for the function of the complex velocity, defined in the first quadrant of the parameter plane, can be formulated as follows: the velocity magnitude is prescribed along the imaginary axis corresponding to the ice–liquid interface, and the argument of the complex velocity is prescribed along the real axis representing the channel bottom

(B1)

\begin{align} v(\eta ) &= \left | \frac {{\rm d}w}{{\rm d}z} \right |, \quad 0 \leqslant \eta \lt \infty , \end{align}

(B2)

\begin{align} \arg \left ( \frac {{\rm d}w}{{\rm d}z} \right ) &= \Im \left [ \ln \left ( \frac {{\rm d}w}{{\rm d}z} \right ) \right ] = \chi (\xi ), \quad 0 \leqslant \xi \lt \infty . \end{align}

By using Chaplygin’s singular point method (Gurevich Reference Gurevich1965, § 5, Chapter 1), the following integral formula can be obtained for solving the mixed boundary-value problem (B1) and (B2) (Semenov & Yoon Reference Semenov and Yoon2009):

(B3)

\begin{align} & \frac {{\rm d}w}{{\rm d}z} = v_\infty \exp \left [ \frac {1}{\pi } \int _0^\infty \frac {{\rm d}\chi }{{\rm d}\xi } \ln \left ( \frac {\zeta + \xi }{\zeta - \xi } \right ) {\rm d}\xi - \frac {1}{\pi } \int _0^\infty \frac {d \ln v}{{\rm d}\eta } \ln \left ( \frac {\zeta - i\eta }{\zeta + i\eta } \right ) {\rm d}\eta + i\chi _\infty \right ],\nonumber \\& \end{align}

where $v_\infty = \lim _{\eta \to \infty } v(\eta )$ and $\chi _\infty = \lim _{\xi \to \infty } \chi (\xi )$ . One can verify that, for $\zeta = \xi$ , (B3) yields the argument $\chi (\xi )$ , while for $\zeta = i\eta$ , it yields the magnitude $v(\eta )$ , thus satisfying the boundary conditions (B1) and (B2). An alternative derivation of this formula is given in Semenov & Cummings (Reference Semenov and Cummings2007).

The argument function $\chi (\xi )$ is related to the velocity direction $\beta _{{b}i}$ along the bottom

(B4)

\begin{equation} \chi (\xi ) = -\beta _{{b}i}, \qquad \frac {{\rm d}\chi }{{\rm d}\xi } \equiv 0, \qquad \xi _{i-1} \lt \xi \leqslant \xi _i, (i=1,\ldots ,n) ,\end{equation}

with step changes at corner points $A_i$ , where $\zeta = a_i$ . These steps are determined as $\varDelta \chi _i = -(\beta _{{b}i} - \beta _{{b}(i-1)}) = -\alpha _i$ . The integral across the step change at points $\zeta = a_i$ can be evaluated as follows:

(B5)

\begin{align} \lim _{\varepsilon \to 0} \int _{a_i - \varepsilon }^{a_i + \varepsilon } \frac {{\rm d}\chi }{{\rm d}\xi } \ln \left ( \frac {\zeta + \xi }{\zeta - \xi } \right ) {\rm d}\xi &= \ln \left ( \frac {\zeta + a_i}{\zeta - a_i} \right ) \varDelta \chi _i \nonumber \\ &= \alpha _i \ln \left ( \frac {a_i - \zeta }{a_i + \zeta } \right ), \quad i = 1, \ldots , n, \end{align}

where the numerator is inverted to ensure that the argument of the logarithm vanishes at $\zeta = \xi \lt a_i$ .

Evaluating the sum over all corner points and setting $v_\infty = v_0 = 1$ and $\chi _\infty = 0$ , we obtain the expression for the complex velocity in the first quadrant of the parameter plane

(B6)

\begin{equation} \frac {{\rm d}w}{{\rm d}z} = \prod _{i=1}^n \left ( \frac {a_i - \zeta }{a_i + \zeta } \right )^{\alpha _i / \pi } \exp \left [ -\frac {i}{\pi } \int _0^\infty \frac {d \ln v}{{\rm d}\eta } \ln \left ( \frac {i\eta - \zeta }{i\eta + \zeta } \right ) {\rm d}\eta \right ]. \end{equation}

Since the expression for the complex potential is simple (2.12a ), it is convenient to replace the parameter $\zeta$ with the complex potential variable $w$ . We introduce the transformations

(B7)

\begin{equation} \zeta = e^{w'}, \quad w' = \phi ' + i\psi ', \quad \text{where} \quad \phi ' = \frac {\pi \phi }{2}, \quad \psi ' = \frac {\pi \psi }{2}, \end{equation}

and determine $\xi$ and $\eta$ from $\zeta = e^{w'}$

(B8)

\begin{equation} \eta = e^{\phi '}, \quad \psi ' = \frac {\pi }{2}, \quad \xi = e^{\phi '}, \quad \psi ' = 0. \end{equation}

Substituting (B7) into (B6), we obtain the expression for the complex velocity as a function of the complex potential $w$

(B9)

\begin{equation} \frac {{\rm d}w}{{\rm d}z} = \prod _{i=1}^n \left ( \frac {a_i - e^{w'}}{a_i + e^{w'}} \right )^{\alpha _i / \pi } \exp \left [ -\frac {i}{\pi } \int _{-\infty }^{\infty } \frac {d \ln v}{{\rm d}\phi '} \ln \left ( \frac {ie^{\phi '} - e^{w'}}{ie^{\phi '} + e^{w'}} \right ) {\rm d}\phi ' \right ], \end{equation}

where $w' = \pi w/2$ and $\phi ' = \pi \phi /2$ .

Appendix C. Numerical determination of the function $\boldsymbol{v} (\boldsymbol{\phi})$

In discrete form, the solution is sought at fixed points $-\phi ^* \lt \phi _j \lt \phi ^*$ , $j = 1,\ldots ,N$ , located on the upper side of the strip in the $w$ -plane (see figure 1 $c$ ). Then, the points $\eta _j$ on the imaginary axis in the $\zeta $ -plane are obtained $ \eta _j = \exp (\pi \phi _j/2)$ (B8). The integral appearing in (2.13) is evaluated over the intervals ( $\eta _{j-1}, \eta _j$ ) analytically using linear interpolation of the function $\ln v(\eta )$

(C1)

\begin{equation} \frac {1}{\pi } \int _{\eta _{j-1}}^{\eta _j} \frac {d \ln v}{{\rm d}\eta } \ln \left ( \frac {i\eta - \zeta }{i\eta + \zeta } \right ) {\rm d}\eta = \varDelta (\ln v)_j \, c_j(\zeta ), \end{equation}

where

(C2)

\begin{align} \varDelta (\ln v)_j &= \ln v_j - \ln v_{j-1} = \ln \left ( \frac {v_j}{v_{j-1}} \right ), \nonumber \\ c_j(\zeta ) &= \frac {1}{\pi \Delta \eta _j} \int _{\eta _{j-1}}^{\eta _j} \ln \left ( \frac {i\eta - \zeta }{i\eta + \zeta } \right ) {\rm d}\eta . \end{align}

This integral can be analytically evaluated. Substituting (B8) into (C2) enables evaluation of the integral in (2.15).

The curvature and its derivatives are determined using spline interpolation of the nodes $y_m = y(\bar {\phi }_m)$ and $s_m = s(\bar {\phi }_m)$ , which are obtained from (2.18) and (2.19), respectively. The nodes $\bar {\phi }_m = \phi _{4m}$ , $m = 1,\ldots ,M$ , with $M = N/4$ , are selected as every fourth point from the discretised set $\phi _j$ , $j = 1,\ldots ,N$ .

A fifth-order spline is used to ensure continuity of the curvature and its derivatives up to the fourth order, as required by the pressure coefficient in (2.21)

(C3)

\begin{align} & y(s) = y_m + a_{1,m}(s - s_{m-1}) + \cdots + a_{n,m}(s - s_{m-1})^n, \quad s_{m-1} \lt s \lt s_m, \quad n = 5,\nonumber \\& \end{align}

with $m = 1,\ldots ,M$ . The curvature and its derivatives are obtained by differentiating (C3).

The dynamic boundary condition (2.8), expressed in terms of pressure coefficients (2.20) and (2.21), is applied at the nodes $\bar {\phi }_m$ , leading to a system of nonlinear equations

(C4)

\begin{equation} G_m(\bar {V}) = c_{\textit{pm}}(\bar {V}) - c_{\textit{pm}}^{\textit{ice}}(\bar {V}) = 0, \quad m = 1,\ldots ,M, \end{equation}

where $\bar {V} = (v_1, \ldots , v_M)^T$ is the vector of unknown nodal velocities $v_m = v(\bar {\phi }_m)$ . The pressure coefficient due to the flow is

(C5)

\begin{equation} c_{\textit{pm}}(\bar {V}) = 1 - v_m^2 - \frac {2(y_m(\bar {V}) - 1)}{F^2} - C_d v_m^2 \left ( \frac {{\rm d}v}{{\rm d}s} \right )_m, \end{equation}

and the pressure coefficient due to ice bending is

(C6)

\begin{equation} c_{\textit{pm}}^{\textit{ice}}(\bar {V}) = 2D \left [ \left ( \frac {d^2 \kappa }{{\rm d}s^2} \right )_m + \frac {1}{2} \kappa _m^3 \right ]. \end{equation}

The nonlinear system (C4) is solved using Newton’s method. The Jacobian matrix is computed numerically via central differences with step size $\varDelta v_m = 10^{-8}$ . Each function evaluation $G_m(\bar {V})$ involves solving the nonlinear system (2.16) using an internal Newton iteration. Intermediate velocities $v_j$ , for $j = 4m-3, 4m-2, 4m-1$ , are computed via linear interpolation of $\ln v_m$ from the solution of system (C4). Convergence is achieved when the residual satisfies

(C7)

\begin{equation} \frac {1}{M} \sum _{m=1}^{M} |G_m(\bar {V}^*)| \lt 10^{-6}, \end{equation}

which is regarded as sufficiently accurate for the nonlinear system.

For supercritical flow regimes, $F\gt F_{{cr}}$ , wave attenuation at the ends of the calculation region is applied. The damping intervals are chosen as $x_{\!P_2} - x_{\!P_1} = 2\lambda _{{gr}}$ and $x_{T_2} - x_{T_1} = 4\lambda _{{gr}}$ , where $\lambda _{{gr}}$ is the gravity wave wavelength determined from the dispersion relation. The damping coefficients $C_{{up}}$ and $C_{{dw}}$ are selected from the interval $0.2\lambda _{{gr}}$ to $0.4\lambda _{{gr}}$ to ensure effective dissipation of both elastic and gravity waves.

In the first iteration, the initial values are taken as $v^{(1)} \equiv 1$ , $s_b(\phi ) = \phi$ and $\phi _{A_i} = s_{bi}$ . The next iteration begins with solving the nonlinear system (2.16).

The number of nodes along the interface is chosen in the range $800 \lt N \lt 4000$ , depending on the required resolution. Typically, 12 nodes per shortest elastic wavelength are used to ensure accuracy. For $N = 800$ , the computation takes a few minutes; for $N = 4000$ , the run time is approximately 15 min on a Precision Tower T7920 workstation.

References

Alam, M.-R., Liu, Y. & Yue, D.K.P. 2009 Bragg resonance of waves in a two-layer fluid propagating over bottom ripples. Part I. Perturbation analysis. J. Fluid Mech. 624, 191–224.10.1017/S0022112008005478CrossRef Google Scholar

Barman, K.K., Chanda, A., Tsai, C.-C. & Hsu, T.-W. 2024 b Bragg scattering of gravity waves by a sea bed of varying depth in the presence of uniform current covered by a floating membrane. Phys. Fluids 36, 012118.10.1063/5.0183629CrossRef Google Scholar

Barman, K.K., Chanda, A., Tsai, C.-C. & Mondale, S. 2024 a Flexural-gravity wave interaction with undulating bottom topography in the presence of uniform current: an asymptotic approach. Appl. Math. Model. 133, 41–64.10.1016/j.apm.2024.05.012CrossRef Google Scholar

Blyth, M.G., Părău, E.I. & Vanden-Broeck, J.-M. 2011 Hydroelastic waves on fluid sheets. J. Fluid Mech. 689, 541–551.10.1017/jfm.2011.451CrossRef Google Scholar

Bonnefoy, F., Meylan, M.H. & Ferrant, P. 2009 Nonlinear wave generation and absorption in a 2D hydrodynamic flume. J. Fluid Mech. 621, 215–242.10.1017/S0022112008004849CrossRef Google Scholar

Boral, S., Meylan, M.H., Sahoo, T. & Ni, B.-Y. 2023 Time-dependent flexural gravity wave scattering due to uneven bottomin the paradigm of blocking dynamics. Phys. Fluids 35, 116603.10.1063/5.0173542CrossRef Google Scholar

Chamberlain, P.G. & Porter, D. 1995 The modified mild-slope equation. J. Fluid Mech. 291, 393–407.10.1017/S0022112095002758CrossRef Google Scholar

Davies, A.G. 1982 The reflection of wave energy by undulations on the seabed. Dyn. Atmos. Oceans 6 (4), 207–232.10.1016/0377-0265(82)90029-XCrossRef Google Scholar

Dias, F. & Vanden-Broeck, J.-M. 1989 Open channel flows with submerged obstructions. J. Fluid Mech. 206, 155–170.10.1017/S0022112089002260CrossRef Google Scholar

Forbes, L.K. 1986 Surface waves of large amplitude beneath an elastic sheet. Part 1. High-order series solution. J. Fluid Mech. 169, 409–428.10.1017/S0022112086000708CrossRef Google Scholar

Forbes, L.K. 1988 Surface waves of large amplitude beneath an elastic sheet. Part 2. Galerkin solution. J. Fluid Mech. 188, 491–508.10.1017/S0022112088000813CrossRef Google Scholar

Gao, W., Wang, Z. & Milewski, P.A. 2019 Nonlinear hydroelastic waves on a linear shear current at finite depth. J. Fluid Mech. 876, 55–86.10.1017/jfm.2019.528CrossRef Google Scholar

Greenhill, A.G. 1887 Wave motion in hydrodynamics. Am. J. Math. 9, 62–112.10.2307/2369499CrossRef Google Scholar

Gurevich, M.I. 1965 Theory of Jets in Ideal Fluids. Academic Press.Google Scholar

Guyenne, P. & Părău, E.I. 2012 Computations of fully nonlinear hydroelastic solitary waves on deep water. J. Fluid Mech. 713, 307–329.10.1017/jfm.2012.458CrossRef Google Scholar

Joukowskii, N.E. 1890 Modification of Kirchhoff’s method for determination of a fluid motion in two directions at a fixed velocity given on the unknown streamline. Math. Sbornik. 15 (1), 121–278.Google Scholar

Karmakar, D., Bhattacharjee, J. & Sahoo, T. 2010 Oblique flexural gravity-wave scattering due to changes in bottom topography. J. Engng Maths 66, 325–341.10.1007/s10665-009-9297-8CrossRef Google Scholar

Kheisin, D.E. 1963 Moving load on an elastic plate which floats on the surface of an ideal fluid (in Russian). Izv. Akad. Nauk SSSR Otd. Tekh. Nauk. Mekh. i Mashinostroenie 1, 178–180.Google Scholar

Kheisin, D.E. 1967 Dynamics of ice cover. Tech. Rep. FSTC-HT-23-485-69. U.S. Army Foreign Science and Technology Center.Google Scholar

Kochin, N.E., Kibel, I.A. & Roze, N.V. 1964 Theoretical Hydromechanics, vol. 1964. Wiley Interscience.Google Scholar

Korobkin, A.A., Malenica, S. & Khabakhpasheva, T.I. 2018 Interaction of flexural-gravity waves in ice cover with vertical walls. Phil. Trans. R. Soc. A Math. Phys. Engng Sci. 376 (2129), 20170347.10.1098/rsta.2017.0347CrossRef Google Scholar PubMed

Korobkin, A.A., Părău, E.I. & Vanden-Broeck, J.-M. 2011 The mathematical challenges and modelling of hydroelasticity. Phil. Trans. R. Soc. Lond. A 369, 2803–2812.Google Scholar PubMed

Mei, C.C. 1985 Resonant reflection of surface water waves by periodic sandbars. J. Fluid Mech. 152, 315–335.10.1017/S0022112085000714CrossRef Google Scholar

Michell, J.H. 1890 On the theory of free streamlines. Phil. Trans. R. Soc. Lond. A 181, 389–431 Google Scholar

Milewski, P.A., Vanden-Broeck, J.-M. & Wang, Z. 2011 Hydroelastic solitary waves in deep water. J. Fluid Mech. 679, 628–640.10.1017/jfm.2011.163CrossRef Google Scholar

Milewski, P.A. & Wang, Z. 2013 Three dimensional flexural–gravity waves. Stud. Appl. Maths 131, 135–148.10.1111/sapm.12005CrossRef Google Scholar

Ni, B.-Y., Han, D.-F., Di, S.-C. & Xue, Y.-Z. 2020 On the development of ice-water-structure interaction. J. Hydrodyn. 32, 629–652.10.1007/s42241-020-0047-8CrossRef Google Scholar

Ni, B.-Y., Khabakhpasheva, T.I. & Semenov, Y.A. 2023 Nonlinear gravity waves in the channel covered by broken ice. Phys. Fluids 35, 102118.10.1063/5.0166567CrossRef Google Scholar

Ni, B.-Y., Semenov, Y.A., Khabakhpasheva, T.I., Părău, E.I. & Korobkin, A.A. 2024 Nonlinear ice sheet/liquid interaction in a channel with an obstruction. J. Fluid Mech. 983, A41.10.1017/jfm.2024.177CrossRef Google Scholar

Page, C. & Părău, E.I. 2014 Hydraulic falls under a floating ice plate due to submerged obstructions. J. Fluid Mech. 745, 208–222.10.1017/jfm.2014.92CrossRef Google Scholar

Părău, E.I. & Dias, F. 2002 Nonlinear effects in the response of a floating ice plate to a moving load. J. Fluid Mech. 460, 281–305.10.1017/S0022112002008236CrossRef Google Scholar

Plotnikov, P.I. & Toland, J.F. 2011 Modelling nonlinear hydroelastic waves. Phil. Trans. R. Soc. Lond. A 369, 2942–2956.Google Scholar PubMed

Porter, D. & Porter, R. 2004 Approximations to wave scattering by an ice sheet of variable thickness over undulating bed topography. J. Fluid Mech. 509, 145–179.10.1017/S0022112004009267CrossRef Google Scholar

Semenov, Y.A. 2021 Nonlinear flexural-gravity waves due to a body submerged in the uniform stream. Phys. Fluids 33 (5), 052115.10.1063/5.0048887CrossRef Google Scholar

Semenov, Y.A. & Cummings, L.J. 2007 Free boundary Darcy flows with surface tension: analytical and numerical study. Eur. J. Appl. Maths 17, 607–631.10.1017/S0956792506006759CrossRef Google Scholar

Semenov, Y.A. & Yoon, B.-S. 2009 Onset of flow separation at oblique water impact of a wedge. Phys. Fluids 21, 112103–112111.10.1063/1.3261805CrossRef Google Scholar

Shishmarev, K., Khabakhpasheva, T.I. & Korobkin, A.A. 2016 The response of ice cover to a load moving along a frozen channel. Appl. Ocean Res. 59, 313–326.10.1016/j.apor.2016.06.008CrossRef Google Scholar

Shishmarev, K., Khabakhpasheva, T.I. & Korobkin, A.A. 2019 Ice response to an underwater body moving in a frozen channel. Appl. Ocean Res. 91, 101877.10.1016/j.apor.2019.101877CrossRef Google Scholar

Squire, V.A. 2020 Ocean wave interactions with sea ice: a reappraisal. Annu. Rev. Fluid Mech. 52, 37–60.10.1146/annurev-fluid-010719-060301CrossRef Google Scholar

Squire, V.A., Dugan, J.P., Wadhams, P., Rottier, P.J. & Liu, A.K. 1995 Of ocean waves and sea ice. Annu. Rev. Fluid Mech. 27, 115–168.10.1146/annurev.fl.27.010195.000555CrossRef Google Scholar

Sturova, I.V. 2009 Time-dependent response of a heterogeneous elastic plate floating on shallow water of variable depth. J. Fluid Mech. 637, 305–325.10.1017/S0022112009990504CrossRef Google Scholar

Xiong, H., Ni, B., Semenov, Y. & Korobkin, A. 2025 Inertial motion of a circular cylinder approaching obliquely an ice cover. Sci Rep. 15, 9412.10.1038/s41598-025-93435-1CrossRef Google Scholar PubMed

Xue, Y.Z., Zeng, L.D., Ni, B.Y., Korobkin, A.A. & Khabakhpasheva, T. 2021 Hydroelastic waves in an ice-covered channel with leads. Phys. Fluids 33, 037109.10.1063/5.0037682CrossRef Google Scholar

Yoon, B.-S. & Semenov, Y.A. 2011 Separated inviscid sheet flows. J. Fluid Mech. 678, 511–534.10.1017/jfm.2011.123CrossRef Google Scholar

Figure 1.($a$) Physical plane, ($b$) parameter, or $\zeta$-plane, and compex potential, or $w$-plane.

Figure 2.Wavenumber vs Froude number for different thicknesses of the ice sheet, $\delta$.

Figure 3.Nonlinear (solid lines) and linear (dash-dotted lines) approximations of the ice/liquid interface, and the local Froude number above the rectangle (dashed lines), for various rectangle heights shown by the solid thick line at the bottom. (a) Froude number $ F = 0.5$ and obstruction heights $ a = 0.04$ (black), $ 0.06$ (red), $ 0.08$ (blue); (b) Froude number $ F = 0.4$ and obstruction heights $ a = 0.08$ (black), $ 0.15$ (red), $ 0.21$ (blue). The rectangle width is $ b = 12$, and the ice thickness is $ \delta = 0.003$.

Figure 4.Effect of rectangle width on the ice/liquid interface at the onset of convergence of the solution for Froude number $F = 0.5$ and ice thickness $\delta = 0.003$. (a) Interface profiles for varying rectangle widths and corresponding heights: $b = 12$, $a = 0.085$ (black); $b = 8$, $a = 0.085$ (red); $b = 4$, $a = 0.09$ (blue); $b = 2$, $a = 0.140$ (olive); $b = 1.5$, $a = 0.095$ (wine); $b = 1.0$, $a = 0.10$ (cyan); $b = 0.5$, $a = 0.130$ (magenta); $b = 0.1$, $a = 0.205$ (orange). (b) Close-up view: $b = 4.0$, $a = 0.09$ (blue); $b = 3.0$, $a = 0.09$ (pink); $b = 2.3$, $a = 0.095$ (cyan); $b = 2.1$, $a = 0.135$ (orange); $b = 2.0$, $a = 0.140$ (olive); $b = 1.9$, $a = 0.140$ (magenta); $b = 1.8$, $a = 0.095$ (green).

Figure 5.(a) Effect of trench depth on the interface shape for a trench width $b = 12$, Froude number $F = 0.5$ and ice thickness $\delta = 0.003$. Results are shown for varying trench depths $a$: 0.2 (orange), 0.5 (blue), 1.0 (olive), 2.0 (magenta), 4.0 (red) and 10 (black). (b) Maximum deformation at infinite trench depth for different trench widths: $b = 12$ (black), 6 (red), 4 (magenta), 2 (olive) and 1 (blue).

Figure 6.(a) Effect of the sawtooth bottom profile on the interface shape for Froude number $F = 0.5$, ice thickness $\delta = 0.003$ and two different heights of the ripples: (a) $a = 0.125$, (b) $a = 0.115$. The ripple base is indicated in the figure.

Figure 7.(a) Ice–water interface and (b) bending moment, $m=-D\kappa$, along the interface for ice thickness $\delta = 0.003$ and various Froude numbers in the range $F_{{cr}} \lt F \lt 1$, as indicated in the figure. The rectangular obstruction has height $a = 0.03$ and width $b = 16$; the critical Froude number is $F_{{cr}} = 0.5748$.

Figure 8.Effect of rectangle width $ b$ on the interface shape for Froude number $ F = 0.7$ and rectangle height $ a = 0.03$. (a) Nonlinear results for widths $ b_{{ab}} = 17.9$ (blue), $ b_{\textit{ag}} = 15.58$ (green) and $ b_{{ar}} = 13.5$ (red); (b), (c) and (d) show comparisons between the nonlinear results (solid lines) and the corresponding linear approximations for widths $ b_{{bb}} = 13.5$ (blue), $ b_{{bg}} = 11.68$ (green) and $ b_{{br}} = 9.5$ (red), respectively.

Figure 9.Effect of rectangle width $ b$ on the free-surface shape for Froude number $ F = 0.5$ and rectangle height $ a = 0.1$. (a) Nonlinear results for widths $ b_{{ab}} = 9.4$ (red), $ b_{\textit{ag}} = 8.8$ (green) and $ b_{{ar}} = 8.4$ (blue); (b) shows comparisons between the nonlinear results from ($a$) (solid lines of the same colour) and the corresponding linear approximations (shot-dashed lines).

Figure 10.Ice/liquid interface for Froude number $F = 0.7$, ice thickness $\delta = 0.003$ and ripple height $a = 0.015$, shown for various ripple periods $b$ indicated in the figure. Resonance is observed at $b = 3.25$, which is close to the FG wavelength $\lambda _{{gr}} = 3.12$.

Figure 11.Soliton-type (solid lines) and perturbed-type (dashed lines) free-surface profiles without an ice sheet for Froude number $F = 1.3$, with a rectangular obstruction of height $a = 0.05$ and two widths: $b = 9.0$ (blue) and $b = 5.0$ (red).

Figure 12.Effect of the rectangular obstruction on the ice–water interface for Froude number $F = 1.3$, ice thickness $\delta = 0.02$ and rectangle length $b = 16$. The nonlinear results (solid lines) and the linear approximations (dot-dashed lines) are shown for various obstruction heights: $a = 0.01$ (red line), $0.03$ (blue line) and $0.05$ (black line).

Figure 13.Effect of the rectangular width on the ice–water interface for Froude number $F = 1.3$ and ice thickness $\delta = 0.02$. The values of the rectangle width are indicated in the figure.

Figure 14.Effect of the undulating bottom on the ice/liquid interface for Froude number $F = 1.3$ and ice thickness $\delta = 0.02$. The ripple periods are indicated in the figure.

Article contents

Nonlinear flexural-gravity waves for flows over bottom topography

Abstract

JFM classification

Information

1. Introduction

2. Formulation of the problem and analysis

2.1. Complex potential

2.2. System of equations for the unknowns $\phi _{A_i}$ and the function $v(\eta )$

2.3. Analysis of the dispersion relation

3. Results and discussion

3.1. Subcritical flows

3.1.1. Rectangular obstruction

3.1.2. Sawtooth patch on the bottom

3.2. Ice-supercritical–channel-subcritical flows, $F_{{cr}} \lt F \lt 1$

3.2.1. Rectangular obstruction

3.2.2. Sawtooth patch on the bottom

3.3. Channel-supercritical flows, $F \gt 1$

3.3.1. Rectangular obstruction

3.3.2. Sawtooth patch on the bottom

4. Conclusions

Funding

Declaration of interests

Appendix A. Dispersion relation

Appendix B. Complex velocity in the parameter and potential planes

Appendix C. Numerical determination of the function $\boldsymbol{v} (\boldsymbol{\phi})$

References

Save article to Kindle

Save article to Dropbox

Save article to Google Drive

Reply to: Submit a response

Your details

You have entered the maximum number of contributors

Conflicting interests